Process airborne PFP merge files
Contents
Process airborne PFP merge files¶
R Program to read in HIPPO, ORCAS, and ATom PFP merge files, massage time variables, add a strat flag, and subset
then filter aircraft data for strong local continental influences, subtract off NOAA in situ SPO, and write out flat text files
have to run process_aircraft_10s.ipynb first
library('ncdf4')
library('yaml')
project_tmpdir_obs = read_yaml('../_config_calc.yml')$project_tmpdir_obs
username = Sys.info()['user']
project_tmpdir_obs = gsub('\\{\\{env\\[\'USER\'\\]\\}\\}', username, project_tmpdir_obs)
# specify aircraft data file names
hippomergefile=paste(project_tmpdir_obs,'/aircraft-merge-products/HIPPO_noaa_flask_allparams_merge_insitu_20121129.tbl',sep='') # this is the official HIPPO merge product
atommergedir=paste(project_tmpdir_obs,'/aircraft-merge-products',sep='') # these are version 2.0 (21-08-26)
atommergefiles=c('MER-PFP_DC8_ATom-2.nc','MER-PFP_DC8_ATom-3.nc','MER-PFP_DC8_ATom-4.nc') # no ATom-1 PFP data
# no ORCAS PFP data
# set strat flag cutoffs for use below based upon Jin et al., 2021 (https://doi.org/10.5194/acp-21-217-2021)
stratcoh2o=50
stratcoo3=150
stratcon2o=319
# read in global N2O for detrending aircraft N2O
glbn2ofile=url('ftp://aftp.cmdl.noaa.gov/products/trends/n2o/n2o_annmean_gl.txt')
hlines=61
glbn2o=read.table(glbn2ofile,skip=hlines,header=F,stringsAsFactors=F)
colnames(glbn2o)=c('year','n2o','unc')
# for HIPPO, need to interpolate prof variable from 10-sec merge - load from process_aircraft_10s.ipynb output before reading PFP data below
load('HIPPO_10s.RData')
hippomerge10s=hippomerge # 'hippomerge' reused below
hippo10sdt=ISOdatetime(hippomerge$year,hippomerge$mon,hippomerge$day,hippomerge$hour,hippomerge$min,hippomerge$sec,tz='UTC')
read in HIPPO file, calc strat flag, and subset¶
# read and add time variables
hippomerge=read.table(hippomergefile,header=T)
hippodt=strptime(paste(hippomerge[,"Year"],hippomerge[,"DOY"]),format='%Y %j',tz='UTC')+hippomerge[,"UTC"] # DOY is day of year of takeoff; UTC is seconds since midnight on day of takeoff
hippomerge$Month=as.POSIXlt(hippodt)$mon+1
hippomerge$Day=as.POSIXlt(hippodt)$mday
hippomerge$Hour=as.POSIXlt(hippodt)$hour
hippomerge$Min=as.POSIXlt(hippodt)$min
hippomerge$Sec=as.POSIXlt(hippodt)$sec
# interpolate prof from 10 sec file
hippomerge$prof=approx(as.POSIXct(hippo10sdt),hippomerge10s$prof,as.POSIXct(hippodt),method='constant',f=0)$y
# add strat flag
hippomerge$strat=rep(0,nrow(hippomerge)) # 0 means trop
h2oref=hippomerge$H2Oppmv_vxl; h2oref[is.na(h2oref)]=hippomerge$H2O_UWV[is.na(h2oref)]
hippomerge$h2oref=h2oref # for output
h2oref[is.na(h2oref)]=0 # if H2O missing treat as if potentially strat
n2oref=hippomerge$N2O_CCG; n2oref[is.na(n2oref)]=hippomerge$N2O_QCLS[is.na(n2oref)]
n2oref=n2oref-(approx(glbn2o$year+0.5,glbn2o$n2o,hippomerge$Year+hippomerge$DOY/365)$y-glbn2o$n2o[glbn2o$year==2009])
hippomerge$n2oref=n2oref # for output
n2oref[is.na(n2oref)]=400 # if N2O missing do not use for filter
o3ref=hippomerge$O3_ppb; o3ref[is.na(o3ref)]=hippomerge$O3_UO3[is.na(o3ref)]
hippomerge$o3ref=o3ref # for output
o3ref[is.na(o3ref)]=0 # if O3 missing do not use for filter
hippomerge$strat[h2oref<stratcoh2o&(o3ref>stratcoo3|n2oref<stratcon2o|(o3ref==0&n2oref==400))]=1 # if either o3 or n2o criteria are met, or if both are missing, consider strat
hippomerge$strat[h2oref==0&o3ref==0&n2oref==400&hippomerge$GGALT<8000]=0 # if all 3 missing assume < 8 km is trop
# select columns
colsel=c('Year','Month','Day','Hour','Min','Sec','H.no','flt','prof','GGLAT','GGLON','GGALT','PSXC','THETA','CO2_CCG','CO2_QCLS','CO2_OMS','CO2_AO2','CH4_CCG','CH4_QCLS','SF6_CCG','strat','h2oref','n2oref','o3ref')
hippomerge=hippomerge[,is.element(colnames(hippomerge),colsel)]
hippomerge=hippomerge[,match(colsel,names(hippomerge))] # reorder
names(hippomerge)=c('year','mon','day','hour','min','sec','camp','flt','prof','lat','lon','alt','pressure','theta','co2','co2qcls','co2oms','co2ao2','ch4pfp','ch4qcls','sf6pfp','strat','h2oref','n2oref','o3ref') ## 'co2' = CO2_CCG = co2pfp
read in ATom files, calc strat flag, and subset¶
# read and add time variables
atomvar=c('time','Flight_Date','DLH-H2O/H2O_DLH','UCATS-H2O/H2O_UWV','UCATS-O3/O3_UCATS','QCLS-CH4-CO-N2O/N2O_QCLS','NOyO3-O3/O3_CL','MMS/G_ALT','RF','prof.no','MMS/P','MMS/POT','MMS/G_LAT','MMS/G_LONG',
'MMS/G_ALT','NOAA-Picarro/CO2_NOAA','QCLS-CO2/CO2_QCLS','AO2/CO2_AO2','CO2.X','NOAA-Picarro/CH4_NOAA','QCLS-CH4-CO-N2O/CH4_QCLS','PFP/CO2_PFP','PFP/CH4_PFP','PFP/N2O_PFP','PFP/SF6_CCGG_PFP')
atommerge=NULL
for(i in c(1:3)){
atomnc=nc_open(paste(atommergedir,'/',atommergefiles[i],sep=''))
count=length(ncvar_get(atomnc,'time'))
campdata=NULL
for(var in atomvar){
campdata=cbind(campdata,ncvar_get(atomnc,var))
}
campdata=cbind(campdata,rep(i+1,count)) # A.no
nc_close(atomnc)
atommerge=rbind(atommerge,campdata)
}
atommerge=data.frame(atommerge,stringsAsFactors=F)
names(atommerge)=c(gsub('.*/','',atomvar),'A.no')
atommerge$YYYYMMDD=atommerge$Flight_Date
atomdt=as.POSIXlt(ISOdatetime(2016,1,1,0,0,0,tz='UTC')+atommerge$time,tz='UTC')
atommerge$Year=atomdt$year+1900
atommerge$Month=as.POSIXlt(atomdt)$mon+1
atommerge$Day=as.POSIXlt(atomdt)$mday
atommerge$Hour=as.POSIXlt(atomdt)$hour
atommerge$Min=as.POSIXlt(atomdt)$min
atommerge$Sec=as.POSIXlt(atomdt)$sec
# add strat flag
atommerge$strat=rep(0,nrow(atommerge)) # 0 means trop
h2oref=atommerge$H2O_DLH; h2oref[is.na(h2oref)]=atommerge$H2O_UWV[is.na(h2oref)]
atommerge$h2oref=h2oref # for output
h2oref[is.na(h2oref)]=0 # if H2O missing treat as if potentially strat
n2oref=atommerge$N2O_QCLS; n2oref[is.na(n2oref)]=atommerge$N2O_PFP[is.na(n2oref)]
n2oref=n2oref-(approx(glbn2o$year+0.5,glbn2o$n2o,atommerge$Year+atomdt$yday/365)$y-glbn2o$n2o[glbn2o$year==2009])
atommerge$n2oref=n2oref # for outptut
n2oref[is.na(n2oref)]=400 # if N2O missing do not use for filter
o3ref=atommerge$O3_CL; o3ref[is.na(o3ref)]=atommerge$O3_UCATS[is.na(o3ref)]
atommerge$o3ref=o3ref
o3ref[is.na(o3ref)]=0 # if O3 missing do not use for filter
atommerge$strat[h2oref<stratcoh2o&(o3ref>stratcoo3|n2oref<stratcon2o|(o3ref==0&n2oref==400))]=1 # if either o3 or n2o criteria are met, or if both are missing, consider strat
atommerge$strat[h2oref==0&o3ref==0&n2oref==400&atommerge$G_ALT<8000]=0 # if all 3 missing assume < 8 km is trop
# select column
colsel=c('Year','Month','Day','Hour','Min','Sec','A.no','RF','prof.no','G_LAT','G_LONG','G_ALT','P','POT','CO2_PFP','CO2_NOAA','CO2_QCLS','CO2_AO2','CH4_PFP','CH4_NOAA','CH4_QCLS','SF6_CCGG_PFP','strat','h2oref','n2oref','o3ref')
atommerge=atommerge[,is.element(colnames(atommerge),colsel)]
atommerge=atommerge[,match(colsel,names(atommerge))] # reorder
names(atommerge)=c('year','mon','day','hour','min','sec','camp','flt','prof','lat','lon','alt','pressure','theta','co2','co2noaa','co2qcls','co2ao2','ch4pfp','ch4noaa','ch4qcls','sf6pfp','strat','h2oref','n2oref','o3ref') ## 'co2' = CO2_PFP
Filter airborne PFP data and subtract SPO¶
# read in NOAA in situ record from SPO
sponc=nc_open(paste(project_tmpdir_obs,'/obspack_co2_1_GLOBALVIEWplus_v6.0_2020-09-11/data/nc/co2_spo_surface-insitu_1_allvalid.nc',sep=''))
spoco2=data.frame(cbind(ncvar_get(sponc,'time_decimal'),t(ncvar_get(sponc,'time_components')),ncvar_get(sponc,'value')*1E6)) ; colnames(spoco2)=c('date','year','mon','day','hour','min','sec','co2')
qcflag=ncvar_get(sponc,'qcflag'); spoco2$co2[substr(qcflag,1,1)!='.']=NA; spoco2$co2[substr(qcflag,2,2)!='.']=NA
spodt=ISOdatetime(spoco2$year,spoco2$mon,spoco2$day,spoco2$hour,spoco2$min,spoco2$sec,tz='UTC')
HIPPO¶
# filter
ints=read.table(paste(project_tmpdir_obs,'/hippo_xsect_filt_datetime.txt',sep=''),header=T)
startdt=ISOdatetime(ints$startyear,ints$startmon,ints$startday,ints$starthour,ints$startmin,ints$startsec,tz='UTC')
stopdt=ISOdatetime(ints$stopyear,ints$stopmon,ints$stopday,ints$stophour,ints$stopmin,ints$stopsec,tz='UTC')
blfilt=rep(T,nrow(hippomerge))
for(i in c(1:nrow(ints))){
blfilt[difftime(hippodt,startdt[i])>=0&difftime(hippodt,stopdt[i])<=0]=F
}
hippodt=hippodt[blfilt]
hippomerge=hippomerge[blfilt,]
print(paste('Filtered ',sum(!blfilt),' of ',length(blfilt),' HIPPO obs (',round(sum(!blfilt)/length(blfilt)*100,1),'%)',sep=''))
[1] "Filtered 14 of 1374 HIPPO obs (1%)"
# calculate differences
hippomerge$co2mspo=round(hippomerge$co2-approx(as.POSIXct(spodt),spoco2$co2,as.POSIXct(hippodt))$y,3) ## co2 = 'CO2_CCG'
hippomerge$co2mqcls=round(hippomerge$co2-hippomerge$co2qcls,3)
hippomerge$co2moms=round(hippomerge$co2-hippomerge$co2oms,3)
hippomerge$co2mao2=round(hippomerge$co2-hippomerge$co2ao2,3)
hippomerge$ch4mqcls=round(hippomerge$ch4pfp-hippomerge$ch4qcls,3)
# write out
write(names(hippomerge),'../data/aircraft-obs/HIPPO_SO_mSPO_pfp.txt',ncol=ncol(hippomerge))
write(t(hippomerge),'../data/aircraft-obs/HIPPO_SO_mSPO_pfp.txt',ncol=ncol(hippomerge),append=T)
print(apply(!is.na(hippomerge),2,sum))
year mon day hour min sec camp flt
1360 1360 1360 1360 1360 1360 1360 1360
prof lat lon alt pressure theta co2 co2qcls
1360 1360 1360 1360 1355 1355 1296 944
co2oms co2ao2 ch4pfp ch4qcls sf6pfp strat h2oref n2oref
995 1125 1348 1010 1347 1360 1360 1355
o3ref co2mspo co2mqcls co2moms co2mao2 ch4mqcls
1350 1296 911 965 1071 1000
ATom¶
# filter
ints=read.table(paste(project_tmpdir_obs,'/atom_xsect_filt_datetime.txt',sep=''),header=T)
startdt=ISOdatetime(ints$startyear,ints$startmon,ints$startday,ints$starthour,ints$startmin,ints$startsec,tz='UTC')
stopdt=ISOdatetime(ints$stopyear,ints$stopmon,ints$stopday,ints$stophour,ints$stopmin,ints$stopsec,tz='UTC')
blfilt=rep(T,nrow(atommerge))
for(i in c(1:nrow(ints))){
blfilt[difftime(atomdt,startdt[i])>=0&difftime(atomdt,stopdt[i])<=0]=F
}
atomdt=atomdt[blfilt]
atommerge=atommerge[blfilt,]
print(paste('Filtered ',sum(!blfilt),' of ',length(blfilt),' ATom obs (',round(sum(!blfilt)/length(blfilt)*100,1),'%)',sep=''))
[1] "Filtered 13 of 875 ATom obs (1.5%)"
# calculate differences
atommerge$co2mspo=round(atommerge$co2-approx(as.POSIXct(spodt),spoco2$co2,as.POSIXct(atomdt))$y,2) ## co2 = 'CO2_PFP'
atommerge$co2mqcls=round(atommerge$co2-atommerge$co2qcls,3)
atommerge$co2mao2=round(atommerge$co2-atommerge$co2ao2,3)
atommerge$co2mnoaa=round(atommerge$co2-atommerge$co2noaa,3)
atommerge$ch4mqcls=round(atommerge$ch4pfp-atommerge$ch4qcls,3)
atommerge$ch4mnoaa=round(atommerge$ch4pfp-atommerge$ch4noaa,3)
# write out
write(names(atommerge),'../data/aircraft-obs/ATOM_SO_mSPO_pfp.txt',ncol=ncol(atommerge))
write(t(atommerge),'../data/aircraft-obs/ATOM_SO_mSPO_pfp.txt',ncol=ncol(atommerge),append=T)
print(apply(!is.na(atommerge),2,sum))
year mon day hour min sec camp flt
862 862 862 862 862 862 862 862
prof lat lon alt pressure theta co2 co2noaa
862 862 862 862 862 862 859 838
co2qcls co2ao2 ch4pfp ch4noaa ch4qcls sf6pfp strat h2oref
592 742 849 838 585 849 862 862
n2oref o3ref co2mspo co2mqcls co2mao2 co2mnoaa ch4mqcls ch4mnoaa
860 862 859 591 739 835 577 826