accomp<-function(start='090601',end='090630',locs=c('NWR','FEF','SPL','HDP','EFS','RBA'),outflag='',ico=0.0,eco=1.0,png=F,nlines=c(3,3,3,3,3,3),c13cor=T){ #start='090601' #end='090627' #locs=c('NWR','FEF','SPL','HDP','EFS','RBA') #outflag='' #ico=0.0 #eco=1.0 #png=F #nlines=c(3,3,3,3,3,3) log<-file(paste("comp_",start,"_",end,outflag,".log",sep=''),open="wt") sink(log) sink(log,type="message") # R can only plot in current time zone Sys.setenv(TZ="GMT") if(outflag!=''){outflag<-paste('_',outflag,sep='')} # start and end cannot span a year change ## CALCULATE DAYS OF YEAR #dpm<-c(31,28,31,30,31,30,31,31,30,31,30,31) #if(as.numeric(substr(start,1,2))/4-trunc(as.numeric(substr(start,1,2))/4)==0) dpm[2]<-29 # can not handle multiple year runs including a Feb. leap year #ico<-ico+c(0,cumsum(dpm)[1:11])[as.numeric(substr(start,3,4))]+as.numeric(substr(start,5,6)) #eco<-eco+c(0,cumsum(dpm)[1:11])[as.numeric(substr(end,3,4))]+as.numeric(substr(end,5,6)) ## add days to account for whole years since Dec. 31, 2004 #year<-as.numeric(substr(start,1,2)) # can not handle multiple year runs #yrdays<-0 #if(year>5){ #for(i in c(5:(year-1))){ #if(i/4-trunc(i/4)==0){ # leap year #yrdays<-yrdays+366 #}else{ #yrdays<-yrdays+365 #} #} #} #ico<-ico+yrdays #eco<-eco+yrdays # Calculate year.frac startdt=strptime(start,format='%y%m%d',tz="GMT") enddt=as.POSIXlt(strptime(end,format='%y%m%d',tz="GMT")+60*60*24) # end of day ye=startdt$year+1900 if(ye/4-trunc(ye/4)==0){ dpy=366 } else { dpy=365 } ico=round(ye+as.numeric(difftime(startdt,strptime(paste(ye,'-01-01 00:00:00',sep=''),format='%Y-%m-%d %H:%M:%S',tz="GMT"),units='days'))/dpy,7) eco=round(ye+as.numeric(difftime(enddt,strptime(paste(ye,'-01-01 00:00:00',sep=''),format='%Y-%m-%d %H:%M:%S',tz="GMT"),units='days'))/dpy,7) # READ IN AND MERGE DATA minyrfrac<-1E36 maxyrfrac<-0 print(c(ico,eco)) for(i in locs){ filename<-'' oldruns<-system(paste('ls -rt ',i,'*.lco',sep=''),T) data<-NULL print(oldruns) if(length(oldruns)>0){ # data<-data.frame(scan(oldruns[1],what=list(yrfrac=0,co2i=0,prs=0,ltmp=0,btmp=0,rh=0,sfl=0,fl1=0,fl2=0,fl3=0,fl4=0, # nlin=0,ngas=0,up=0,sdco2=0,sdprs=0,sdsfl=0,hs2=0,hs1=0,ls1=0,ls2=0, # shpbn=0,pcorr=0,co2p=0,co2t=0,co2f=0,lin0=0,lin1=0,cur0=0,cur1=0,cur2=0,co2c=0,co2n=0,co2final=0))) data<-data.frame(scan(oldruns[1],what=list(yrfrac=0,nlin=0,ngas=0,co2final=0),skip=1)) ## temporary fix for change to year.frac from day since dec. 31, 2004 data$yrfrac[data$yrfrac<366]=(data$yrfrac[data$yrfrac<366]-1)/365+2005 data$yrfrac[data$yrfrac<731]=(data$yrfrac[data$yrfrac<731]-365-1)/365+2006 data$yrfrac[data$yrfrac<1096]=(data$yrfrac[data$yrfrac<1096]-365-365-1)/365+2007 data$yrfrac[data$yrfrac<1462]=(data$yrfrac[data$yrfrac<1462]-365-365-365-1)/366+2008 data$yrfrac[data$yrfrac<2000]=(data$yrfrac[data$yrfrac<2000]-365-365-365-366-1)/365+2009 data$yrfrac[data$yrfraceco]<-NA data<-na.omit(data) if(length(attributes(data)$na.action)>=1){ print(paste('Removed ',length(unique(attributes(data)$na.action)),' record(s) > eco or < ico from',i)) } if(length(data$yrfrac)==0){ nlines<-nlines[locs!=i] locs<-locs[locs!=i] } attributes(data)$na.action<-NULL assign(paste(i,'data',sep=''),data) minyrfrac<-min(c(minyrfrac,data$yrfrac),na.rm=T) maxyrfrac<-max(c(maxyrfrac,data$yrfrac),na.rm=T) } else { nlines<-nlines[locs!=i] locs<-locs[locs!=i] } if(minyrfrac==1E36&maxyrfrac==0){print(paste('NO DATA FOR INTERVAL FROM',i))} } # times do not line up perfectly because of different tubing lag times and possible clock corrections, so approx all to even 5 minutes #intdoyc<-floor(mindoyc)+c(1:((ceiling(maxdoyc)-floor(mindoyc))*288))/288 mindatetime=as.POSIXlt(strptime(paste(ye,'-01-01 00:00:00',sep=''),format='%Y-%m-%d %H:%M:%S',tz="GMT")+60*60*24*dpy*(minyrfrac-ye)) maxdatetime=as.POSIXlt(strptime(paste(ye,'-01-01 00:00:00',sep=''),format='%Y-%m-%d %H:%M:%S',tz="GMT")+60*60*24*dpy*(maxyrfrac-ye)) mindayfrac=as.numeric(as.POSIXct(mindatetime))/60/60/24 # as.POSIXct is in seconds, dayfrac is in days since 1970 maxdayfrac=as.numeric(as.POSIXct(maxdatetime))/60/60/24 # as.POSIXct is in seconds, dayfrac is in days since 1970 intdayfrac=seq(trunc(mindayfrac),trunc(maxdayfrac)+1,5/60/24) # every 5 min between start and end intdatetime=as.POSIXlt(intdayfrac*60*60*24,origin='1970-01-01',tz="GMT") intyrfrac=round(intdatetime$year+1900+as.numeric(difftime(intdatetime,strptime(paste(intdatetime$year+1900,'-01-01 00:00:00',sep=''),format='%Y-%m-%d %H:%M:%S',tz="GMT"),units='days'))/dpy,7) allco21<-NULL allco22<-NULL allco23<-NULL allco24<-NULL allco25<-NULL allco26<-NULL print('Interpolating records to common 5-min timeseries') #sink(type="message") #sink() #close(log) #browser() for(i in locs){ print(i) lnum<-c(1:length(locs))[i==locs] tempyrfrac<-get(paste(i,'data',sep=''))$yrfrac tempco2<-get(paste(i,'data',sep=''))$co2final tempngas<-get(paste(i,'data',sep=''))$ngas tempnlin<-get(paste(i,'data',sep=''))$nlin if(sum(tempnlin==1&tempngas==1)>1){ intco21<-approx(tempyrfrac[tempnlin==1&tempngas==1],tempco2[tempnlin==1&tempngas==1],intyrfrac,rule=2)$y for(j in c(1:length(intyrfrac))){ # screen for gaps > 5 min #if(min(abs(intdoyc[j]-tempdoyc[tempnlin==1&tempngas==1]),na.rm=T)>0.00347){ intco21[j]<-NA} if(min(abs(intyrfrac[j]-tempyrfrac[tempnlin==1&tempngas==1]),na.rm=T)>5/60/24/365){ intco21[j]<-NA} # 5-minutes # |all(intdoyc[j]-tempdoyc[tempnlin==1&tempngas==1]>0,na.rm=T)| # would also screen for unbounded extrapolation # all(intdoyc[j]-tempdoyc[tempnlin==1&tempngas==1]<0,na.rm=T) # but I do not think this is necessary (why in conc?) } }else{ intco21<-rep(NA,length(intyrfrac)) } if(sum(tempnlin==2&tempngas==1)>1){ intco22<-approx(tempyrfrac[tempnlin==2&tempngas==1],tempco2[tempnlin==2&tempngas==1],intyrfrac,rule=2)$y for(j in c(1:length(intyrfrac))){ # screen for gaps > 5 min if(min(abs(intyrfrac[j]-tempyrfrac[tempnlin==2&tempngas==1]),na.rm=T)>5/60/24/365){ intco22[j]<-NA} } }else{ intco22<-rep(NA,length(intyrfrac)) } if(sum(tempnlin==3&tempngas==1)>1){ intco23<-approx(tempyrfrac[tempnlin==3&tempngas==1],tempco2[tempnlin==3&tempngas==1],intyrfrac,rule=2)$y for(j in c(1:length(intyrfrac))){ # screen for gaps > 5 min if(min(abs(intyrfrac[j]-tempyrfrac[tempnlin==3&tempngas==1]),na.rm=T)>5/60/24/365){ intco23[j]<-NA} } }else{ intco23<-rep(NA,length(intyrfrac)) } if(nlines[lnum]==3){ if(sum(tempnlin==4)>1){ intco24<-approx(tempyrfrac[tempnlin==4&tempngas==1],tempco2[tempnlin==4&tempngas==1],intyrfrac,rule=2)$y for(j in c(1:length(intyrfrac))){ # screen for gaps > 2.5000 min to try to get only get 1 value for each LT run if(min(abs(intyrfrac[j]-tempyrfrac[tempnlin==4]),na.rm=T)>2.5/60/24/365){ intco24[j]<-NA} } }else{ intco24<-rep(NA,length(intyrfrac)) } }else if(nlines[lnum]==5){ if(sum(tempnlin==4&tempngas==1)>1){ intco24<-approx(tempyrfrac[tempnlin==4&tempngas==1],tempco2[tempnlin==4&tempngas==1],intyrfrac,rule=2)$y for(j in c(1:length(intyrfrac))){ # screen for gaps > 5 min if(min(abs(intyrfrac[j]-tempyrfrac[tempnlin==4&tempngas==1]),na.rm=T)>5/60/24/365){ intco24[j]<-NA} } }else{ intco24<-rep(NA,length(intyrfrac)) } if(sum(tempnlin==5&tempngas==1)>1){ intco25<-approx(tempyrfrac[tempnlin==5&tempngas==1],tempco2[tempnlin==5&tempngas==1],intyrfrac,rule=2)$y for(j in c(1:length(intyrfrac))){ # screen for gaps > 5 min if(min(abs(intyrfrac[j]-tempyrfrac[tempnlin==5&tempngas==1]),na.rm=T)>5/60/24/365){ intco25[j]<-NA} } }else{ intco25<-rep(NA,length(intyrfrac)) } if(sum(tempnlin==6)>1){ intco26<-approx(tempyrfrac[tempnlin==6&tempngas==1],tempco2[tempnlin==6&tempngas==1],intyrfrac,rule=2)$y for(j in c(1:length(intyrfrac))){ # screen for gaps > 2.5000 min to try to get only get 1 value for each LT run if(min(abs(intyrfrac[j]-tempyrfrac[tempnlin==6]),na.rm=T)>2.5/60/24/365){ intco26[j]<-NA} } }else{ intco26<-rep(NA,length(intyrfrac)) } }else{ print("INVALID NUMBER OF INLETS (NLINES = 3 OR 5)") # allow 2 or 4? } allco21<-cbind(allco21,intco21) allco22<-cbind(allco22,intco22) allco23<-cbind(allco23,intco23) allco24<-cbind(allco24,intco24) if(nlines[lnum]==5){ allco25<-cbind(allco25,intco25) allco26<-cbind(allco26,intco26) } # allow 2 or 4? } # READ IN INLET DATA print('Reading Inlet Data') inlets<-scan('../../aircoa_inlets.txt',what=list(loc='',unit='',date=0,time='',L1=0,L2=0,L3=0,L4=0,L5=0,hz1=0,hz2=0,hz3=0,hz4=0,hz5=0,upl=0),skip=1) # line heights (1-5), horizontal distance to tower, and # of units per line inlets$dayf<-as.numeric(substr(inlets$time,1,2))/24+as.numeric(substr(inlets$time,3,4))/60/24 topco2<-NULL ltco2<-NULL hghts<-matrix(c('Line 1','Line 2','Line 3','Line 4','Line 5'),ncol=5,nrow=length(locs),byrow=T) toplin<-rep(NA,length(locs)) for(i in c(1:length(locs))){ for(j in c(1:length(inlets$loc))){ if(inlets$loc[j]==locs[i]&inlets$date[j]+inlets$dayf[j]<(as.numeric(start)+ico)){ hghts[i,]<-c(paste(inlets$L1[j],'m'),paste(inlets$L2[j],'m'),paste(inlets$L3[j],'m'),paste(inlets$L4[j],'m'),paste(inlets$L5[j],'m')) temp<-c(1:5)[c(inlets$L1[j],inlets$L2[j],inlets$L3[j],inlets$L4[j],inlets$L5[j])==max(c(inlets$L1[j],inlets$L2[j],inlets$L3[j],inlets$L4[j],inlets$L5[j]))] toplin[i]<-temp[length(temp)] # if there is a tie for heighest inlet, break with inlet number } # assumes inlet heights are not changing within comp period } topco2<-cbind(topco2,get(paste('allco2',toplin[i],sep=''))[,i]) ltco2<-cbind(ltco2,get(paste('allco2',(nlines[i]+1),sep=''))[,i]) } # READ IN LT DATA ##caldir<-'/h/eol/stephens/NCAR/CALFAC/' ##source(paste(caldir,'calfac_lkup.r',sep='')) cflkup<-function(cylid='JJ1348'){ system(paste('grep -s ',cylid,' ../../aircoa_caldat.txt > temp.txt'),intern=T) #calrec<-scan('temp.txt',what=list(id='',date='',run='',pos='',gmdco2=0,siggmd=0,nco2=0),quiet=T,flush=T,fill=T) ### work with new columns? calrec<-scan('temp.txt',what=list(id='',date='',run='',pos='',gmdco2=0,siggmd=0,slpgmd=0,nco2=0),quiet=T,flush=T,fill=T) calrec<-data.frame(calrec) calrec<-calrec[!is.element(calrec$run,letters[14:26]),] # excludes experimental runs (n-z) if(nrow(calrec)>0){ row.names(calrec)<-as.character(c(1:length(calrec$id))) } write(t(calrec),'cflkup.out',ncol=7) system('rm temp.txt') } cylrec<-scan('../../aircoa_cylrec.txt',what=list(loc='',unit='',date=0,time=0,HS2='',HS1='',LS1='',LS2='',LT=''),skip=1) # clump and ave all aircoa cyls source('../../aircoa_ave_cals.r') ltval<-rep(NA,length(locs)) #ltval=NULL # will have columns for each loc and nrow=length(intyrfrac) crpos=1 for(i in c(1:length(locs))){ foundlt<-F for(j in c(1:length(cylrec$loc))){ if(cylrec$loc[j]==locs[i]&cylrec$date[j]<=as.numeric(start)){ # for now just use LT online at 'start' #### improve to allow LT changes mid-month crpos<-j foundlt<-T } } if(!foundlt|cylrec$LT[crpos]=='NoneUsed'){ # if comp run starts before site or LT was installed use first real LT for(j in c(length(cylrec$loc):1)){ if(cylrec$loc[j]==locs[i]&cylrec$LT[j]!='NoneUsed'){ crpos<-j foundlt<-T } } } if(foundlt){ cflkup(cylrec$LT[crpos]) if(F){ # replacing with aircoa_ave_cals.r calrec<-scan('cflkup.out',what=list(id='',date='',run='',pos='',gmdco2=0,siggmd=0,slpgmd=0,nco2=0)) ### work with new columns? # calrec<-scan('cflkup.out',what=list(file='',pos='',id='',co2=0,sigma=0,gmdco2=0,siggmd=0,nco2=0,o2=0,o2sigma=0,no2=0,date='',run='')) } else { #write('id date gmdco2 siggmd avesig nrun ntot','aircoa_caldat_ave.txt') calrec=read.table('aircoa_caldat_ave.txt',header=T,colClasses=c('character','character',NA,NA,NA,NA,NA)) } foundlt<-F for(j in c(1:length(calrec$id))){ if(as.numeric(calrec$date[j])as.numeric(ambient13cend[i])]=(1-0.028*0.011*1/3)*calrec$gmdco2[calrec$id==ambient13cid[i]&as.numeric(calrec$date)>as.numeric(ambient13cend[i])] # was once an ambient 13C gas but since refilled by Scott-Marrin (now has light 13C) print(calrec$id[calrec$id==ambient13cid[i]&as.numeric(calrec$date)>as.numeric(ambient13cend[i])]) print(calrec$date[calrec$id==ambient13cid[i]&as.numeric(calrec$date)>as.numeric(ambient13cend[i])]) } } # DO PLOTTING #lcols<-c('Red','Green','Blue','Purple') lcols<-c('Red','Green','Blue','Cyan','Orange') scols<-c('Dark Blue','Dark Green','Brown','Orange','Cyan','Dark Red') # loops after 6 sites # PLOT CONCENTRATIONS BY SITE print('Plotting Concentrations by Site') if(png){png(paste("comp_",start,"_",end,outflag,"_con.png",sep=''),width=950,height=550)}else{ bitmap(paste("comp_",start,"_",end,outflag,"_con.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)} # for 1-3 sites want 2 x 2, for 4-5 sites want 3 x 2, for 6-8 sites want 3 x 3, for 9-11 sites want 4 x 3 if(length(locs)<4){ par(mfrow=c(2,2)) }else if(length(locs)<6){ par(mfrow=c(3,2)) }else if(length(locs)<9){ par(mfrow=c(3,3)) }else if(length(locs)<12){ par(mfrow=c(4,3)) }else{ par(mfrow=c(4,4)) } par(mar=c(3,3,2,1)+.1) #defaults 5,4,4,2 (bot,lef,top,rig) par(mgp=c(2,1,0)) #defaults 3,1,0 par(oma=c(0,0,2,0)) #sink(type="message") #sink() #close(log) #browser() for(j in c(1:length(locs))){ #plot(intyrfrac,allco21[,j],main=locs[j],xlab='Year',ylab='Measured CO2 (ppm)',type='n',xlim=c(ico,eco),ylim=c(360,480)) plot(intdatetime,allco21[,j],main=locs[j],xlab='Date',ylab='Measured CO2 (ppm)',type='n',ylim=c(360,480)) for (i in c(1:nlines[j])){ # step through inlet lines and do LT separately if(any(!is.na(get(paste('allco2',i,sep=''))[,j]))){ lines(intdatetime[!is.na(get(paste('allco2',i,sep=''))[,j])],get(paste('allco2',i,sep=''))[,j][!is.na(get(paste('allco2',i,sep=''))[,j])],col=lcols[i]) } } if(any(!is.na(get(paste('allco2',(nlines[j]+1),sep=''))[,j]))){ lines(intdatetime[!is.na(get(paste('allco2',(nlines[j]+1),sep=''))[,j])],get(paste('allco2',(nlines[j]+1),sep=''))[,j][!is.na(get(paste('allco2',(nlines[j]+1),sep=''))[,j])],col='Purple') } legend(mean(intdatetime,na.rm=T),470,c(hghts[j,(1:nlines[j])],'LT'),lty=rep(1,(nlines[j]+1)),col=c(lcols[1:nlines[j]],'Purple'),ncol=2) } plot(intdatetime,topco2[,1],main='Highest Level at Each Site',xlab='Date',ylab='Measured CO2 (ppm)',type='n',ylim=c(360,480)) for(j in c(1:length(locs))){ if(any(!is.na(topco2[,j]))){ #print(c(j,length(intyrfrac[!is.na(topco2[,j])]),length(topco2[,j][!is.na(get(paste('allco2',nlines[j],sep=''))[,j])]),length(intyrfrac),length(intyrfrac[is.na(topco2[,j])]))) lines(intdatetime[!is.na(topco2[,j])],topco2[,j][!is.na(topco2[,j])],col=scols[j]) ##lines(intdoyc[!is.na(topco2[,j])],topco2[,j][!is.na(get(paste('allco2',nlines[j],sep=''))[,j])],col=scols[j]) } } if(length(locs)>1){ legend(mean(intdatetime,na.rm=T),470,paste(locs,' ',diag(hghts[,toplin])),lty=rep(1,length(locs)),col=scols[1:length(locs)],ncol=2,cex=0.85) }else{ legend(mean(intdatetime,na.rm=T),470,paste(locs,' ',hghts[,toplin]),lty=rep(1,length(locs)),col=scols[1:length(locs)],ncol=2,cex=0.85) } mtext(paste(start,' to ',end,'CO2 Concentrations by Site'),3,0.5,T) dev.off() # PLOT CONCENTRATIONS BY SITE ZOOMED print('Plotting Concentrations by Site Zoomed') if(png){png(paste("comp_",start,"_",end,outflag,"_con_zm.png",sep=''),width=950,height=550)}else{ bitmap(paste("comp_",start,"_",end,outflag,"_con_zm.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)} # for 1-3 sites want 2 x 2, for 4-5 sites want 3 x 2, for 6-8 sites want 3 x 3, for 9-11 sites want 4 x 3 if(length(locs)<4){ par(mfrow=c(2,2)) }else if(length(locs)<6){ par(mfrow=c(3,2)) }else if(length(locs)<9){ par(mfrow=c(3,3)) }else if(length(locs)<12){ par(mfrow=c(4,3)) }else{ par(mfrow=c(4,4)) } par(mar=c(3,3,2,1)+.1) #defaults 5,4,4,2 (bot,lef,top,rig) par(mgp=c(2,1,0)) #defaults 3,1,0 par(oma=c(0,0,2,0)) for(j in c(1:length(locs))){ if(nlines[j]==3){ plot(rep(intdatetime,3),c(allco21[,j],allco22[,j],allco23[,j]),main=locs[j],xlab='Date',ylab='Measured CO2 (ppm)',type='n') }else if(nlines[j]==5){ plot(rep(intdatetime,5),c(allco21[,j],allco22[,j],allco23[,j],allco24[,j],allco25[,j]),main=locs[j],xlab='Date',ylab='Measured CO2 (ppm)',type='n') }else{ print("INVALID NUMBER OF INLETS (NLINES = 3 OR 5)") # allow 2 or 4? } for (i in c(1:nlines[j])){ # step through inlet lines if(any(!is.na(get(paste('allco2',i,sep=''))[,j]))){ lines(intdatetime[!is.na(get(paste('allco2',i,sep=''))[,j])],get(paste('allco2',i,sep=''))[,j][!is.na(get(paste('allco2',i,sep=''))[,j])],col=lcols[i]) } } legend(mean(intdatetime,na.rm=T),(par('usr')[4]-par('usr')[3])*.9+par('usr')[3],hghts[j,(1:nlines[j])],lty=rep(1,nlines[j]),col=lcols[1:nlines[j]],ncol=2) } plot(rep(intdatetime,length(locs)),topco2,main='Highest Level at Each Site',xlab='Date',ylab='Measured CO2 (ppm)',type='n') for(j in c(1:length(locs))){ #if(any(!is.na(get(paste('allco2',nlines[j],sep=''))[,j]))){ if(any(!is.na(topco2[,j]))){ lines(intdatetime[!is.na(topco2[,j])],topco2[,j][!is.na(topco2[,j])],col=scols[j]) } } if(length(locs)>1){ legend(mean(intdatetime,na.rm=T),(par('usr')[4]-par('usr')[3])*.9+par('usr')[3],paste(locs,' ',diag(hghts[,toplin])),lty=rep(1,length(locs)),col=scols[1:length(locs)],ncol=2,cex=0.85) }else{ legend(mean(intdatetime,na.rm=T),(par('usr')[4]-par('usr')[3])*.9+par('usr')[3],paste(locs,' ',hghts[,toplin]),lty=rep(1,length(locs)),col=scols[1:length(locs)],ncol=2,cex=0.85) } mtext(paste(start,' to ',end,'CO2 Concentrations by Site Zoomed'),3,0.5,T) dev.off() # PLOT CONCENTRATIONS BY SITE ZOOMED 1 PANEL if(png){png(paste("comp_",start,"_",end,outflag,"_con_zm_1p.png",sep=''),width=950,height=550)}else{ bitmap(paste("comp_",start,"_",end,outflag,"_con_zm_1p.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)} par(mfrow=c(1,1)) par(mar=c(3,3,2,1)+.1) #defaults 5,4,4,2 (bot,lef,top,rig) par(mgp=c(2,1,0)) #defaults 3,1,0 par(oma=c(0,0,2,0)) plot(rep(intdatetime,length(locs)),topco2,main='Highest Level at Each Site',xlab='Date',ylab='Measured CO2 (ppm)',type='n') for(j in c(1:length(locs))){ #if(any(!is.na(get(paste('allco2',nlines[j],sep=''))[,j]))){ if(any(!is.na(topco2[,j]))){ lines(intdatetime[!is.na(topco2[,j])],topco2[,j][!is.na(topco2[,j])],col=scols[j]) } } legend(mean(intdatetime,na.rm=T),(par('usr')[4]-par('usr')[3])*.9+par('usr')[3],paste(locs,' ',hghts[,toplin]),lty=rep(1,length(locs)),col=scols[1:length(locs)],ncol=2) mtext(paste(start,' to ',end,'CO2 Concentrations for Highest Level at Each Site'),3,0.5,T) dev.off() # PLOT PM CONCENTRATIONS BY SITE print('Plotting PM Concentrations by Site') if(png){png(paste("comp_",start,"_",end,outflag,"_con_pm.png",sep=''),width=950,height=550)}else{ bitmap(paste("comp_",start,"_",end,outflag,"_con_pm.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)} par(mfrow=c(2,1)) par(mar=c(3,3,2,1)+.1) #defaults 5,4,4,2 (bot,lef,top,rig) par(mgp=c(2,1,0)) #defaults 3,1,0 par(oma=c(0,0,2,0)) c1=0.75 # 1800 GMT c2=0.8333 # 2000 GMT c3=0.9167 # 2200 GMT intdayf<-intdayfrac-trunc(intdayfrac) # fraction of each day #if(any(intdayf>=c1&intdayf<=c3&!is.na(topco2[,1]))){ # first loc must have data intdays<-trunc(intdayfrac) # days since 1970 pmdays<-sort(unique(intdays)) #sink(type="message") #sink() #close(log) #browser() plot(intdatetime[intdayf>=c1&intdayf<=c3],topco2[,1][intdayf>=c1&intdayf<=c3],main='1800-2000 UTC Averages',xlab='Date',ylab='Measured CO2 (ppm)',type='n',ylim=c(min(topco2[intdayf>=c1&intdayf<=c3,],na.rm=T),max(topco2[intdayf>=c1&intdayf<=c3,],na.rm=T))) # setting up plot using first loc for(j in c(1:length(locs))){ if(any(!is.na(topco2))){ # why check all topco2 not just loc? #pmyrfrac1<-rep(NA,length(pmdays)) pmyrfrac1<-rep(strptime(start,format='%y%m%d',tz="GMT"),length(pmdays)) pmval1<-rep(NA,length(pmdays)) for(i in c(1:length(pmdays))){ #pmyrfrac1[i]<-mean(intyrfrac[intdays==pmdays[i]&intdayf>=c1&intdayf<=c2],na.rm=T) pmyrfrac1[i]<-mean(intdatetime[intdays==pmdays[i]&intdayf>=c1&intdayf<=c2],na.rm=T) pmval1[i]<-mean(topco2[,j][intdays==pmdays[i]&intdayf>=c1&intdayf<=c2],na.rm=T) } points(pmyrfrac1,pmval1,type='b',col=scols[j]) } } plot(intdatetime[intdayf>=c1&intdayf<=c3],topco2[,1][intdayf>=c1&intdayf<=c3],main='2000-2200 UTC Averages',xlab='Date',ylab='Measured CO2 (ppm)',type='n',ylim=c(min(topco2[intdayf>=c1&intdayf<=c3,],na.rm=T),max(topco2[intdayf>=c1&intdayf<=c3,],na.rm=T))) # setting up plot using first loc for(j in c(1:length(locs))){ if(any(!is.na(topco2))){ # why check all topco2 not just loc? #pmyrfrac2<-rep(NA,length(pmdays)) pmyrfrac2<-rep(strptime(start,format='%y%m%d',tz="GMT"),length(pmdays)) pmval2<-rep(NA,length(pmdays)) for(i in c(1:length(pmdays))){ #pmyrfrac2[i]<-mean(intyrfrac[intdays==pmdays[i]&intdayf>=c2&intdayf<=c3],na.rm=T) pmyrfrac2[i]<-mean(intdatetime[intdays==pmdays[i]&intdayf>=c2&intdayf<=c3],na.rm=T) pmval2[i]<-mean(topco2[,j][intdays==pmdays[i]&intdayf>=c2&intdayf<=c3],na.rm=T) } points(pmyrfrac2,pmval2,type='b',col=scols[j]) } } legend(mean(intdatetime[intdayf>=c1&intdayf<=c3],na.rm=T),max(topco2[intdayf>=c1&intdayf<=c3,],na.rm=T),paste(locs,' ',hghts[,toplin]),lty=rep(1,length(locs)),pch=rep(1,length(locs)),col=scols[1:length(locs)],ncol=2) mtext(paste(start,' to ',end,'Highest Level at Each Site'),3,0.5,T) #} # if any afternoon dev.off() # CALCULATE AVERAGE DIURNAL CYCLES AT 5-MIN RESOLUTION # do 15-min as in aircoa_conc_days.r? print("Calculating diurnal cycles") diurdayf<-intdayfrac-trunc(intdayfrac) # fraction of each day mindiur<-1000 maxdiur<-0 for(j in c(1:length(locs))){ # diur<-list(dayf=NULL,co2diur1=NULL,co2diur2=NULL,co2diur3=NULL,co2diur2=NULL,co2diur3=NULL) diur<-list(dayf=NULL,co2diur=NULL) diur$dayf<-sort(unique(round(diurdayf,4))) for(i in c(1:length(diur$dayf))){ diur$co2diur[i]<-mean(topco2[,j][round(diurdayf,4)==diur$dayf[i]],na.rm=T) if(!is.na(diur$co2diur[i])){ if(diur$co2diur[i]maxdiur){maxdiur<-diur$co2diur[i]} } # diur$co2diur1[i]<-mean(allco21[,j][round(diurdayf,4)==diur$dayf[i]],na.rm=T) # diur$co2diur2[i]<-mean(allco22[,j][round(diurdayf,4)==diur$dayf[i]],na.rm=T) # diur$co2diur3[i]<-mean(allco23[,j][round(diurdayf,4)==diur$dayf[i]],na.rm=T) # if(!is.na(diur$co2diur3[i])){ # if(diur$co2diur3[i]maxdiur){maxdiur<-diur$co2diur3[i]} # } } assign(paste(locs[j],'diur',sep=''),diur) } # PLOT DIURNAL CYCLES print("Plotting diurnal cycles") if(png){png(paste("comp_",start,"_",end,outflag,"_diur.png",sep=''),width=950,height=550)}else{ bitmap(paste("comp_",start,"_",end,outflag,"_diur.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)} par(mfrow=c(1,1)) par(mar=c(3,3,2,1)+.1) #defaults 5,4,4,2 (bot,lef,top,rig) par(mgp=c(2,1,0)) #defaults 3,1,0 par(oma=c(0,0,2,0)) plot(diur$dayf*24,get(paste(locs[1],'diur',sep=''))$co2diur,type='n',ylab='Average CO2 (ppm)',xlab='Hour of Day (UTC)',ylim=c(mindiur,maxdiur)) for(j in c(1:length(locs))){ if(any(!is.na(get(paste(locs[j],'diur',sep=''))$co2diur))){ lines(diur$dayf[!is.na(get(paste(locs[j],'diur',sep=''))$co2diur)]*24,get(paste(locs[j],'diur',sep=''))$co2diur[!is.na(get(paste(locs[j],'diur',sep=''))$co2diur)],col=scols[j]) } } if(length(locs)>1){ legend(15,maxdiur,paste(locs,' ',diag(hghts[,toplin])),lty=rep(1,length(locs)),col=scols[1:length(locs)],ncol=2) }else{ legend(15,maxdiur,paste(locs,' ',hghts[,toplin]),lty=rep(1,length(locs)),col=scols[1:length(locs)],ncol=2) } mtext(paste(start,' to ',end,'Diurnal CO2 Cycles for Highest Level at Each Site'),3,0.5,T) dev.off() # PLOT LT COMPARISON print('Plotting LT Comparison') if(png){png(paste("comp_",start,"_",end,outflag,"_lt.png",sep=''),width=950,height=550)}else{ bitmap(paste("comp_",start,"_",end,outflag,"_lt.png",sep=''),width=9.36,height=5.4,res=100,type='png16',pointsize=10)} # for 1-3 sites want 2 x 2, for 4-5 sites want 3 x 2, for 6-8 sites want 3 x 3, for 9-11 sites want 4 x 3 if(length(locs)<4){ par(mfrow=c(2,2)) }else if(length(locs)<6){ par(mfrow=c(3,2)) }else if(length(locs)<9){ par(mfrow=c(3,3)) }else if(length(locs)<12){ par(mfrow=c(4,3)) }else{ par(mfrow=c(4,4)) } par(mar=c(3,3,2,1)+.1) #defaults 5,4,4,2 (bot,lef,top,rig) par(mgp=c(2,1,0)) #defaults 3,1,0 par(oma=c(0,0,2,0)) for(j in c(1:length(locs))){ if(any(!is.na(ltco2[,j]))){ plot(intdatetime,ltco2[,j],main=locs[j],xlab='Date',ylab='Measured CO2 (ppm)',type='n',ylim=ltval[j]+c(-1,1)) abline(h=ltval[j]) lines(intdatetime[!is.na(ltco2[,j])],ltco2[,j][!is.na(ltco2[,j])],col='Purple',type='b') legend(mean(intdatetime,na.rm=T),ltval[j]+0.9,c('Field Measured','Laboratory Assigned'),lty=c(1,1),col=c('Purple','Black')) } } plot(intdatetime,ltco2[,1],main='Subtracted from Assigned at Each Site',xlab='Date',ylab='Measured CO2 (ppm)',type='n',ylim=c(-1,1)) abline(h=0) for(j in c(1:length(locs))){ if(any(!is.na(ltco2[,j]))){ lines(intdatetime[!is.na(ltco2[,j])],ltco2[,j][!is.na(ltco2[,j])]-ltval[j],col=scols[j],type='b') } } if(any(!is.na(ltco2))){ legend(mean(intdatetime,na.rm=T),.9,locs,lty=rep(1,length(locs)),col=scols[1:length(locs)],ncol=2) mtext(paste(start,' to ',end,'LT Concentrations by Site'),3,0.5,T) } dev.off() ## Plot differences vs. T and P ? # OUTPUT MERGED FILE #outmat<-cbind(intdoyc,allco21,allco22,allco23,allco24,allco25,allco26,topco2) outmat<-cbind(intyrfrac,round(topco2,3),round(ltco2,3)) if(length(locs)>1){ # t1<-paste(locs,diag(hghts[,toplin]),collapse=' ') t1<-paste(locs,diag(hghts[,toplin])) }else{ # t1<-paste(locs,hghts[,toplin],collapse=' ') t1<-paste(locs,hghts[,toplin]) } #t2<-paste(locs,'LT',collapse=' ') t2<-paste(locs,'LT') #write(paste('DOY ',t1,' ',t2),paste('comp_',start,'_',end,outflag,".out",sep='')) #write(t(outmat),paste('comp_',start,'_',end,outflag,".out",sep=''),ncol=ncol(outmat),append=T) outmat2=rbind(c("DOY",t1,t2),outmat) # binding text somehow ensures 7 digits in yrfrac output write(t(outmat2),paste('comp_',start,'_',end,outflag,".out",sep=''),ncol=ncol(outmat2)) # CALCULATE OFFSETS write('SITE(S) LINE N MEANDIF MEDIANDIF STDEVDIF',paste('comp_',start,'_',end,outflag,".res",sep='')) write('',paste('comp_',start,'_',end,outflag,".res",sep=''),append=T) #done<-NULL #for(i in c(1:length(locs))){ # done<-c(done,i) # for(j in c(1:length(locs))){ # if(!any(j==done)){ # for(k in c(1:3)){ # diff<-get(paste('allco2',k,sep=''))[,i]-get(paste('allco2',k,sep=''))[,j] # write(paste(locs[i],'-',locs[j],k,' ',length(diff[!is.na(diff)]), # format(round(mean(diff,na.rm=T),3),nsmall=3), # format(round(median(diff,na.rm=T),3),nsmall=3), # format(round(sqrt(var(diff,na.rm=T)),3),nsmall=3)), # paste('comp_',start,'_',end,outflag,".res",sep=''),append=T) # } # k<-4 # diff<-get(paste('allco2',k,sep=''))[,i]-get(paste('allco2',k,sep=''))[,j] # chance of missing some if approx put at different time # write(paste(locs[i],'-',locs[j],k,' ',length(diff[!is.na(diff)]), # format(round(mean(diff,na.rm=T),3),nsmall=3), # format(round(median(diff,na.rm=T),3),nsmall=3), # format(round(sqrt(var(diff,na.rm=T)),3),nsmall=3)), # paste('comp_',start,'_',end,outflag,".res",sep=''),append=T) # write('',paste('comp_',start,'_',end,outflag,".res",sep=''),append=T) # } # } #} for(i in c(1:length(locs))){ write(paste(locs[i],'- cal LT ',length(ltco2[,i][!is.na(ltco2[,i])]), format(round(mean(ltco2[,i]-ltval[i],na.rm=T),3),nsmall=3), format(round(median(ltco2[,i]-ltval[i],na.rm=T),3),nsmall=3), format(round(sqrt(var(ltco2[,i]-ltval[i],na.rm=T)),3),nsmall=3)), paste('comp_',start,'_',end,outflag,".res",sep=''),append=T) } sink(type="message") sink() close(log) }