From b0039716d4e4333f99944196da960e5aee380003 Mon Sep 17 00:00:00 2001 From: John-C-Lin Date: Mon, 7 Mar 2022 15:28:14 -0700 Subject: [PATCH] Initial commit --- Fch4_simple.r | 782 +++++++++++++++++++++++++ Fch4_simple_multimettype.r | 288 +++++++++ README | 82 +++ Trend_foot_dCH4.r | 200 +++++++ average_foot_monthly.r | 140 +++++ compare_foots_HYSPLIT-STILT_vs_STILT.r | 150 +++++ compare_trajs_HYSPLIT-STILT_vs_STILT.r | 168 ++++++ convolve_foot_wellinfo.r | 194 ++++++ dCH4_vs_foot_wellinfo.r | 524 +++++++++++++++++ export_Fch4simple.r | 33 ++ export_Uintah_tseries.r | 29 + grid_well_info_foot.r | 149 +++++ merge_CH4_UVwind_wellinfo.r | 79 +++ plot_CH4tseries_diurnal_seasonal.r | 230 ++++++++ plot_foot_and_wells.r | 155 +++++ plot_foot_wellinfo_convolved.r | 116 ++++ plot_foot_wellinfo_tseries_trend.r | 263 +++++++++ plot_gridded_wellinfoV1.r | 115 ++++ plot_tseries_and_COVID.r | 245 ++++++++ 19 files changed, 3942 insertions(+) create mode 100644 Fch4_simple.r create mode 100644 Fch4_simple_multimettype.r create mode 100644 README create mode 100644 Trend_foot_dCH4.r create mode 100644 average_foot_monthly.r create mode 100644 compare_foots_HYSPLIT-STILT_vs_STILT.r create mode 100644 compare_trajs_HYSPLIT-STILT_vs_STILT.r create mode 100644 convolve_foot_wellinfo.r create mode 100644 dCH4_vs_foot_wellinfo.r create mode 100644 export_Fch4simple.r create mode 100644 export_Uintah_tseries.r create mode 100644 grid_well_info_foot.r create mode 100644 merge_CH4_UVwind_wellinfo.r create mode 100644 plot_CH4tseries_diurnal_seasonal.r create mode 100644 plot_foot_and_wells.r create mode 100644 plot_foot_wellinfo_convolved.r create mode 100644 plot_foot_wellinfo_tseries_trend.r create mode 100644 plot_gridded_wellinfoV1.r create mode 100644 plot_tseries_and_COVID.r diff --git a/Fch4_simple.r b/Fch4_simple.r new file mode 100644 index 0000000..89f4a3a --- /dev/null +++ b/Fch4_simple.r @@ -0,0 +1,782 @@ +# Simple STILT analyses of Uintah Basin methane emissions +# Idea is to calculate AFTERNOON enhancements (CH4_HPL - CH4_FRU), & then divide by total footprint strength: +# F_CH4 = (CH4_HPL - CH4_FRU)/sum(foot_i,j) +# +# Then see how stable F_CH4 is over different seasons and years. +# V2(200911): use new HYSPLIT-STILT executable instead of STILT executable +# V3(200920): instead of annual production data, use monthly production data +# V4(210320): add CSP as receptor, save results that are site-specific +# V5(210412): add option to fill in gaps in FRU time series and merge in wind obs + sim to filter out times when HRRR-STILT was off using quadrant method +# V6(210416): filter out days based on within 45o of observed wind direction (better than quadrant method above, since not throw way data close to quadrant boundaries) +# V7(210422): calculate Energy-Normalized Methane Average (ENMA) emissions, make CH4 volume content explicit variable, combine monthly and annual emissions (for paper figure), create function to convert from [MCF NatGas] => [kg CH4] +# V8(210423): introduce 'mettype' to distinguish between met files driving STILT (to assess transport error)--e.g., "HRRR", "NARR", "NAM12" +# V9(211004): sensitivity analysis to change model domain +# May 29th, 2019 by John C. Lin (John.Lin@utah.edu) + +require(ncdf4); require(fields); require(geosphere) +#################### +YEARs <- 2015:2020 +MONsub <- 4:9 #subset of months to calculate fluxes and leak rates +#MONsub <- 6:8 #subset of months to calculate fluxes and leak rates +HRs<-20:23 #[UTC] only analyze afternoon hours, following Foster papers +#HRs<-17:23 #[UTC] only analyze afternoon hours, following Foster papers +SITE<-"HPL" #HPL is site to focus on for calculating long-term F_CH4 +#SITE<-"ROO" +#SITE<-"CSP" +if(SITE=="CSP")YEARs <- c(2016,2019,2020) +if(SITE=="ROO")YEARs <- 2015:2019 + +obsdir <- "/uufs/chpc.utah.edu/common/home/lin-group4/jcl/SimCity/" +winddir <- "/uufs/chpc.utah.edu/common/home/lin-group7/jcl/Transporterr_stiltread/Output" #where wind obs and sim values are found + +mettype <- "HRRR" +#mettype <- "HRRR_2000particles" +#mettype <- "NAM12" +#mettype <- "WRF27" +# Where STILT output is found +# a) STILT executable +# STILTdir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_STILT/out/" +# b) HYSPLIT-STILT executable +#mettype <- toupper(mettype) +if(mettype=="HRRR")STILTdir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_HYSPLIT-STILT_HRRR/out/" +if(mettype=="HRRR_2000particles")STILTdir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_HYSPLIT-STILT_HRRR/out/" +if(mettype=="NAM12")STILTdir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_HYSPLIT-STILT_NAM12/out/" +if(mettype=="WRF27")STILTdir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_HYSPLIT-STILT_WRF27km/out/" +print(paste("STILTdir=",STILTdir)) + +resultname<-paste("Fch4_",SITE,"_daily_",mettype,".rds",sep="") + +Nday.min <- 10 # minimum days necessary for a monthly average to be retained (otherwise assigned NA) + +preparedata.TF <- TRUE # run lines to prepare data? +fillFRU.TF <- TRUE # fill in gaps in background (FRU) time series? +filterUV.TF <- TRUE # filter times based on U/V (filter out times when HRRR is off--i.e., large transport errors) ? +domainsens.TF <- FALSE # sensitivity analysis varyiing the domain size? +CH4.VOLFRAC <- 0.89 # volume fraction of CH4 in natural gas of 0.89 [Karion et al. 2013] +XLIMS<-c(-110.6,-109) # lon range of Uintah Basin--default +YLIMS<-c(39.9,40.5) # lat range of Uintah Basin--default + +if(domainsens.TF){ + print("vary domain from default; remember to set preparedata.TF to TRUE!") + # preparedata.TF <- TRUE + XLIMS<-c(-110,-109) # lon range of Uintah Basin--change for sensitivity analysis (in response to Reviewer #1) +} # if(domainsens.TF){ + +#################### + +if(filterUV.TF&substring(mettype,1,4)!="HRRR")stop(paste0("Only extracted HRRR wind vectors; filterUV.TF needs to be set to FALSE for mettype=",mettype)) + +if(SITE=="HPL"){LAT <- 40.1434; LON <- -109.4680; zagl <- 4.06; WINDsite <- "UBHSP"; WINDsite.2015 <- "A1633"} # UBHSP site not available in 2015, so use A1633 (RedWash) site instead in 2015 +if(SITE=="ROO"){LAT <- 40.2941; LON <- -110.0090; zagl <- 4.06; WINDsite <- "QRS"} +if(SITE=="CSP"){LAT <- 40.0509; LON <- -110.0194; zagl <- 4.06; WINDsite <- "UBCSP"} + +# Function that converts from natural gas volume in thousands of [MCF] to [kg of CH4] +MCF2kg <- function(MCF,CH4.VOLFRAC=0.89){ + NatGas<-MCF*CH4.VOLFRAC #[MCF natural gas] => [MCF CH4], using volume fraction of CH4 in natural gas of 0.89 [Karion et al. 2013] + NatGas<-NatGas*1000 #[MCF CH4] => [ft^3 CH4] + NatGas<-NatGas*0.0283 #[ft^3 CH4] => [m^3 CH4] + R.CH4<-8.3143*1000/16.043 #Ideal Gas constant of CH4 [J/K/kg]: (Rg/mCH4), where Rg is universal gas constant and mCH4 is molar mass of CH4 + rho.CH4<-(101300)/(R.CH4*288.7) #density of CH4 [kg/m3], using P=101.3kPa and T=288.7K [Karion et al. 2013] + NatGas<-NatGas*rho.CH4 #[m^3 CH4] => [kg CH4] + return(NatGas) +} # MCF2kg <- function(MCF,CH4.VOLFRAC=0.89){ + +if(preparedata.TF){ +######################################################### +#I. Calculate CH4 enhancements over background site (FRU) +dat.all<-NULL +for(i in 1:length(YEARs)){ + objname<-paste0("SimCity_CH4_allsites_hrly_",YEARs[i]) + print(paste("Reading in.....",objname)) + tmp<-getr(objname,path=obsdir)[,c("Time",paste0("CH4_",c(SITE,"FRU")))] + print(colnames(tmp)) + dat.all<-rbind(dat.all,tmp) + gc() +} #for(i in 1:length(YEARs)){ + +# filter for specified hrs [UTC] +YYYYMMDDHH <- format(dat.all$Time,"%Y%m%d%H") +dat.all <- data.frame(YYYYMMDDHH,dat.all) +sel <- as.numeric(substring(dat.all$YYYYMMDDHH,9,10))%in%HRs +dat.all <- dat.all[sel,] + +# fill in gaps in background (FRU) time series +if(fillFRU.TF){ + dat <- dat.all + YYYYMM <- substring(dat$YYYYMMDDHH,1,6) + dat <- data.frame(YYYYMM,dat) + # calculate % of NAs in each month + N <- tapply(dat$YYYYMM,dat$YYYYMM,length) + sumNA <- function(x)return(sum(is.na(x))) + N.NA <- tapply(dat$CH4_FRU,dat$YYYYMM,sumNA) + print("% of data with NAs at FRU in each month:") + print(round(100*N.NA/N,2)) + + # calculate CV (coefficient of variation)--what % is variability in FRU compared to magnitude of CH4 enhancement? + xsigma <- sqrt(tapply(dat$CH4_FRU,dat$YYYYMM,var,na.rm=T)) + dCH4.ave <- tapply(dat[,paste0('CH4_',SITE)]-dat$CH4_FRU,dat$YYYYMM,mean,na.rm=T) + print('% stdev(FRU.hrly)/mean(SITE-FRU):') + print(round(100*xsigma/dCH4.ave,2)) + # calculate CV for DAILY-AVERAGED CH4 (for selected UThrs) + FRU.day <- tapply(dat$CH4_FRU,substring(dat$YYYYMMDDHH,1,8),mean,na.rm=T) + dCH4.day <- tapply(dat[,paste0('CH4_',SITE)]-dat$CH4_FRU,substring(dat$YYYYMMDDHH,1,8),mean,na.rm=T) + YYYYMMDD <- names(FRU.day) + xsigma <- sqrt(tapply(FRU.day,substring(YYYYMMDD,1,6),var,na.rm=T)) + dCH4.ave <- tapply(dCH4.day,substring(YYYYMMDD,1,6),mean,na.rm=T) + print('% stdev(FRU.day)/mean(SITE-FRU):') + print(round(100*xsigma/dCH4.ave,2)) + + # fill in gaps in FRU time series with monthly average + print("fill in gaps in FRU time series with monthly average") + f <- function(x){ + x.ave <- mean(x,na.rm=T) + x[is.na(x)] <- x.ave + return(x) + } # f <- function(x){ + FRU.gapfilled <- unlist(tapply(dat$CH4_FRU,dat$YYYYMM,f)) + if(TRUE){ + FRU.monave <- tapply(dat$CH4_FRU,dat$YYYYMM,mean,na.rm=T) + tt <- as.numeric(substring(names(FRU.monave),1,4))+as.numeric(substring(names(FRU.monave),5,6))/12 + dev.new(); plot(tt,FRU.monave,pch=16,type="o") + dev.new(); plot(dat[,c("Time","CH4_FRU")],pch=16) + isNA <- is.na(dat$CH4_FRU) + points(dat$Time[isNA],FRU.gapfilled[isNA],pch=16,col="orange") + legend("topright",c("all","gap-filled"),pch=16,col=c("black","orange")) + } #if(FALSE){ + dat$CH4_FRU <- FRU.gapfilled + dat.all <- dat +} # if(fillFRU.TF){ + +# calculate CH4 enhancement over FRU [ppm] +dCH4 <- dat.all[,paste0('CH4_',SITE)]-dat.all$CH4_FRU +dat.all <- data.frame(dCH4,dat.all) + +#calculate DAILY dCH4 +dCH4 <- data.frame(Time=dat.all[,"Time"],dCH4) #[ppm] +Time.day<-format(dCH4[,"Time"],format="%Y%m%d") +dCH4.ave<-tapply(dCH4[,"dCH4"],Time.day,mean,na.rm=T) +Time<-strptime(names(dCH4.ave),"%Y%m%d",tz="GMT") +dCH4.ave<-data.frame(YYYYMMDD=names(dCH4.ave),Time,dCH4.ave,stringsAsFactors=FALSE) +dev.new();par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3) +plot(dCH4.ave[,c("Time","dCH4.ave")],main=paste(SITE,"\nUThrs:",paste(HRs,collapse=" ")),pch=16, + xlab="Time",ylab="Daily Averge Enhancement over Background [ppm]") +dCH4.ave<-data.frame(dCH4.ave,footsum=NA,foot.basin=NA) +sel<-substring(dCH4.ave[,"YYYYMMDD"],1,6)%in%c("201501","201502","201503","201504","201505","201506") #lack of obs--remove +dCH4.ave<-dCH4.ave[!sel,] + + +dat.all <- data.frame(dat.all,footsum=NA,foot.basin=NA) +######################################################### +#II. Calculate total footprint strengths and subset of footprint strength that just fall within Uintah Basin +footdir<-paste0(STILTdir,"/",SITE,"/footprints/") +if(mettype=="HRRR_2000particles") footdir<-paste0(STILTdir,"/",SITE,"_2000particles/footprints/") #sensitive analysis with 10x base simulation particle number +for(i in 1:nrow(dat.all)){ + YYYYMMDDHH<-dat.all[i,"YYYYMMDDHH"] + footfile<-list.files(footdir,pattern=YYYYMMDDHH) + if(length(footfile)==0){print(paste("Cannot find any STILT footprints for",YYYYMMDDHH,"; skip"));next} + if(length(footfile)> 1){print(paste("More than ONE STILT footprint file found for",YYYYMMDDHH,"; skip"));next} + xfile <- paste0(footdir,footfile) + print(paste("Reading in....",footfile)) + xnc <- nc_open(xfile) + xlon <- ncvar_get(xnc,varid="lon"); xlat<-ncvar_get(xnc,varid="lat") + xfoot <- apply(ncvar_get(xnc,varid="foot"),c(1,2),sum) + # total footprint ALL doomain + footsum <- sum(xfoot) #[ppm/(umole/m2/s)] + dat.all[i,"footsum"] <- footsum + # total footprint within Uintah Basin + sel.x <- xlon>=XLIMS[1]&xlon<=XLIMS[2] + sel.y <- xlat>=YLIMS[1]&xlat<=YLIMS[2] + foot.basin <- sum(xfoot[sel.x,sel.y]) + dat.all[i,"foot.basin"]<-foot.basin + + saveRDS(dat.all,file=paste0("Fch4_",SITE,"_",mettype,".rds")) + nc_close(xnc) + gc() +} #for(i in 1:nrow(dCH4)){ + + +######################################################### +#III. Merge with wind data (obs and HRRR simulated winds) to enable times when HRRR are off to be filtered out +dat.all <- readRDS(paste0("Fch4_",SITE,"_",mettype,".rds")) +wind.all <- NULL +for(yy in 1:length(YEARs)){ + # read in HRRR winds + observed values + windname <- paste0("HRRR_obs_",YEARs[yy],"01010000to",YEARs[yy],"12312300.RDS") + if(YEARs[yy]=="2018")windname <- paste0("HRRR_obs_",YEARs[yy],"01010000to",YEARs[yy],"12310000.RDS") # typo during 2018 + colnms <- c("Time",paste0(c("Usim.","Vsim.","Tsim.","Uobs.","Vobs.","Tobs."),WINDsite)) + # UBHSP site not available in 2015, so use A1633 (RedWash) site instead in 2015 + if(SITE=="HPL")colnms <- c("Time",paste0(c("Usim.","Vsim.","Tsim.","Uobs.","Vobs.","Tobs."),WINDsite), + paste0(c("Usim.","Vsim.","Tsim.","Uobs.","Vobs.","Tobs."),WINDsite.2015)) + if(file.exists(paste0(winddir,"/",windname))) { + print(paste("Reading in.....",windname)) + tmp0 <- readRDS(paste0(winddir,"/",windname))$dat + tmp <- matrix(NA,nrow=nrow(tmp0),ncol=length(colnms)) + colnames(tmp) <- colnms + tmp <- data.frame(tmp) + colnms.sel <- colnames(tmp0)%in%colnms + tmp[,colnames(tmp0)[colnms.sel]] <- tmp0[,colnms.sel] + } else { + print(paste("Doesn't exists:",windname,"; replaces with NAs")) + tmp.dat <- dat.all[substring(dat.all$YYYYMM,1,4)==YEARs[yy],] + tmp <- matrix(NA,nrow=nrow(tmp.dat),ncol=length(colnms)) + colnames(tmp) <- colnms + tmp <- data.frame(tmp) + tmp$Time <- tmp.dat$Time + } # if(file.exists(paste0(winddir,"/",windname)){ + wind.all <- rbind(wind.all,tmp) + + gc() +} # for(yy in 1:length(YEARs)){ + + # merge dat.all with wind data + result <- merge(x=dat.all,y=wind.all,by='Time',all.x=TRUE,all.y=FALSE) + dat.all <- result + saveRDS(dat.all,file=paste0("Fch4_",SITE,"_",mettype,".rds")) + print(paste0("Fch4_",SITE,"_",mettype,".rds"," written out")) + write.csv(dat.all,file=paste0("Fch4_",SITE,"_",mettype,".csv"),row.names=FALSE) + print(paste0("Fch4_",SITE,"_",mettype,".csv"," written out")) + +} # if(preparedata.TF){ + +######################################################### +#IV. Filter times based on U/V (filter out times when HRRR is off--i.e., large transport errors) ? +dat.all <- readRDS(paste0("Fch4_",SITE,"_",mettype,".rds")) +if(filterUV.TF){ + dat <- dat.all + Usim <- dat[,paste0("Usim.",WINDsite)]; Vsim <- dat[,paste0("Vsim.",WINDsite)] + Uobs <- dat[,paste0("Uobs.",WINDsite)]; Vobs <- dat[,paste0("Vobs.",WINDsite)] + if(SITE=="HPL"){ + # UBHSP site not available in 2015, so use A1633 (RedWash) site instead in 2015 + sel <- substring(dat$YYYYMM,1,4)=="2015" + Usim[sel] <- dat[sel,paste0("Usim.",WINDsite.2015)] + Vsim[sel] <- dat[sel,paste0("Vsim.",WINDsite.2015)] + Uobs[sel] <- dat[sel,paste0("Uobs.",WINDsite.2015)] + Vobs[sel] <- dat[sel,paste0("Vobs.",WINDsite.2015)] + } # if(SITE=="HPL"){ + isNA <- is.na(Usim)|is.na(Vsim)|is.na(Uobs)|is.na(Vobs) + # a) U/V filter: both simulated U & V have to have the same sign as the observed (i.e., wind vector has to be same QUADRANT as observed) + # SEL <- sign(Usim)!=sign(Uobs) + # SEL <- sel|(sign(Vsim)!=sign(Vobs)) + + # b) U/V filter: filter out days based on +/-45o of observed (better than quadrant method above, since not throw way data close to quadrant boundaries + WSPD.sim <- sqrt(Usim^2 + Vsim^2); WSPD.obs <- sqrt(Uobs^2 + Vobs^2) + print("Observed windspeeds (quantile):") + print(signif(quantile(WSPD.obs,c(0,0.05,0.1,0.25,0.5,0.75,1.0),na.rm=T),3)) + # see: http://tornado.sfsu.edu/geosciences/classes/m430/Wind/WindDirection.html + WDIR.sim <- (180/pi)*atan2(y=-1*Vsim,x=-1*Usim) # wind direction (-180o to +180o); direction wind is blowing FROM + WDIR.obs <- (180/pi)*atan2(y=-1*Vobs,x=-1*Uobs) # wind direction (-180o to +180o); direction wind is blowing FROM + dWDIR <- WDIR.sim - WDIR.obs + sel <- abs(dWDIR) > 180 + dWDIR[sel&!isNA] <- sign(dWDIR[sel&!isNA])*(360 - abs(dWDIR[sel&!isNA])) + SEL <- abs(dWDIR) > 45 # simulated wind direction has to be within 45o of observed + SEL[WSPD.obs < 1.0] <- FALSE # not remove data if windspeed is too low (wind anemometer not reliable) + + print(paste("Out of",sum(!isNA),"non-NA U/V times,",sum(SEL[!isNA]),"failed to pass U/V filter: ",signif(100*sum(SEL[!isNA])/sum(!isNA),3),"%")) + # assign NAs to footsum and foot.basin when failed to pass U/V filter + print(paste("Before filtering, footsum has",sum(!is.na(dat$footsum)),"non-NA values")) + dat$footsum[SEL] <- NA + print(paste(" After filtering, footsum has",sum(!is.na(dat$footsum)),"non-NA values")) + dat$foot.basin[SEL] <- NA + dat.all <- dat +} # if(filterUV.TF){ + +######################################################### +#V. Calculate daily averages +dat <- dat.all +Time.day<-format(dat$Time,format="%Y%m%d") +dCH4.ave<-tapply(dat$dCH4,Time.day,mean,na.rm=T) +Time<-strptime(names(dCH4.ave),"%Y%m%d",tz="GMT") +footsum<-tapply(dat$footsum,Time.day,mean,na.rm=T) +foot.basin<-tapply(dat$foot.basin,Time.day,mean,na.rm=T) + +result<-data.frame(YYYYMMDD=names(dCH4.ave),Time,dCH4.ave,footsum,foot.basin,stringsAsFactors=FALSE) + +######################################################### +#VI. Divide CH4 enhancement by footprint strength to calculate CH4 fluxes +# e.g., F_CH4 = (CH4_HPL - CH4_FRU)/sum(foot_i,j) +Fch4.sum<-result[,"dCH4.ave"]/result[,"footsum"] #Fch4 calculated with total footprint +sel<-result[,"footsum"][umole/hr] +Ech4.sum<-Ech4.sum*(12.0111+4*1.008)*1E-6/1000/(1000) #[umole/hr]=>[10^3 kg/hr] +Ech4.basin<-Fch4.basin*AREA*3600 #[umole/m2/s]=>[umole/hr] +Ech4.basin<-Ech4.basin*(12.0111+4*1.008)*1E-6/1000/(1000) #[umole/hr]=>[10^3 kg/hr] +Ech4.basin.sd<-Fch4.basin.sd*AREA*3600 #[umole/m2/s]=>[umole/hr] +Ech4.basin.sd<-Ech4.basin.sd*(12.0111+4*1.008)*1E-6/1000/(1000) #[umole/hr]=>[10^3 kg/hr] + +######################################################### +#IX. Plot all months +dev.new();par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3) +XMAIN <- paste0(SITE,"\nfillFRU.TF=",fillFRU.TF,"; filterUV.TF=",filterUV.TF) +plot(frYr,Ech4.sum,main=XMAIN,pch=16,type="o",xlab="Year",ylab="Emissions of CH4 from Uintah Basin [10^3 kg/hr]",ylim=c(0,100),sub=XSUB) +lines(frYr,Ech4.basin,type="o",pch=16,col="blue") +#abline(h=55,col="red",lty=2) #Uintah Basin-wide CH4 emissions [10^3 kg/hr], as reported by Karion et al. [2013] +legend("topright",c("sum(foot)","sum(foot.basin)"),col=c("black","blue"),text.col=c("black","blue"),pch=16,cex=1.3,lty=1) +dev.copy(png,"Fch4_simple_0.png");dev.off() + +#Generate plots used for NCGG-8 presentation +Ech4.basin.seas<-Fch4.basin.seas*AREA*3600 #[umole/m2/s]=>[umole/hr] +Ech4.basin.seas<-Ech4.basin.seas*(12.0111+4*1.008)*1E-6/1000/(1000) #[umole/hr]=>[10^3 kg/hr] +dev.new();par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3,mar=c(5,5,4,5)) +XMAIN <- paste0(SITE,"; UThrs: ",paste(HRs,collapse=","),"\nfillFRU.TF=",fillFRU.TF,"; filterUV.TF=",filterUV.TF) +plot(as.numeric(names(Ech4.basin.seas)),Ech4.basin.seas,main=XMAIN,sub=XSUB, + pch=16,type="o",xlab="Month",ylab="Emissions of CH4 from Uintah Basin [10^3 kg/hr]") +#abline(h=55,col="red",lty=2,lwd=2) #Uintah Basin-wide CH4 emissions [10^3 kg/hr], as reported by Karion et al. [2013] +par(new=T) +plot(as.numeric(names(Fch4.basin.seas)),Fch4.basin.seas,xlab="",ylab="",pch=16,type="o",axes=F) +axis(4) +mtext("Average CH4 Flux [umole/m2/s]",side=4,cex=1.3,line=2) + +######################################################### +#X. Look at time series of Basin-Wide Emissions only for a subset of months +YYYYMM<-names(Ech4.basin) +frYr<-as.numeric(substring(YYYYMM,1,4))+(as.numeric(substring(YYYYMM,5,6))-0.5)/12 +SEL<-as.numeric(substring(YYYYMM,5,6))%in%MONsub +#calculate error bars, sample numbers +dat<-readRDS(resultname) +isNA<-is.na(dat[,"Fch4.basin"]); dat2<-dat[!isNA,] +YYYYMM<-substring(dat2[,"YYYYMMDD"],1,6) +NN<-tapply(YYYYMM,YYYYMM,length) +Ech4.basin.stderr<-Ech4.basin.sd +Ech4.basin.stderr[names(Ech4.basin.stderr)]<-Ech4.basin.sd[names(Ech4.basin.stderr)]/sqrt(NN[names(Ech4.basin.stderr)]) +Ech4.basin[names(NN)][NN=min(frYr1) & frYr <=max(frYr1) +plot(frYr[SEL2],Oil[SEL2],pch=16,type="o",xlab="",ylab="",col="orange",axes=F,lwd=2) +axis(4,col="orange",col.axis="orange") +mtext(ylab,side=4,cex=1.3,line=2,col="orange") +dev.copy(png,"Fch4_simple_1.5.png");dev.off() + + +#Look at time series ONLY for a subset of months--plot DEVIATION from annual average emissinos +YYYYMM<-names(Ech4.basin) +frYr<-as.numeric(substring(YYYYMM,1,4))+(as.numeric(substring(YYYYMM,5,6))-0.5)/12 +SEL<-as.numeric(substring(YYYYMM,5,6))%in%MONsub +#calculate error bars, sample numbers +dat<-readRDS(resultname) +isNA<-is.na(dat[,"Fch4.basin"]); dat2<-dat[!isNA,] +YYYYMM<-substring(dat2[,"YYYYMMDD"],1,6) +NN<-tapply(YYYYMM,YYYYMM,length) +Ech4.basin.stderr<-Ech4.basin.sd +Ech4.basin.stderr[names(Ech4.basin.stderr)]<-Ech4.basin.sd[names(Ech4.basin.stderr)]/sqrt(NN[names(Ech4.basin.stderr)]) +Ech4.basin[names(NN)][NN=min(frYr1) & frYr <=max(frYr1) +plot(frYr[SEL2],Oil[SEL2],pch=16,type="o",xlab="",ylab="",col="orange",axes=F,lwd=2) +axis(4,col="orange",col.axis="orange") +mtext(ylab,side=4,cex=1.3,line=2,col="orange") +dev.copy(png,"Fch4_simple_1.6.png");dev.off() + R <- cor(dEch4.basin[SEL], Oil[names(dEch4.basin[SEL])], use="na.or.complete") + + + +######################################################### +#XI. Plot emissions as % of CH4 from Natural Gas PRODUCTION (monthly) +#convert units of natural gas & oil production +tmp<-readRDS("State_of_Utah_oil_gas_DATA/gas.oil.production_monthly.rds") +counties<-unique(tmp$County) +NatGas <- tapply(tmp[,"Natural.Gas..MCF."],list(tmp[,"Time"]),sum,na.rm=TRUE) #sum monthly values over counties +NatGas <- MCF2kg(MCF=NatGas,CH4.VOLFRAC=CH4.VOLFRAC) #[MCF Natural Gas per Month] => [kg CH4 per Month] +NatGas <- NatGas/1000 #[kg CH4 per Month] => [10^3 kg CH4 per Month] +NatGas <- NatGas/(24*(365/12)) # [10^3 kg CH4 per Month] => [10^3 kg CH4 per Hour] +ylab<-"Natural Gas [10^3 kg/hr]" +YYYYMM<-paste0(substring(names(NatGas),1,4),substring(names(NatGas),6,7)) +names(NatGas)<-YYYYMM +production<-NatGas[names(Ech4.basin)] +YYYYMM<-names(Ech4.basin) +Year<-as.numeric(substring(YYYYMM,1,4)) +MM<-as.numeric(substring(YYYYMM,5,6)) +percEprod<-100*Ech4.basin/production +SEL<-MM%in%MONsub +frYr<-Year+(MM-0.5)/12 +XMAIN <- paste0(SITE,"; UThrs: ",paste(HRs,collapse=","),"\nfillFRU.TF=",fillFRU.TF,"; filterUV.TF=",filterUV.TF) +XSUB <- paste0("mettype=",mettype,"; Mons: ",paste(MONsub,collapse=","),"; Nday.min=",Nday.min,"; CH4.VOLFRAC=",CH4.VOLFRAC) +dev.new();par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3,mar=c(5,5,4,5)) +plot(frYr[SEL],percEprod[SEL],main=XMAIN,sub=XSUB,pch=16, + xlab="Year",ylab="Emissions of CH4 from Uintah Basin [% Production]") +dev.copy(png,"Fch4_simple_2.png");dev.off() + +#Look at time series for a subset of years, and months +dat<-readRDS(resultname) +isNA<-is.na(dat[,"Fch4.basin"]); dat2<-dat[!isNA,] +dat<-dat2[as.numeric(substring(dat2[,"YYYYMMDD"],5,6))%in%MONsub,] +#plot time series of monthly averages +Time.mon<-format(dat[,"Time"],format="%Y%m") +Fch4.basin<-tapply(dat[,"Fch4.basin"],Time.mon,mean,na.rm=T) +Fch4.basin[is.nan(Fch4.basin)] <- NA +# calculate how many valid days there are in each month +Ntmp <- tapply(dat[!is.na(dat$Fch4.basin),"Fch4.basin"],Time.mon[!is.na(dat$Fch4.basin)],length) +N <- Fch4.basin; N[1:length(N)] <- NA; N[names(Ntmp)] <- Ntmp +Fch4.basin[N[10^3 kg/hr] +Fch4.yr<-tapply(Fch4.basin,YYYY,mean,na.rm=T) +Fch4.yr.sd<-tapply(Fch4.basin,YYYY,sd,na.rm=T) +Ech4.yr<-Fch4.yr*CONV +Ech4.yr.sd<-Fch4.yr.sd*CONV +Ech4.yr.stderr<-Ech4.yr.sd/sqrt(NN) +YYYY<-names(Ech4.yr) +frYr<-as.numeric(substring(YYYY,1,4))+0.5 +XMAIN <- paste0(SITE,"; UThrs: ",paste(HRs,collapse=","),"\nfillFRU.TF=",fillFRU.TF,"; filterUV.TF=",filterUV.TF) +XSUB <- paste0("mettype=",mettype,"; Mons: ",paste(MONsub,collapse=","),"; Nday.min=",Nday.min) +XSUB <- paste0(XSUB,"\n XLIMS=",round(XLIMS[1],2),",",round(XLIMS[2],2),"; YLIMS=",round(YLIMS[1],2),",",round(YLIMS[2],2)) +dev.new();par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3,mar=c(5,5,4,5)) +ylims <- c(18,60) +plot(frYr,Ech4.yr,main=XMAIN,sub=XSUB,pch=16,type="o",ylim=ylims, + xlab="Year",ylab="Emissions of CH4 from Uintah Basin [10^3 kg/hr]") +#abline(h=55,col="red",lty=2,lwd=2) #Uintah Basin-wide CH4 emissions [10^3 kg/hr], as reported by Karion et al. [2013] +segments(x0=frYr,y0=Ech4.yr-Ech4.yr.stderr,x1=frYr,y1=Ech4.yr+Ech4.yr.stderr) +dev.copy(png,"Fch4_simple_4.png");dev.off() + +#Plot emissions as % of CH4 from Natural Gas PRODUCTION (annual) +production<-readRDS("State_of_Utah_oil_gas_DATA/gas.oil.production_annual.rds") +rownames(production)<-production$Year +NatGas <- production[as.character(production$Year) %in% rownames(Ech4.yr),"Natural.Gas..MCF."] # [MCF NatGas/Year] +NatGas <- MCF2kg(MCF=NatGas,CH4.VOLFRAC=CH4.VOLFRAC) #[MCF NatGas/Year] => [kg CH4/Year] +NatGas <- NatGas/1000 #[kg CH4/Year] => [10^3 kg CH4/Year] +NatGas <- NatGas/(24*365) # [10^3 kg CH4/Year] => [10^3 kg CH4/Hour] +Ech4.yr.perc<-100*Ech4.yr/NatGas #emission as % of CH4 in natural gas production +Ech4.yr.perc.stderr<-100*Ech4.yr.stderr/NatGas #stderr as % of CH4 in natural gas production +XMAIN <- paste0(SITE,"; UThrs: ",paste(HRs,collapse=","),"\nfillFRU.TF=",fillFRU.TF,"; filterUV.TF=",filterUV.TF) +XSUB <- paste0("mettype=",mettype,"; Mons: ",paste(MONsub,collapse=","),"; Nday.min=",Nday.min,"; CH4.VOLFRAC=",CH4.VOLFRAC) +dev.new();par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3,mar=c(5,5,4,5)) +ylims<-c(0,10) +plot(frYr,Ech4.yr.perc,main=XMAIN,sub=XSUB,pch=16,type="o", + xlab="Year",ylab="Emissions of CH4 from Uintah Basin [% Production]",ylim=ylims) +segments(x0=frYr,y0=Ech4.yr.perc-Ech4.yr.perc.stderr,x1=frYr,y1=Ech4.yr.perc+Ech4.yr.perc.stderr) +dev.copy(png,"Fch4_simple_5.png");dev.off() + + +######################################################### +#XII. Plot emissions in ENERGY UNITS from Basin (as % of energy in CH4 emissions from Natural Gas + Oil PRODUCTION (annual)) +# This is called the Energy-Normalized Methane Average (ENMA) emissions in Robertson et al. (2017) ES&T paper +production <- readRDS("State_of_Utah_oil_gas_DATA/gas.oil.production_annual.rds") +Oil <- production[as.character(production$Year) %in% rownames(Ech4.yr),"Oil..BBLs."] +names(Oil) <- rownames(Ech4.yr) +Oil.BTU <- Oil * 5.8E6 # oil production in energy units [Barrels/Yr] => [BTU/Yr]; from Robertson et al. (2017) +NatGas <- production[as.character(production$Year) %in% rownames(Ech4.yr),"Natural.Gas..MCF."] # [MCF NatGas/Year] +NatGas.BTU <- NatGas * 1E6 # gas production in energy units [MCF Natural Gas/Year] => [BTU/Year] from Robertson et al. (2017) +# convert methane emissions back from [10^3 kg CH4 per Hour] => [MCF (thousands of cubic feet) natural gas/Year] +Egas <- data.frame(Ech4.yr,Ech4.yr.stderr) # [10^3 kg CH4/hr] +Egas <- Egas*1000*24*365 # [10^3 kg CH4/hr] => [kg CH4/Year] +CONV <- MCF2kg(MCF=1,CH4.VOLFRAC=CH4.VOLFRAC) #[MCF Natural Gas] => [kg CH4] conversion factor +Egas <- Egas/CONV # [kg CH4/Year] => [MCF Natural Gas/Year] +Egas.BTU <- Egas * 1E6 # gas emissions in energy units [MCF Natural Gas/Year] => [BTU/Year] from Robertson et al. (2017) +ENBA <- 100*Egas.BTU/(Oil.BTU+NatGas.BTU) +dev.new();par(cex.axis=1.3,cex.lab=1.2,cex.main=1.3,mar=c(5,5,4,5)) +ylims <- c(0,10) +#ylims <- NULL +plot(frYr,ENBA$Ech4.yr,main=XMAIN,sub=XSUB,pch=16,type="o", + xlab="Year",ylab="Energy-Normalized Methane Average (ENMA) emissions\n [% Production]",ylim=ylims) +segments(x0=frYr,y0=ENBA$Ech4.yr-ENBA$Ech4.yr.stderr,x1=frYr,y1=ENBA$Ech4.yr+ENBA$Ech4.yr.stderr) +dev.copy(png,"Fch4_simple_6.png");dev.off() + +######################################################### +#XIII. Combine ENMA emissions as % of total NatGas+Oil production of energy with CH4 emissions as % of NatGas production [both %] in same plot +YYYY<-names(Ech4.yr.perc) +frYr<-as.numeric(substring(YYYY,1,4))+0.5 +dev.new();par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3,mar=c(5,5,4,5)) +ylims<-c(0,10) +#ylims<-c(0,18) +#ylims<-NULL +XMAIN <- paste0(SITE,"; UThrs: ",paste(HRs,collapse=","),"\nfillFRU.TF=",fillFRU.TF,"; filterUV.TF=",filterUV.TF) +XSUB <- paste0("mettype=",mettype,"; Mons: ",paste(MONsub,collapse=","),"; Nday.min=",Nday.min,"; CH4.VOLFRAC=",CH4.VOLFRAC) +plot(frYr,Ech4.yr.perc,main=XMAIN,sub=XSUB,pch=16,type="o", + xlab="Year",ylab="CH4 Emissions [% NatGas Production or % Energy Produced]",ylim=ylims,lwd=2) +segments(x0=frYr,y0=Ech4.yr.perc-Ech4.yr.perc.stderr,x1=frYr,y1=Ech4.yr.perc+Ech4.yr.perc.stderr) +lines(frYr,ENBA$Ech4.yr,main=XMAIN,sub=XSUB,pch=16,type="o",col="darkgray",lwd=2) +segments(x0=frYr,y0=ENBA$Ech4.yr-ENBA$Ech4.yr.stderr,x1=frYr,y1=ENBA$Ech4.yr+ENBA$Ech4.yr.stderr,col="darkgray") +legend("bottomleft",c("% NatGas Production","% Energy (NatGas+Oil)"),col=c("black","darkgray"),pch=16,lwd=2,text.col=c("black","darkgray"),cex=1.2,bty="n") +dev.copy(png,"Fch4_simple_7.png");dev.off() + + +######################################################### +#XIV. Superimpose time series of monthly Basin-Wide Emissions with ANNUAL values +YYYYMM<-names(Ech4.basin) +frYr<-as.numeric(substring(YYYYMM,1,4))+(as.numeric(substring(YYYYMM,5,6))-0.5)/12 +SEL<-as.numeric(substring(YYYYMM,5,6))%in%MONsub +# calculate error bars, sample numbers +dat<-readRDS(resultname) +isNA<-is.na(dat[,"Fch4.basin"]); dat2<-dat[!isNA,] +YYYYMM<-substring(dat2[,"YYYYMMDD"],1,6) +NN<-tapply(YYYYMM,YYYYMM,length) +Ech4.basin.stderr<-Ech4.basin.sd +Ech4.basin.stderr[names(Ech4.basin.stderr)]<-Ech4.basin.sd[names(Ech4.basin.stderr)]/sqrt(NN[names(Ech4.basin.stderr)]) +Ech4.basin[names(NN)][NN [kg CH4 per Month] + NatGas <- NatGas/1000 #[kg CH4 per Month] => [10^3 kg CH4 per Month] + NatGas <- NatGas/(24*(365/12)) # [10^3 kg CH4 per Month] => [10^3 kg CH4 per Hour] + YYYYMM <- paste0(substring(names(NatGas),1,4),substring(names(NatGas),6,7)) + names(NatGas) <- YYYYMM + NatGas <- NatGas[names(Ech4.basin)] + YYYYMMsel <- paste0(rep(YEARs,each=length(MONsub)),rep(formatC(MONsub,width=2,flag="0"),length(YEARs))) + NatGas <- NatGas[YYYYMMsel] + # average production during selected months each year + NatGas.yr <- tapply(NatGas,substring(names(NatGas),1,4),mean,na.rm=T) # [10^3 kg CH4/Hr] + NatGas.yr.bycounty[[cc]] <- NatGas.yr + } # for(cc in 1:length(counties)){ + names(NatGas.yr.bycounty) <- counties + + #2) scale gas production in Duchesne county based on same emission characteristics as Uintah County using ratio in gas production + # Ech4.basin.yr.domainsens <- Ech4.basin.yr*(1+NatGas.yr.bycounty$DUCHESNE/NatGas.yr.bycounty$UINTAH) + + #3) calculate emissions from Duchesne county as solely due to leakage from gas production from oil wells (using leakage rate from oil wells of 14.86%) + Ech4.basin.yr.domainsens <- Ech4.basin.yr+NatGas.yr.bycounty$DUCHESNE*0.1486 + ############################################################################################# +lines(frYr,Ech4.basin.yr,pch=16,type="o",lwd=2,col="gray",lty=2) +lines(frYr,Ech4.basin.yr.domainsens,pch=16,type="o",lwd=3,col="gray") +segments(x0=frYr,y0=Ech4.basin.yr.domainsens-Ech4.basin.yr.stderr,x1=frYr,y1=Ech4.basin.yr.domainsens+Ech4.basin.yr.stderr,col="gray") +# add in default emissions (without modifying model domain--output from "Fch4_simple_multimettype.r" +default <- readRDS(paste0(mettype,".rds")) +Ech4.yr <- default$Ech4.yr; YYYY<-names(Ech4.yr) +frYr<-as.numeric(substring(YYYY,1,4))+0.5 +lines(frYr,Ech4.yr,lwd=2,type="o",pch=16) +legend("topright",c("Base case(Uintah&Duchesne)","Domain(Uintah)","Domain(Uintah) + Duchesne"), + lty=c(1,2,1),lwd=c(3,2,2),col=c("black","gray","gray"),text.col=c("black","gray","gray"),cex=1.1) +dev.copy(png,"Fch4_simple_8.5.png");dev.off() +} # if(domainsens.TF){ + +######################################################### +#XV. Plot emissions as % of CH4 from Natural Gas PRODUCTION during SELECTED MONTHS, instead of using ANNUAL PRODUCTION values +#convert units of natural gas & oil production +tmp<-readRDS("State_of_Utah_oil_gas_DATA/gas.oil.production_monthly.rds") +counties<-unique(tmp$County) +NatGas <- tapply(tmp[,"Natural.Gas..MCF."],list(tmp[,"Time"]),sum,na.rm=TRUE) #sum monthly values over counties +NatGas <- MCF2kg(MCF=NatGas,CH4.VOLFRAC=CH4.VOLFRAC) #[MCF Natural Gas per Month] => [kg CH4 per Month] +NatGas <- NatGas/1000 #[kg CH4 per Month] => [10^3 kg CH4 per Month] +NatGas <- NatGas/(24*(365/12)) # [10^3 kg CH4 per Month] => [10^3 kg CH4 per Hour] +YYYYMM <- paste0(substring(names(NatGas),1,4),substring(names(NatGas),6,7)) +names(NatGas) <- YYYYMM +NatGas <- NatGas[names(Ech4.basin)] +YYYYMMsel <- paste0(rep(YEARs,each=length(MONsub)),rep(formatC(MONsub,width=2,flag="0"),length(YEARs))) +NatGas <- NatGas[YYYYMMsel] +# average production during selected months each year +NatGas.yr <- tapply(NatGas,substring(names(NatGas),1,4),mean,na.rm=T) # [10^3 kg CH4/Hr] +frYr<-as.numeric(names(NatGas.yr))+0.5 +# Ech4.yr is in units of [10^3 kg CH4/hr] +Ech4.yr.perc <- Ech4.yr*100/NatGas.yr +XMAIN <- paste0(SITE,"; UThrs: ",paste(HRs,collapse=","),"\nfillFRU.TF=",fillFRU.TF,"; filterUV.TF=",filterUV.TF) +XSUB <- paste0("mettype=",mettype,"; Mons: ",paste(MONsub,collapse=","),"; Nday.min=",Nday.min,"; CH4.VOLFRAC=",CH4.VOLFRAC) +dev.new();par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3,mar=c(5,5,4,5)) +ylims<-c(1,10) +plot(frYr,Ech4.yr.perc,main=XMAIN,sub=XSUB,pch=16,type="o", + xlab="Year",ylab="Emissions of CH4 from Uintah Basin [% Production]",ylim=ylims) +segments(x0=frYr,y0=Ech4.yr.perc-Ech4.yr.perc.stderr,x1=frYr,y1=Ech4.yr.perc+Ech4.yr.perc.stderr) +dev.copy(png,"Fch4_simple_9.png");dev.off() + +######################################################### +#XVI. Combine ENMA emissions as % of total NatGas+Oil production of energy with CH4 emissions as % of NatGas production [both %] in same plot +tmp<-readRDS("State_of_Utah_oil_gas_DATA/gas.oil.production_monthly.rds") +counties<-unique(tmp$County) +Oil <- tapply(tmp[,"Oil..BBLs."],list(tmp[,"Time"]),sum,na.rm=TRUE) #sum monthly values over counties +Oil.BTU <- Oil * 5.8E6 # oil production in energy units [Barrels/Month] => [BTU/Month]; from Robertson et al. (2017) +NatGas <- tapply(tmp[,"Natural.Gas..MCF."],list(tmp[,"Time"]),sum,na.rm=TRUE) #sum monthly values over counties +NatGas.BTU <- NatGas*1E6 # gas productioon in energy units [MCF Natural Gas/Month] => {BTU/Month]; from Robertson et al. (2017) +# select subset based on Year/Month +YYYYMMsel <- paste0(rep(YEARs,each=length(MONsub)),"-",rep(formatC(MONsub,width=2,flag="0"),length(YEARs)),"-01") +NatGas.BTU <- NatGas.BTU[YYYYMMsel] ; Oil.BTU <- Oil.BTU[YYYYMMsel] +# average each year +NatGas.BTU.monave <- tapply(NatGas.BTU,substring(names(NatGas.BTU),1,4),mean,na.rm=T) +Oil.BTU.monave <- tapply(Oil.BTU,substring(names(Oil.BTU),1,4),mean,na.rm=T) +NatGas.BTU.hr <- NatGas.BTU.monave/(24*(365/12)) # [BTU per Month] => [BTU per Hour] +Oil.BTU.hr <- Oil.BTU.monave/(24*(365/12)) # [BTU per Month] => [BTU per Hour] +# convert methane emissions back from [10^3 kg CH4 per Hour] => [MCF (thousands of cubic feet) natural gas/Year] +Egas <- data.frame(Ech4.yr,Ech4.yr.stderr) # [10^3 kg CH4/hr] +Egas <- Egas*1000 # [10^3 kg CH4/hr] => [kg CH4/hr] +CONV <- MCF2kg(MCF=1,CH4.VOLFRAC=CH4.VOLFRAC) #[MCF Natural Gas] => [kg CH4] conversion factor +Egas <- Egas/CONV # [kg CH4/hr] => [MCF Natural Gas/hr] +Egas <- Egas*CH4.VOLFRAC # [MCF Natural Gas/hr] => [MCF methane/hr] +Egas.BTU.hr <- Egas * 1E6 # gas emissions in energy units [MCF Natural Gas/hr] OR [MCF methane/hr]=> [BTU/hr] from Robertson et al. (2017) +ENBA <- 100*Egas.BTU.hr/(Oil.BTU.hr+NatGas.BTU.hr) + +YYYY<-names(Ech4.yr.perc) +frYr<-as.numeric(substring(YYYY,1,4))+0.5 +dev.new();par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3,mar=c(5,5,4,5)) +ylims<-c(1,10) +#ylims<-NULL +XMAIN <- paste0(SITE,"; UThrs: ",paste(HRs,collapse=","),"\nfillFRU.TF=",fillFRU.TF,"; filterUV.TF=",filterUV.TF) +XSUB <- paste0("mettype=",mettype,"; Mons: ",paste(MONsub,collapse=","),"; Nday.min=",Nday.min,"; CH4.VOLFRAC=",CH4.VOLFRAC) +plot(frYr,Ech4.yr.perc,main=XMAIN,sub=XSUB,pch=16,type="o", + xlab="Year",ylab="CH4 Emissions [% NatGas Production or % Energy Produced]",ylim=ylims,lwd=2) +segments(x0=frYr,y0=Ech4.yr.perc-Ech4.yr.perc.stderr,x1=frYr,y1=Ech4.yr.perc+Ech4.yr.perc.stderr) +lines(frYr,ENBA$Ech4.yr,main=XMAIN,sub=XSUB,pch=16,type="o",col="darkgray",lwd=2) +segments(x0=frYr,y0=ENBA$Ech4.yr-ENBA$Ech4.yr.stderr,x1=frYr,y1=ENBA$Ech4.yr+ENBA$Ech4.yr.stderr,col="darkgray") +legend("bottomleft",c("% NatGas Production","% Energy (NatGas+Oil)"),col=c("black","darkgray"),pch=16,lwd=2,text.col=c("black","darkgray"),cex=1.2,bty="n") +dev.copy(png,"Fch4_simple_10.png");dev.off() + + diff --git a/Fch4_simple_multimettype.r b/Fch4_simple_multimettype.r new file mode 100644 index 0000000..41ae8da --- /dev/null +++ b/Fch4_simple_multimettype.r @@ -0,0 +1,288 @@ +# Simple STILT analyses of Uinta Basin methane emissions +# Idea is to calculate AFTERNOON enhancements (CH4_HPL - CH4_FRU), & then divide by total footprint strength: +# F_CH4 = (CH4_HPL - CH4_FRU)/sum(foot_i,j) +# +# Based on 'Fch4_simpleV8.r', but plot results from MULTIPLE MET FIILES driving STILT +# V2(210522): for paper submission, add Karion number, change oil/gas colors, and remove some of the infoormational labels +# May 3rd, 2021 by John C. Lin (John.Lin@utah.edu) + +require(ncdf4); require(fields); require(geosphere) +#################### +YEARs <- 2015:2020 +MONsub <- 4:9 #subset of months to calculate fluxes and leak rates +HRs<-20:23 #[UTC] only analyze afternoon hours, following Foster papers +#HRs<-17:23 #[UTC] only analyze afternoon hours, following Foster papers +SITE<-"HPL" #HPL is site to focus on for calculating long-term F_CH4 +#SITE<-"ROO" +#SITE<-"CSP" +if(SITE=="CSP")YEARs <- c(2016,2019,2020) +if(SITE=="ROO")YEARs <- 2015:2019 + +mettypes <- c("HRRR","NAM12","WRF27") + +XLIMS<-c(-110.6,-109) # lon range of Uinta Basin +YLIMS<-c(39.9,40.5) # lat range of Uinta Basin + + +Nday.min <- 10 # minimum days necessary for a monthly average to be retained (otherwise assigned NA) +fillFRU.TF <- TRUE # fill in gaps in background (FRU) time series? +filterUV.TF <- TRUE # filter times based on U/V (filter out times when HRRR is off--i.e., large transport errors) ? +CH4.VOLFRAC <- 0.89 # volume fraction of CH4 in natural gas of 0.89 [Karion et al. 2013] + +col.gas <- "#B91C1C" # red +col.oil <- "#1D4ED8" # blue +col.tot <- "black" # total/black +#################### + + +if(SITE=="HPL"){LAT <- 40.1434; LON <- -109.4680; zagl <- 4.06; WINDsite <- "UBHSP"; WINDsite.2015 <- "A1633"} # UBHSP site not available in 2015, so use A1633 (RedWash) site instead in 2015 +if(SITE=="ROO"){LAT <- 40.2941; LON <- -110.0090; zagl <- 4.06; WINDsite <- "QRS"} +if(SITE=="CSP"){LAT <- 40.0509; LON <- -110.0194; zagl <- 4.06; WINDsite <- "UBCSP"} + +# Function that converts from natural gas volume in thousands of [MCF] to [kg of CH4] +MCF2kg <- function(MCF,CH4.VOLFRAC=0.89){ +NatGas<-MCF*CH4.VOLFRAC #[MCF natural gas] => [MCF CH4], using volume fraction of CH4 in natural gas of 0.89 [Karion et al. 2013] +XLIMS<-c(-110.6,-109) # lon range of Uinta Basin +YLIMS<-c(39.9,40.5) # lat range of Uinta Basin + +NatGas<-NatGas*1000 #[MCF CH4] => [ft^3 CH4] + NatGas<-NatGas*0.0283 #[ft^3 CH4] => [m^3 CH4] + R.CH4<-8.3143*1000/16.043 #Ideal Gas constant of CH4 [J/K/kg]: (Rg/mCH4), where Rg is universal gas constant and mCH4 is molar mass of CH4 + rho.CH4<-(101300)/(R.CH4*288.7) #density of CH4 [kg/m3], using P=101.3kPa and T=288.7K [Karion et al. 2013] + NatGas<-NatGas*rho.CH4 #[m^3 CH4] => [kg CH4] + return(NatGas) +} # MCF2kg <- function(MCF,CH4.VOLFRAC=0.89){ + + +######################################################### +# Calculate average gridcell area +dx1<-distHaversine(c(XLIMS[1],YLIMS[1]),c(XLIMS[2],YLIMS[1])) #[m] +dx2<-distHaversine(c(XLIMS[1],YLIMS[2]),c(XLIMS[2],YLIMS[2])) #[m] +dy1<-distHaversine(c(XLIMS[1],YLIMS[1]),c(XLIMS[1],YLIMS[2])) #[m] +dy2<-distHaversine(c(XLIMS[2],YLIMS[1]),c(XLIMS[2],YLIMS[2])) #[m] +AREA<-mean(c(dx1,dx2))*mean(c(dy1,dy2)) #[m^2] + +for(i in 1:length(mettypes)){ + #mettype <- "HRRR" + mettype <- mettypes[i] + resultname<-paste("Fch4_",SITE,"_daily_",mettype,".rds",sep="") + +######################################################### +# Read in data generated by "Fch4_simple.r" and calculate monthly-averaged emissions +dat<-readRDS(resultname) +#plot time series of monthly averages +Time.mon<-format(dat[,"Time"],format="%Y%m") +Fch4.sum<-tapply(dat[,"Fch4.sum"],Time.mon,mean,na.rm=T) +Fch4.basin<-tapply(dat[,"Fch4.basin"],Time.mon,mean,na.rm=T) +Fch4.sum[is.nan(Fch4.sum)] <- NA; Fch4.basin[is.nan(Fch4.basin)] <- NA +# calculate how many valid days there are in each month +Ntmp <- tapply(dat[!is.na(dat$Fch4.basin),"Fch4.basin"],Time.mon[!is.na(dat$Fch4.basin)],length) +N <- Fch4.basin; N[1:length(N)] <- NA; N[names(Ntmp)] <- Ntmp +Fch4.basin[N[umole/hr] +Ech4.sum<-Ech4.sum*(12.0111+4*1.008)*1E-6/1000/(1000) #[umole/hr]=>[10^3 kg/hr] +Ech4.basin<-Fch4.basin*AREA*3600 #[umole/m2/s]=>[umole/hr] +Ech4.basin<-Ech4.basin*(12.0111+4*1.008)*1E-6/1000/(1000) #[umole/hr]=>[10^3 kg/hr] +Ech4.basin.sd<-Fch4.basin.sd*AREA*3600 #[umole/m2/s]=>[umole/hr] +Ech4.basin.sd<-Ech4.basin.sd*(12.0111+4*1.008)*1E-6/1000/(1000) #[umole/hr]=>[10^3 kg/hr] + + + +######################################################### +# Calculate annual emissions +dat<-readRDS(resultname) +isNA<-is.na(dat[,"Fch4.basin"]); dat2<-dat[!isNA,] +dat<-dat2[as.numeric(substring(dat2[,"YYYYMMDD"],5,6))%in%MONsub,] +#plot time series of monthly averages +Time.mon<-format(dat[,"Time"],format="%Y%m") +Fch4.basin<-tapply(dat[,"Fch4.basin"],Time.mon,mean,na.rm=T) +Fch4.basin[is.nan(Fch4.basin)] <- NA +# calculate how many valid days there are in each month +Ntmp <- tapply(dat[!is.na(dat$Fch4.basin),"Fch4.basin"],Time.mon[!is.na(dat$Fch4.basin)],length) +N <- Fch4.basin; N[1:length(N)] <- NA; N[names(Ntmp)] <- Ntmp +Fch4.basin[N[10^3 kg/hr] +Fch4.yr<-tapply(Fch4.basin,YYYY,mean,na.rm=T) +Fch4.yr.sd<-tapply(Fch4.basin,YYYY,sd,na.rm=T) +Ech4.yr<-Fch4.yr*CONV +Ech4.yr.sd<-Fch4.yr.sd*CONV +Ech4.yr.stderr<-Ech4.yr.sd/sqrt(NN) + + +######################################################### +# Calculatet monthly Basin-Wide Emissions +YYYYMM<-names(Ech4.basin) +frYr<-as.numeric(substring(YYYYMM,1,4))+(as.numeric(substring(YYYYMM,5,6))-0.5)/12 +SEL<-as.numeric(substring(YYYYMM,5,6))%in%MONsub +# calculate error bars, sample numbers +dat<-readRDS(resultname) +isNA<-is.na(dat[,"Fch4.basin"]); dat2<-dat[!isNA,] +YYYYMM<-substring(dat2[,"YYYYMMDD"],1,6) +NN<-tapply(YYYYMM,YYYYMM,length) +Ech4.basin.stderr<-Ech4.basin.sd +Ech4.basin.stderr[names(Ech4.basin.stderr)]<-Ech4.basin.sd[names(Ech4.basin.stderr)]/sqrt(NN[names(Ech4.basin.stderr)]) +Ech4.basin[names(NN)][NN [kg CH4 per Month] +NatGas <- NatGas/1000 #[kg CH4 per Month] => [10^3 kg CH4 per Month] +NatGas <- NatGas/(24*(365/12)) # [10^3 kg CH4 per Month] => [10^3 kg CH4 per Hour] +YYYYMM <- paste0(substring(names(NatGas),1,4),substring(names(NatGas),6,7)) +names(NatGas) <- YYYYMM +NatGas <- NatGas[names(Ech4.basin)] +YYYYMMsel <- paste0(rep(YEARs,each=length(MONsub)),rep(formatC(MONsub,width=2,flag="0"),length(YEARs))) +NatGas <- NatGas[YYYYMMsel] +# average production during selected months each year +NatGas.yr <- tapply(NatGas,substring(names(NatGas),1,4),mean,na.rm=T) # [10^3 kg CH4/Hr] +frYr<-as.numeric(names(NatGas.yr))+0.5 +# Ech4.yr is in units of [10^3 kg CH4/hr] +Ech4.yr.perc <- Ech4.yr*100/NatGas.yr +Ech4.yr.perc.stderr<-100*Ech4.yr.stderr/NatGas.yr #stderr as % of CH4 in natural gas production + + +#XMAIN <- paste0(SITE,"; UThrs: ",paste(HRs,collapse=","),"\nfillFRU.TF=",fillFRU.TF,"; filterUV.TF=",filterUV.TF) +#XSUB <- paste0("mettype=",mettype,"; Mons: ",paste(MONsub,collapse=","),"; Nday.min=",Nday.min,"; CH4.VOLFRAC=",CH4.VOLFRAC) +#dev.new();par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3,mar=c(5,5,4,5)) +#ylims<-c(0,10) +#plot(frYr,Ech4.yr.perc,main=XMAIN,sub=XSUB,pch=16,type="o", +# xlab="Year",ylab="Emissions of CH4 from Uinta Basin [% Production]",ylim=ylims) +#segments(x0=frYr,y0=Ech4.yr.perc-Ech4.yr.perc.stderr,x1=frYr,y1=Ech4.yr.perc+Ech4.yr.perc.stderr) +#dev.copy(png,"Fch4_simple_9.png");dev.off() + +######################################################### +# Calculate Energy-Normalized Methane Average (ENMA) emissions--i.e., CH4 emissions as % of total NatGas+Oil production of energy +tmp<-readRDS("State_of_Utah_oil_gas_DATA/gas.oil.production_monthly.rds") +counties<-unique(tmp$County) +Oil <- tapply(tmp[,"Oil..BBLs."],list(tmp[,"Time"]),sum,na.rm=TRUE) #sum monthly values over counties +Oil.BTU <- Oil * 5.8E6 # oil production in energy units [Barrels/Month] => [BTU/Month]; from Robertson et al. (2017) +NatGas <- tapply(tmp[,"Natural.Gas..MCF."],list(tmp[,"Time"]),sum,na.rm=TRUE) #sum monthly values over counties +NatGas.BTU <- NatGas*1E6 # gas productioon in energy units [MCF Natural Gas/Month] => {BTU/Month]; from Robertson et al. (2017) +# select subset based on Year/Month +YYYYMMsel <- paste0(rep(YEARs,each=length(MONsub)),"-",rep(formatC(MONsub,width=2,flag="0"),length(YEARs)),"-01") +NatGas.BTU <- NatGas.BTU[YYYYMMsel] ; Oil.BTU <- Oil.BTU[YYYYMMsel] +# average each year +NatGas.BTU.monave <- tapply(NatGas.BTU,substring(names(NatGas.BTU),1,4),mean,na.rm=T) +Oil.BTU.monave <- tapply(Oil.BTU,substring(names(Oil.BTU),1,4),mean,na.rm=T) +NatGas.BTU.hr <- NatGas.BTU.monave/(24*(365/12)) # [BTU per Month] => [BTU per Hour] +Oil.BTU.hr <- Oil.BTU.monave/(24*(365/12)) # [BTU per Month] => [BTU per Hour] +# convert methane emissions back from [10^3 kg CH4 per Hour] => [MCF (thousands of cubic feet) natural gas/Year] +Egas <- data.frame(Ech4.yr,Ech4.yr.stderr) # [10^3 kg CH4/hr] +Egas <- Egas*1000 # [10^3 kg CH4/hr] => [kg CH4/hr] +CONV <- MCF2kg(MCF=1,CH4.VOLFRAC=CH4.VOLFRAC) #[MCF Natural Gas] => [kg CH4] conversion factor +Egas <- Egas/CONV # [kg CH4/hr] => [MCF Natural Gas/hr] +Egas <- Egas*CH4.VOLFRAC # [MCF Natural Gas/hr] => [MCF methane/hr] +Egas.BTU.hr <- Egas * 1E6 # gas emissions in energy units [MCF Natural Gas/hr] OR [MCF methane/hr]=> [BTU/hr] from Robertson et al. (2017) +ENBA <- 100*Egas.BTU.hr/(Oil.BTU.hr+NatGas.BTU.hr) + +######################################################### +# Save all of the relevant objects as a LIST for each met type +result <- list(Ech4.basin=Ech4.basin,Ech4.basin.sd=Ech4.basin.sd,Ech4.basin.stderr=Ech4.basin.stderr, + Ech4.yr=Ech4.yr,Ech4.yr.sd=Ech4.yr.sd,Echr.yr.stderr=Ech4.yr.stderr, + Ech4.yr.perc=Ech4.yr.perc,Ech4.yr.perc.stderr,ENBA=ENBA) +finalresultfile <- paste0(mettype,".rds") +saveRDS(result,file=finalresultfile) +print(paste(finalresultfile,"generated")) + +} # for(i in 1:length(mettypes)){ + + +######################################################### +# Plot annual emissions +dev.new();par(cex.axis=1.5,cex.lab=1.5,cex.main=1.6,mar=c(5,5,4,3)) +#XMAIN <- paste0("CH4 Emissions from Uinta Basin\nUThrs: ",paste(HRs,collapse=",")) +#XSUB <- paste0("Mons: ",paste(MONsub,collapse=","),"; Nday.min=",Nday.min) +XMAIN <- expression("CH"[4]*" Emissions from Uinta Basin") +xlims <- c(2012,2021) +XSUB <- "" +plot(0,0,main=XMAIN,sub=XSUB,pch=16,type="o",xlim=xlims, + ylim=c(18,60),xlab="Year",ylab=expression(paste("CH"[4]," Emissions [10"^3," kg hr"^-1,"]"))) +Ech4.yr.sum <- 0 +for(mm in 1:length(mettypes)){ + mettype <- mettypes[mm] + dat <- readRDS(paste0(mettype,".rds")) + Ech4.yr <- dat$Ech4.yr + YYYY<-names(Ech4.yr) + frYr<-as.numeric(substring(YYYY,1,4))+0.5 + lines(frYr,Ech4.yr,lty=mm+1,lwd=2) + #segments(x0=frYr,y0=Ech4.yr-Ech4.yr.stderr,x1=frYr,y1=Ech4.yr+Ech4.yr.stderr) + Ech4.yr.sum <- Ech4.yr.sum + Ech4.yr +} # for(mm in 1:length(mettypes)){ +axis(1,at=seq(xlims[1],xlims[2],1),cex=1.5,labels=FALSE) +Ech4.yr.ave <- Ech4.yr.sum/length(mettypes) +lines(frYr,Ech4.yr.ave,lty=1,lwd=4) +# add Karion et al. [2013] emissions result based on Feb. 2012 observations +x.Karion <- 2012 + 1.5/12; y.Karion <- 55; y.sigma <- 15 +points(x.Karion,y.Karion,pch=17,cex=1.5) +arrows(x0=x.Karion,y0=y.Karion-y.sigma,x1=x.Karion,y1=y.Karion+y.sigma,code=3,angle=90,length=0.1,lwd=2) +text(x=2013.5,y=y.Karion,"Karion et al.\n(2013)",cex=1.5) +legend("topright",c("Model Mean",mettypes),lty=1:(length(mettypes)+1),lwd=c(4,rep(2,length(mettypes))),cex=1.3) +dev.copy(png,"Fch4_annual_ALLmet.png");dev.off() + + +dev.new();par(cex.axis=1.5,cex.lab=1.5,cex.main=1.5,mar=c(5,6,5,3)) +ylims <- c(1,10) +xlims <- c(2012,2021) +#XMAIN <- paste0("CH4 Emissions as % of Natural Gas Production or\n Total Energy Production from Uinta Basin\nUThrs: ",paste(HRs,collapse=",")) +#XSUB <- paste0("Mons: ",paste(MONsub,collapse=","),"; Nday.min=",Nday.min) +XMAIN <- expression(paste("CH"[4]," Emissions as % of Production from Uinta Basin")) +XSUB <- "" +plot(0,0,main=XMAIN,sub=XSUB,pch=16,, + xlab="Year",ylab="% NatGas Production or\n% Energy Production",xlim=xlims,ylim=ylims) +Ech4.yr.perc.sum <- 0 +ENBA.sum <- 0 +for(mm in 1:length(mettypes)){ + mettype <- mettypes[mm] + dat <- readRDS(paste0(mettype,".rds")) + Ech4.yr.perc <- dat$Ech4.yr.perc + print(paste("----------",mettype,"%NatGas Prod----------")) + print(paste(round(mean(Ech4.yr.perc),2),"+/-",round(sd(Ech4.yr.perc),2),"[%]")) + ENBA <- dat$ENBA + YYYY<-names(Ech4.yr.perc) + frYr<-as.numeric(substring(YYYY,1,4))+0.5 + lines(frYr,Ech4.yr.perc,lty=mm+1,lwd=2,col=col.gas) + #segments(x0=frYr,y0=Ech4.yr.perc-Ech4.yr.perc.stderr,x1=frYr,y1=Ech4.yr.perc+Ech4.yr.perc.stderr) + lines(frYr,ENBA$Ech4.yr,,col=col.tot,lty=mm+1,lwd=2) + #segments(x0=frYr,y0=ENBA$Ech4.yr-ENBA$Ech4.yr.stderr,x1=frYr,y1=ENBA$Ech4.yr+ENBA$Ech4.yr.stderr,col="darkgray") + Ech4.yr.perc.sum <- Ech4.yr.perc.sum + Ech4.yr.perc + ENBA.sum <- ENBA.sum + ENBA$Ech4.yr +} # for(mm in 1:length(mettypes)){ +axis(1,at=seq(xlims[1],xlims[2],1),cex=1.5,labels=FALSE) +Ech4.yr.perc.ave <- Ech4.yr.perc.sum/length(mettypes) +ENBA.ave <- ENBA.sum/length(mettypes) +lines(frYr,Ech4.yr.perc.ave,,col=col.gas,lty=1,lwd=4) +lines(frYr,ENBA.ave,,col=col.tot,lty=1,lwd=4) +# add Karion et al. [2013] leakage rate result based on Feb. 2012 observations +x.Karion <- 2012 + 1.5/12; y.Karion <- mean(c(6.2,11.7)); y.sigma <- y.Karion-6.2 +points(x.Karion,y.Karion,pch=17,cex=1.5,col=col.gas) +arrows(x0=x.Karion,y0=y.Karion-y.sigma,x1=x.Karion,y1=y.Karion+y.sigma,code=3,angle=90,length=0.1,lwd=2,col=col.gas) +text(x=2013.5,y=y.Karion,"Karion et al.\n(2013)",cex=1.5,col=col.gas) +legend("bottomleft",c("Model Mean",mettypes),lty=1:(length(mettypes)+1),lwd=c(4,rep(2,length(mettypes))),cex=1.3,bty="n") +legend(x=2015.5,y=3,legend=c("% NatGas Production","% Energy Production\n(NatGas+Oil)"), + col=c(col.gas,col.tot),lty=1,lwd=2,text.col=c(col.gas,col.tot),cex=1.3,bty="n") +dev.copy(png,"Fch4_percGasOil_ALLmet.png");dev.off() + + diff --git a/README b/README new file mode 100644 index 0000000..6cfcd3b --- /dev/null +++ b/README @@ -0,0 +1,82 @@ +March 2021 +0) plot_CH4tseries_diurnal_seasonal.r +1) Fch4_simple.r +1.5) Fch4_simple_multimettype.r +1.6) export_Fch4simple.r +2) Trend_foot_dCH4.r +3) average_foot_monthly.r +4) plot_foot_and_wells.r +5) grid_well_info_foot.r +6) convolve_foot_wellinfo.r +7) plot_foot_wellinfo_convolved.r +8) merge_CH4_UVwind_wellinfo.r +9) dCH4_vs_foot_wellinfo.r +10) plot_foot_wellinfo_tseries_trend.r + +---------------------------------- +May 28th, 2019 + +Simple STILT analyses of Uintah Basin methane emissions +Idea is to calculate AFTERNOON enhancements (CH4_HPL - CH4_FRU), and then divide by total footprint strength: + F_CH4 = (CH4_HPL - CH4_FRU)/sum(foot_i,j) + +Then see how stable F_CH4 is over different seasons and years. + +"Fch4_simple.r" + + +-------- Forwarded Message -------- +Subject: Re: Uintah Basin CH4 inversion time period +Date: Tue, 12 Feb 2019 16:41:32 -0700 +From: John C. Lin +Reply-To: John.Lin@utah.edu +To: Ryan Bares +CC: Ben Fasoli , Lewis Kunik + + +Thanks Ryan! Let's focus on March~May 2016 for the Uintah Basin analyses then. + +John + +On 2019-02-12 4:25 p.m., Ryan Bares wrote: +> See my comments below regarding each site. +> +> FRU: A small 3 day data gap at the beginning of March. Other than that the data looks good. +> ROO: Two small 1 to 2 day gaps. The rest of the data looks great. +> HPL: One two day data gap that overlaps one of the ROO gaps. Data otherwise looks great. +> CSP: Data looks great with no gaps. +> +> The uncertainty of the data during this period should be low and the few missing days of data shouldn't impact / bias the results of the inversion analysis. If you want I can generate actual uncertainty numbers using Logan's interpolated method. All in all I think this is a good time period for the analysis especially given the CSP data. +> +> Ryan +> +> On Tue, Feb 12, 2019 at 4:05 PM John C. Lin wrote: +> +> Ryan, +> +> How do the data availability/quality look at HPL, ROO, and FRU between +> March~May 2016? +> +> John +> +> On 2019-02-11 10:27 a.m., Ben Fasoli wrote: +> > Nov 2015 - May 2016 were the LGR installation dates at CSP so I’d recommend that window (Mar/Apr if looking for a shorter period). We had a period later with the Nerdmobile parked next to the CSP shed but that data hasn’t been QAQC’ed yet. +> > +> > Welcome back Lewis - happy to see some eyes on the Uinta Basin dataset! +> > +> > Ben +> > +> > +> >> On Feb 11, 2019, at 10:20 AM, John C. Lin wrote: +> >> +> >> Ryan, Ben: +> >> +> >> Thinking about the time period for Lewis to carry out the CH4 inversion for the Uintah Basin, I was reminded of the importance of the Castlepeak (CSP) site. +> >> +> >> What would be a spring or fall period in which we have the most amount of quality Uintah Basin CH4 data? +> >> Would it be spring 2016? +> >> +> >> At most a month to two months would be enough for this analysis, I think. +> >> +> >> John +> >> diff --git a/Trend_foot_dCH4.r b/Trend_foot_dCH4.r new file mode 100644 index 0000000..cf53992 --- /dev/null +++ b/Trend_foot_dCH4.r @@ -0,0 +1,200 @@ +# Look at trends in footprint strength to make sure trends in Fch4 is not due to trend in footprint, as well as trend in dCH4 +# Takes output from 'Fch4_simple.r' +# V2(210505): add SITE, mettype, filterUV.TF +# 3/16/2021 by John C. Lin (John.Lin@utah.edu) + +require(ggplot2) + +################### +SITE<-"HPL" #HPL is site to focus on for calculating long-term F_CH4 +#SITE<-"ROO" +#SITE<-"CSP" +MONsub <- 4:9 # subset of months to examine (filter out wintertime) +HRs<-20:23 #[UTC] only analyze afternoon hours, following Foster papers + +mettype <- "HRRR" +#mettype <- "NAM12" +#mettype <- "WRF27" +filterUV.TF <- TRUE # filter times based on U/V (filter out times when HRRR is off--i.e., large transport errors) ? +################### + + +if(filterUV.TF&mettype!="HRRR")stop(paste0("Only extracted HRRR wind vectors; filterUV.TF needs to be set to FALSE for mettype=",mettype)) + +if(SITE=="HPL"){LAT <- 40.1434; LON <- -109.4680; zagl <- 4.06; WINDsite <- "UBHSP"; WINDsite.2015 <- "A1633"} # UBHSP site not available in 2015, so use A1633 (RedWash) site instead in 2015 +if(SITE=="ROO"){LAT <- 40.2941; LON <- -110.0090; zagl <- 4.06; WINDsite <- "QRS"} +if(SITE=="CSP"){LAT <- 40.0509; LON <- -110.0194; zagl <- 4.06; WINDsite <- "UBCSP"} + +datfilenm <- paste("Fch4_",SITE,"_daily_",mettype,".rds",sep="") # output from 'Fch4_simple.r' +dat.all <- readRDS(datfilenm) +#IV. Filter times based on U/V (filter out times when HRRR is off--i.e., large transport errors) ? +dat.all <- readRDS(paste0("Fch4_",SITE,"_",mettype,".rds")) +if(filterUV.TF){ + dat <- dat.all + Usim <- dat[,paste0("Usim.",WINDsite)]; Vsim <- dat[,paste0("Vsim.",WINDsite)] + Uobs <- dat[,paste0("Uobs.",WINDsite)]; Vobs <- dat[,paste0("Vobs.",WINDsite)] + if(SITE=="HPL"){ + # UBHSP site not available in 2015, so use A1633 (RedWash) site instead in 2015 + sel <- substring(dat$YYYYMM,1,4)=="2015" + Usim[sel] <- dat[sel,paste0("Usim.",WINDsite.2015)] + Vsim[sel] <- dat[sel,paste0("Vsim.",WINDsite.2015)] + Uobs[sel] <- dat[sel,paste0("Uobs.",WINDsite.2015)] + Vobs[sel] <- dat[sel,paste0("Vobs.",WINDsite.2015)] + } # if(SITE=="HPL"){ + isNA <- is.na(Usim)|is.na(Vsim)|is.na(Uobs)|is.na(Vobs) + # a) U/V filter: both simulated U & V have to have the same sign as the observed (i.e., wind vector has to be same QUADRANT as observed) + # SEL <- sign(Usim)!=sign(Uobs) + # SEL <- sel|(sign(Vsim)!=sign(Vobs)) + + # b) U/V filter: filter out days based on +/-45o of observed (better than quadrant method above, since not throw way data close to quadrant boundaries + WSPD.sim <- sqrt(Usim^2 + Vsim^2); WSPD.obs <- sqrt(Uobs^2 + Vobs^2) + print("Observed windspeeds (quantile):") + print(signif(quantile(WSPD.obs,c(0,0.05,0.1,0.25,0.5,0.75,1.0),na.rm=T),3)) + # see: http://tornado.sfsu.edu/geosciences/classes/m430/Wind/WindDirection.html + WDIR.sim <- (180/pi)*atan2(y=-1*Vsim,x=-1*Usim) # wind direction (-180o to +180o); direction wind is blowing FROM + WDIR.obs <- (180/pi)*atan2(y=-1*Vobs,x=-1*Uobs) # wind direction (-180o to +180o); direction wind is blowing FROM + dWDIR <- WDIR.sim - WDIR.obs + sel <- abs(dWDIR) > 180 + dWDIR[sel&!isNA] <- sign(dWDIR[sel&!isNA])*(360 - abs(dWDIR[sel&!isNA])) + SEL <- abs(dWDIR) > 45 # simulated wind direction has to be within 45o of observed + SEL[WSPD.obs < 1.0] <- FALSE # not remove data if windspeed is too low (wind anemometer not reliable) + + print(paste("Out of",sum(!isNA),"non-NA U/V times,",sum(SEL[!isNA]),"failed to pass U/V filter: ",signif(100*sum(SEL[!isNA])/sum(!isNA),3),"%")) + # assign NAs to footsum and foot.basin when failed to pass U/V filter + print(paste("Before filtering, footsum has",sum(!is.na(dat$footsum)),"non-NA values")) + dat$footsum[SEL] <- NA + print(paste(" After filtering, footsum has",sum(!is.na(dat$footsum)),"non-NA values")) + dat$foot.basin[SEL] <- NA + dat.all <- dat +} # if(filterUV.TF){ + + +dat <- subset(dat.all, as.numeric(substring(dat.all$YYYYMMDD,5,6))%in%MONsub) + +#dev.new() +#theme_set(theme_bw()) +#g <- ggplot(dat, aes(x=Time,y=foot.basin)) + geom_point(,size=1.5) + geom_smooth(method = "lm") +#g <- g + theme(plot.title=element_text(size=25), axis.title.y=element_text(size=20), +# axis.title.x=element_text(size=20), axis.text.x=element_text(size=18,angle = 30, vjust=.5), axis.text.y=element_text(size=18)) # Y axis text +#plot(g) + +# Custom function to create regression line from fit (https://sejohnston.com/2012/08/09/a-quick-and-easy-function-to-plot-lm-results-in-r/) +ggplotRegression <- function (fit,title='',caption='',ylims=NULL,xlims=NULL) { + require(ggplot2) + dev.new() + theme_set(theme_bw()) + + xslope <- signif(fit$coef[[2]], 3) + xp <- signif(summary(fit)$coef[2,4], 3) + adjR2 <- signif(summary(fit)$adj.r.squared, 3) + N <- nobs(fit) + #xsubtitle <- paste0("slope=",signif(fit$coef[[2]], 3),"; p=",signif(summary(fit)$coef[2,4], 3),"; adjR2=",signif(summary(fit)$adj.r.squared, 3)) + xsubtitle <- substitute(paste("slope=",xslope,"; p=",xp,"; adj-",R^2,"=",adjR2,"; N=",N),list(xslope=xslope,xp=xp,adjR2=adjR2,N=N)) + + g <- ggplot(fit$model, aes_string(x = names(fit$model)[2], y = names(fit$model)[1])) + + geom_point() + + stat_smooth(method = "lm", col = "blue") + + labs(title=title) + + labs(subtitle = xsubtitle) + #labs(subtitle = paste0("slope=",signif(fit$coef[[2]], 3),"; p=",signif(summary(fit)$coef[2,4], 3), "; adjR2=",signif(summary(fit)$adj.r.squared, 3))) + g <- g + theme(plot.title=element_text(size=18,hjust=0.5), plot.subtitle=element_text(size=17,hjust=0.5,color="blue"), + plot.caption=element_text(size=12), + axis.title.y=element_text(size=18), axis.text.y=element_text(size=18), + axis.title.x=element_text(size=18), axis.text.x=element_text(size=18,angle = 0, vjust=.5)) + g <- g + ylim(ylims) + xlim(xlims) + plot(g) + return(g) +} # ggplotRegression <- function (fit) { + + + +#################################################### +# I. Trends in Basin-summed footprint # +#################################################### + +xtitle <- paste0(SITE,': Trend in footprint summed over Uinta Basin\nmons: ',paste(MONsub,collapse=","), + '; UThrs: ',paste(HRs,collapse=","), '\n mettype=',mettype) +xcaption <- paste0('filterUV.TF=',filterUV.TF) + +######################### +# a) Daily time series +#xfit <- lm(foot.basin ~ Time, data=dat) +#ggplotRegression(xfit) + + + +######################### +# b) Monthly time series +foot.basin.mon <- tapply(dat$foot.basin, dat$YYYYMM, mean, na.rm=T) +Time <- as.POSIXct(paste0(names(foot.basin.mon),"01"), format("%Y%m%d"),tz="GMT") +xfit <- lm(foot.basin.mon ~ Time) +xfit$coef[[2]] <- xfit$coef[[2]]*(365*24*60*60) # slope is in [1/s]; convert to [1/year] +ylims <- c(0.3,0.9) +xlims <- as.POSIXct(strptime(x=paste0(c(2015,2021),"-01-01"), format="%Y-%m-%d"),tz="GMT") +xbreaks <- as.POSIXct(strptime(x=paste0(2015:2021,"-01-01"), format="%Y-%m-%d"),tz="GMT") +g <- ggplotRegression(xfit,title=xtitle,caption=xcaption,ylim=ylims,xlim=xlims) +g <- g + labs(x="Year",y=expression(paste("Total Footprint in Basin [",mu,"mole m"^-2," s"^-1,"]"))) +g <- g + scale_x_continuous(name="Year", breaks=xbreaks, labels=2015:2021) +plot(g) +ggsave(paste0('Trend_foot_month_',mettype,'.png')) + + + +######################### +# c) Annual time series +foot.basin.yr <- tapply(dat$foot.basin, substring(dat$YYYYMM,1,4), mean, na.rm=T) +Year <- as.numeric(names(foot.basin.yr)) +ylims <- c(0.3,0.9) +xlims <- c(2015,2021) +g <- ggplotRegression(lm(foot.basin.yr ~ Year),title=xtitle,caption=xcaption,ylim=ylims,xlim=xlims) +g <- g + labs(x="Year",y=expression(paste("Total Footprint in Basin [",mu,"mole m"^-2," s"^-1,"]"))) +g <- g + scale_x_continuous(name="Year", breaks=2015:2021, labels=2015:2021) +plot(g) +ggsave(paste0('Trend_foot_annual_',mettype,'.png')) + + + +#################################################### +# II. Trends in Methane Enhancement (dCH4) # +#################################################### + +# for published version, turn off the additional information and caption +#xtitle <- paste0('SITE =',SITE,'; Trend in CH4 enhancement over FRU\nmons: ',paste(MONsub,collapse=","),'; UThrs: ',paste(HRs,collapse=",")) +#xcaption <- paste0('filterUV.TF=',filterUV.TF) +xtitle <- substitute(paste(SITE,": Trend in CH"[4]," enhancement over FRU"),list(SITE=SITE)) +xcaption <- paste0('') +######################### +# b) Monthly time series +dCH4 <- 1000*dat$dCH4 # CH4 enhancement [ppb] +dCH4.mon <- tapply(dCH4, dat$YYYYMM, mean, na.rm=T) +Time <- as.POSIXct(paste0(names(dCH4.mon),"01"), format("%Y%m%d"),tz="GMT") +xfit <- lm(dCH4.mon ~ Time) +xfit$coef[[2]] <- xfit$coef[[2]]*(365*24*60*60) # slope is in [1/s]; convert to [1/year] +ylims <- c(0,100) +xlims <- as.POSIXct(strptime(x=paste0(c(2015,2021),"-01-01"), format="%Y-%m-%d"),tz="GMT") +xbreaks <- as.POSIXct(strptime(x=paste0(2015:2021,"-01-01"), format="%Y-%m-%d"),tz="GMT") +g <- ggplotRegression(xfit,title=xtitle,caption=xcaption,ylim=ylims,xlim=xlims) +g <- g + labs(x="Year",y=expression(paste(Delta,"CH"[4]," [ppb]"))) +g <- g + scale_x_continuous(name="Year", breaks=xbreaks, labels=2015:2021) +#g <- g + xlim(xlims) +plot(g) +ggsave(paste0('Trend_dCH4_month_',mettype,'.png')) + + + +######################### +# c) Annual time series +dCH4 <- 1000*dat$dCH4 # CH4 enhancement [ppb] +dCH4.yr <- tapply(dCH4, substring(dat$YYYYMM,1,4), mean, na.rm=T) +Year <- as.numeric(names(dCH4.yr)) +xfit <- lm(dCH4.yr ~ Year) +ylims <- c(0,100) +xlims <- c(2015,2022) +g <- ggplotRegression(xfit,title=xtitle,caption=xcaption,ylim=ylims,xlim=xlims) +g <- g + labs(x="Year",y=expression(paste(Delta,"CH4 [ppb]"))) +g <- g + scale_x_continuous(name="Year", breaks=2015:2021, labels=2015:2021) +plot(g) +ggsave(paste0('Trend_dCH4_annual_',mettype,'.png')) + + + diff --git a/average_foot_monthly.r b/average_foot_monthly.r new file mode 100644 index 0000000..027a4b9 --- /dev/null +++ b/average_foot_monthly.r @@ -0,0 +1,140 @@ +# Generate MONTHLY-averaged footprint MAPS +# V2(210428): introduce 'mettype' to distinguish between met files driving STILT (to assess transport error)--e.g., "HRRR", "NARR", "NAM12", "WRF27" +# April 18th, 2021 by John C. Lin (John.Lin@utah.edu) + + +################ +YEARs <- 2015:2020 +#SITE<-"HPL" +#SITE<-"ROO" +SITE<-"CSP" +if(SITE=="CSP")YEARs <- c(2016,2019,2020) +if(SITE=="ROO")YEARs <- 2015:2019 + +# meteorology driving STILT +mettype <- "HRRR" +#mettype <- "NAM12" +#mettype <- "WRF27" +if(mettype=="HRRR")STILTdir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_HYSPLIT-STILT_HRRR/out/" +if(mettype=="NARR")STILTdir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_HYSPLIT-STILT_NARR/out/" +if(mettype=="NAM12")STILTdir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_HYSPLIT-STILT_NAM12/out/" +if(mettype=="WRF27")STILTdir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_HYSPLIT-STILT_WRF27km/out/" + +obsdir<-"/uufs/chpc.utah.edu/common/home/lin-group4/jcl/SimCity/" #where CH4 observations reside (to filter out times when SITE measured NAs) + +# subset of hours [UTC]; if NULL, then all hours +# UThrs <- NULL +UThrs <- c(20,21,22,23) #[UTC] only analyze afternoon hours, following Foster papers + +finalresultname <- paste0(SITE,'_foot_monave_',mettype,'.rds') +################ + +outputdir<-paste0(STILTdir,"/",SITE,"/by-id/") + +require(ncdf4);require(ggplot2) +# Footprint grid settings from 'run_stilt_hysplit-stiltV1.r' +xmn <- -112.5; xmx <- -108 +ymn <- 38.5; ymx <- 41.5 +xres <- 0.01; yres <- xres +FLAT <- seq(ymn+yres/2,ymx-yres/2,yres) +FLON <- seq(xmn+xres/2,xmx-xres/2,xres) + +XFILES <- list.files(outputdir) + +# Months and Years to grid +MMs <- formatC(1:12,flag="0",width=2) +YYYYMMs <- paste0(rep(YEARs,each=length(MMs)),rep(MMs,length(YEARs))) +RESULT.MON <- array(data=0, dim=c(length(FLON),length(FLAT),length(YYYYMMs))) +dimnames(RESULT.MON) <- list(FLON,FLAT,YYYYMMs) + +for(yy in 1:length(YEARs)){ + YEAR <- YEARs[yy] + # Simulation timing, yyyy-mm-dd HH:MM:SS (UTC) + t_start <- paste0(YEAR,'-01-01 00:00:00') + t_end <- paste0(YEAR,'-12-31 23:00:00') + run_times <- seq(from = as.POSIXct(t_start, tz = 'UTC'), to = as.POSIXct(t_end, tz = 'UTC'), by = 'hour') + if(!is.null(UThrs)){ + sel <- substring(run_times,12,13)%in%UThrs #only analyze afternoon hours, following Foster papers + run_times <- run_times[sel] + } # if(!is.null(UThrs)){ + + # read in observations + objname <- paste0("SimCity_CH4_allsites_hrly_",YEARs[yy]) + print(paste("Reading in.....",objname)) + obs <- getr(objname,path=obsdir)[,c("Time",paste0("CH4_",SITE))] + + # skip footprint if CH4_SITE is NA + run_times.sub <- subset(run_times,run_times%in%obs$Time) + Times2rm <- obs[is.na(obs[,paste0("CH4_",SITE)]),'Time'] + sel <- run_times.sub%in%Times2rm + print(paste('Skipping',sum(sel),'hrs out of total of',length(run_times.sub),'hrs')) + run_times.sub <- run_times.sub[!sel] + + #!!! skip footprint if transport error is too large !!! + +for(mm in 1:12){ + YYYYMM <- paste0(YEAR,formatC(mm,flag="0",width=2)) + run_times.tt <- run_times.sub[format(run_times.sub,"%Y%m")==YYYYMM] + foot.sum <- 0 + N <- 0 #start counter +for(tt in 1:length(run_times.tt)){ + YYYYMMDDHH <- format(run_times.tt[tt],"%Y%m%d%H") + sel <- substring(XFILES,1,10)==YYYYMMDDHH + if(sum(is.na(XFILES[sel])>0)){print(paste("XFILES[sel] is NA for:",YYYYMMDDHH));next} + ncfilenm <- paste0(outputdir,"/",XFILES[sel],"/",XFILES[sel],"_foot.nc") + if(length(ncfilenm)>1){print(paste("more than 1 ncfilenm:",paste(ncfilenm,collapse="; ")))} + if(!file.exists(ncfilenm)){print(paste("no footprint output for:",YYYYMMDDHH));next} + # retrieve footprint + footfile <- nc_open(ncfilenm) + # print(paste(N,"Reading in:",ncfilenm)) + foot.all <- ncvar_get(footfile,"foot") + foot <- apply(foot.all,c(1,2),sum) #sum over back times + foot.sum <- foot.sum + foot + N <- N+1 + print(paste(N,"Processed",YYYYMMDDHH)) + nc_close(footfile) +} # for(tt in 1:length(run_times.tt)){ + RESULT.MON[,,YYYYMM] <- foot.sum/N + saveRDS(RESULT.MON,finalresultname) + print(paste(finalresultname,"saved")) + gc() +} # for(mm in 1:12){ + gc() +} # for(yy in 1:length(YEARs)){ + + +if(FALSE){ + require(fields) + + SITE <- "CSP" + if(SITE=="HPL"){LAT <- 40.1434; LON <- -109.4680; zagl <- 4.06} + if(SITE=="ROO"){LAT <- 40.2941; LON <- -110.0090; zagl <- 4.06} + if(SITE=="CSP"){LAT <- 40.0509; LON <- -110.0194; zagl <- 4.06} + + #YYYYMM <- "201604" + YYYYMM <- "202004" + resultnm <- paste0(SITE,'_foot_monave.rds') + dat.all <- readRDS(resultnm) + dat <- dat.all[,,YYYYMM] + dat[dat==0] <- NA + + #a) log-10 plot + dat.log10 <- log10(dat) + zlims <- c(-7,-1) + #zlims <- NULL + flon <- as.numeric(rownames(dat)) + flat <- as.numeric(colnames(dat)) + if(!is.null(zlims)){dat.log10[dat.log10zlims[2]]<-zlims[2]} + dev.new();par(cex.lab=1.3,cex.axis=1.3,cex.main=1.3) + image.plot(flon,flat,dat.log10,main=paste("Site=",SITE,"\n",YYYYMM),zlim=zlims) + points(LON,LAT,pch=17,col="black") + + #b) linear plot + #zlims <- NULL + zlims <- c(0.00001,0.0001) + if(!is.null(zlims)){dat[datzlims[2]]<-zlims[2]} + dev.new();par(cex.lab=1.3,cex.axis=1.3,cex.main=1.3) + image.plot(flon,flat,dat,main=paste("Site=",SITE,"\n",YYYYMM),zlim=zlims) + points(LON,LAT,pch=17,col="black") + +} # if(FALSE){ diff --git a/compare_foots_HYSPLIT-STILT_vs_STILT.r b/compare_foots_HYSPLIT-STILT_vs_STILT.r new file mode 100644 index 0000000..b89915d --- /dev/null +++ b/compare_foots_HYSPLIT-STILT_vs_STILT.r @@ -0,0 +1,150 @@ +#Compare footprints from newly merged HYSPLIT-STILT executable versus STILT executable +#Based on "plot_foot_CH4_inversion_UintahV1.r" +#September 6th, 2020 by John C. Lin (John.Lin@utah.edu) + +######################### +#SITE<-"INDY" +#SITE<-"HPL" +#SITE<-"ROO" +SITE<-"CSP" +#SITE<-"SLC" + +# HYSPLIT-STILT output directory +outputdir<-paste0("CH4_inversion_Uintah_HYSPLIT-STILT/out/",SITE,"/") +# STILT output directory +outputdir.s<-paste0("CH4_inversion_Uintah_STILT/out/",SITE,"/") + +#outputdir<-paste0("out/",SITE,"/by-id/") +pngfileTF<-FALSE; overwritepngfileTF<-FALSE +GoogleMapTF<-FALSE + +#lat/lon range of plot +XLIMS<-c(-90,-80);YLIMS<-c(38,42) +#range of footprint plot (log10 limits) +ZLIMS<-c(-7,-1) + +# Simulation timing, yyyy-mm-dd HH:MM:SS (UTC) +# t_start <- '2015-07-01 00:00:00' +# t_end <- '2020-08-31 23:00:00' +t_start <- '2016-01-01 00:00:00' +t_end <- '2016-12-31 23:00:00' +run_times <- seq(from = as.POSIXct(t_start, tz = 'UTC'), + to = as.POSIXct(t_end, tz = 'UTC'), + by = 'hour') +sel<-substring(run_times,12,13)%in%c("20","21","22","23") #only analyze afternoon hours, following Foster papers +run_times<-run_times[sel] + +calc.footTOT.TF <- FALSE +######################### + + +require(fields);require(ncdf4) +require(maps);require(RgoogleMaps) + +#if(!SITE%in%c("WBB","BOS","INDY","JES","SLC","PORT","LA","SF"))stop(paste("Not valid site:",SITE)) +if(SITE=="HPL"){lati <- 40.1434; long <- -109.4680; zagl <- 4.06} +if(SITE=="ROO"){lati <- 40.2941; long <- -110.0090; zagl <- 4.06} +if(SITE=="CSP"){lati <- 40.0509; long <- -110.0194; zagl <- 4.06} +if(SITE=="WBB"){lati <- 40.766 ; long <- -111.849; zagl <- 35} +if(SITE=="BOS"){lati <- 42.3470; long <- -71.0840; zagl <- 20} +if(SITE=="INDY"){lati <- 39.7833; long <- -86.1652; zagl <- 20} +if(SITE=="JES"){lati <- 39.1723; long <- -76.7765; zagl <- 20} +if(SITE=="SLC"){lati <- 40.6539; long <- -111.8878; zagl <- 20} +if(SITE=="LA"){lati <- 33.8700; long <- -118.2800; zagl <- 20} +if(SITE=="PORT"){lati <- 45.5100; long <- -122.6000; zagl <- 20} +if(SITE=="SF"){lati <- 37.7956; long <- -122.2794; zagl <- 20} + +xmn <- long-2; xmx <- long+2 +ymn <- lati-2; ymx <- lati+2 +xres <- 0.01 +yres <- xres + +if(calc.footTOT.TF){ +#loop over all receptor times +YYYYMMDDHHMMs <- format(run_times,format="%Y%m%d%H%M") +RESULT <- data.frame(Time=run_times,HYSPLIT.STILT=NA,STILT=NA) +for(i in 1:length(run_times)){ + YYYYMMDDHHMM<-YYYYMMDDHHMMs[i] + print(paste(" Date:",YYYYMMDDHHMM)) + ############################### + #I. HYSPLIT-STILT results + receptor <- paste0(YYYYMMDDHHMM,"_",long,"_",lati,"_",zagl) + ncfilenm<-paste0(outputdir,"/by-id/",receptor,"/",receptor,"_foot.nc") # e.g., 201802042100_-109.468_40.1434_4.06_foot.nc + if(file.exists(ncfilenm)){ + #read in footprint and calculate total + footfile<-nc_open(ncfilenm) + flat<-ncvar_get(footfile,"lat");flon<-ncvar_get(footfile,"lon") + foot.all<-ncvar_get(footfile,"foot") + footsum<-apply(foot.all,c(1,2),sum) #sum over backward time, but preserve spatial structure + footTOT<-sum(footsum) #TOTAL footprint + # foot.log<-log10(footsum);foot.log[footsum==0]<-NA + RESULT[i,"HYSPLIT.STILT"] <- footTOT + } else { + print(paste("no HYSPLIT-STILT footprint output for:",receptor)) + } # if(file.exists(ncfilenm)){ + + ############################### + #II. STILT results + receptor <- paste0(substring(YYYYMMDDHHMM,1,10),"_",long,"_",lati,"_",zagl) + ncfilenm<-paste0(outputdir.s,"/by-id/",receptor,"/",receptor,"_foot.nc") # e.g., 2018020421_-109.468_40.1434_4.06_foot.nc + if(file.exists(ncfilenm)){ + #read in footprint and calculate total + footfile<-nc_open(ncfilenm) + flat<-ncvar_get(footfile,"lat");flon<-ncvar_get(footfile,"lon") + foot.all<-ncvar_get(footfile,"foot") + footsum<-apply(foot.all,c(1,2),sum) #sum over backward time, but preserve spatial structure + footTOT<-sum(footsum) #TOTAL footprint + # foot.log<-log10(footsum);foot.log[footsum==0]<-NA + RESULT[i,"STILT"] <- footTOT + } else { + print(paste("no STILT footprint output for:",receptor)) + } # if(file.exists(ncfilenm)){ + + saveRDS(RESULT,"footTOT_HYSPLIT-STILT_STILT.RDS") + gc() +} # for(i in 1:length(run_times)){ +} # if(calc.footTOT.TF){ + +dat <- readRDS("footTOT_HYSPLIT-STILT_STILT.RDS") +# plot time series +dev.new(); par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3) +plot(dat[,c("Time","HYSPLIT.STILT")],pch=16,cex=0.5,xlab="Time",ylab="Total Footprint [ppm/(umole/m2/s)]",main=paste("Site:",SITE)) +points(dat[,c("Time","STILT")],pch=16,cex=0.5,col="darkgray") +legend("topright",c("HYSPLIT-STILT","STILT"),text.col=c("black","darkgray"),col=c("black","darkgray"),pch=16) + +# plot differences +dfootTOT <- dat$HYSPLIT.STILT-dat$STILT +dev.new(); par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3) +plot(x=dat$Time,y=dfootTOT,pch=16,cex=0.5,xlab="Time",ylab="Total Footprint Difference (HYSPLIT.STILT minus STILT)", + main=paste("Site:",SITE)) + +# plot differences--subset of months +MONs <- c(4:10) +sel <- as.numeric(format(dat$Time,format="%m"))%in%MONs +dev.new(); par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3) +xmain <- paste("Site:",SITE,"\nMonths:",paste(MONs,collapse=" ")) +plot(x=dat$Time[sel],y=dfootTOT[sel],pch=16,cex=0.5,xlab="Time", + ylab="Total Footprint Diff (HYSPLIT.STILT minus STILT)",main=xmain) +# plot distribution of differences +dev.new(); par(cex.axis=1.3,cex.lab=1.2,cex.main=1.3) +hist(dfootTOT[sel],xlab="Total Footprint Diff (HYSPLIT.STILT minus STILT)\n[ppm/(umole/m2/s)]",main=xmain) +xbias <- mean(dfootTOT[sel],na.rm=T) +xstdev <- stdev(dfootTOT[sel],na.rm=T) +abline(v=xbias,col="red") +legend("topleft",c(paste("mean(HYSPLIT.STILT - STILT):",signif(xbias,3)), + paste("stdev(HYSPLIT.STILT - STILT):",signif(xstdev,3))),cex=1.1,bty="n") +# plot distribution of PERCENTAGE differences +dfootTOT.perc <- 100*(dat$HYSPLIT.STILT-dat$STILT)/dat$HYSPLIT.STILT +dev.new(); par(cex.axis=1.3,cex.lab=1.2,cex.main=1.3) +hist(dfootTOT.perc[sel],xlab="% Footprint Diff: 100*(HYSPLIT.STILT - STILT)/HYSPLIT.STILT",main=xmain) +xbias <- mean(dfootTOT.perc[sel],na.rm=T) +xstdev <- stdev(dfootTOT.perc[sel],na.rm=T) +abline(v=xbias,col="red") +legend("topleft",c(paste("mean:",signif(xbias,3),"%"), + paste("stdev:",signif(xstdev,3),"%")),cex=1.2,bty="n") + + + + + + diff --git a/compare_trajs_HYSPLIT-STILT_vs_STILT.r b/compare_trajs_HYSPLIT-STILT_vs_STILT.r new file mode 100644 index 0000000..6e753f9 --- /dev/null +++ b/compare_trajs_HYSPLIT-STILT_vs_STILT.r @@ -0,0 +1,168 @@ +#Compare trajectories from newly merged HYSPLIT-STILT executable versus STILT executable +#Based on "plot_foot_CH4_inversion_UintahV1.r" +#September 6th, 2020 by John C. Lin (John.Lin@utah.edu) + +######################### +#SITE<-"INDY" +SITE<-"HPL" +#SITE<-"ROO" +#SITE<-"CSP" +#SITE<-"SLC" + +# HYSPLIT-STILT output directory +outputdir<-paste0("CH4_inversion_Uintah_HYSPLIT-STILT/out/",SITE,"/by-id/") +# STILT output directory +outputdir.s<-paste0("CH4_inversion_Uintah_STILT/out/",SITE,"/by-id/") + +#outputdir<-paste0("out/",SITE,"/by-id/") +pngfileTF<-FALSE; overwritepngfileTF<-FALSE +GoogleMapTF<-FALSE + +#lat/lon range of plot +#XLIMS<-c(-115,-109);YLIMS<-c(37,42.5) +#XLIMS<-c(-113.5,-111.2);YLIMS<-c(40,41.7) +#XLIMS<-c(-112.3,-111.5);YLIMS<-c(40,41.3) +#XLIMS<-c(-112.7,-111.5);YLIMS<-c(40,41.3) +XLIMS<-c(-90,-80);YLIMS<-c(38,42) +#range of footprint plot (log10 limits) +ZLIMS<-c(-7,-1) +######################### + +require(fields);require(ncdf4) +require(maps);require(RgoogleMaps) + +#if(!SITE%in%c("WBB","BOS","INDY","JES","SLC","PORT","LA","SF"))stop(paste("Not valid site:",SITE)) +if(SITE=="HPL"){lati <- 40.1434; long <- -109.4680; zagl <- 4.06} +if(SITE=="ROO"){lati <- 40.2941; long <- -110.0090; zagl <- 4.06} +if(SITE=="CSP"){lati <- 40.0509; long <- -110.0194; zagl <- 4.06} +if(SITE=="WBB"){lati <- 40.766 ; long <- -111.849; zagl <- 35} +if(SITE=="BOS"){lati <- 42.3470; long <- -71.0840; zagl <- 20} +if(SITE=="INDY"){lati <- 39.7833; long <- -86.1652; zagl <- 20} +if(SITE=="JES"){lati <- 39.1723; long <- -76.7765; zagl <- 20} +if(SITE=="SLC"){lati <- 40.6539; long <- -111.8878; zagl <- 20} +if(SITE=="LA"){lati <- 33.8700; long <- -118.2800; zagl <- 20} +if(SITE=="PORT"){lati <- 45.5100; long <- -122.6000; zagl <- 20} +if(SITE=="SF"){lati <- 37.7956; long <- -122.2794; zagl <- 20} + +xmn <- long-2; xmx <- long+2 +ymn <- lati-2; ymx <- lati+2 +xres <- 0.01 +yres <- xres + +#loop over the YYYYMMDD folders +YYYYMMDDs<-list.files(outputdir,pattern="2018") #e.g., "201801012000_-109.468_40.1434_4.06" +YYYYMMDDs.s<-list.files(outputdir.s,pattern="2018") #e.g., "2018010120_-109.468_40.1434_4.06" NOTE: does NOT have MINUTE in receptor time +YYYYMMDDs<-YYYYMMDDs[substring(YYYYMMDDs,1,10)%in%substring(YYYYMMDDs.s,1,10)] #make sure that both STILT & HYSPLIT-STILT results are available +# YYYYMMDDs<-sample(YYYYMMDDs,3) #randomly sample a few dates +YYYYMMDDs<-YYYYMMDDs[substring(YYYYMMDDs,1,10)%in%c("2018063020","2018010123","2018053121")] +for(i in 1:length(YYYYMMDDs)){ +#for(i in 875:length(YYYYMMDDs)){ +#for(i in c(1,50)){ +#for(i in 13:length(YYYYMMDDs)){ + ncfilenm<-paste0(outputdir,"/",YYYYMMDDs[i],"/",YYYYMMDDs[i],"_foot.nc") + if(!file.exists(ncfilenm)){print(paste("no footprint output for:",YYYYMMDDs[i]));next} + + #Generate footprint + footfile<-nc_open(ncfilenm) + lat<-ncvar_get(footfile,"lat");lon<-ncvar_get(footfile,"lon") + foot.all<-ncvar_get(footfile,"foot") + flat<-lat;flon<-lon + footsum<-apply(foot.all,c(1,2),sum) #sum over backward time + foot.log<-log10(footsum);foot.log[footsum==0]<-NA + + #Read in trajectories + trajfile<-paste0(outputdir,"/",YYYYMMDDs[i],"/",YYYYMMDDs[i],"_traj.rds") + if(!file.exists(trajfile)){print(paste("cannot find HYSPLIT-STILT trajectory file",trajfile));next} + # change from YYYYMMDDHHmm to YYYYMMDDHH for STILT receptor + trajname.s <- paste0(substring(YYYYMMDDs[i],1,10),substring(YYYYMMDDs[i],13,nchar(YYYYMMDDs[i]))) + trajfile.s<-paste0(outputdir.s,"/",trajname.s,"/",trajname.s,"_traj.rds") + if(!file.exists(trajfile.s)){print(paste("cannot find STILT trajectory file",trajfile.s));next} + # HYSPLIT-STILT trajectories + dat.all<-readRDS(file=trajfile) + receptor<-dat.all$receptor + lon0<-receptor$long;lat0<-receptor$lati;agl0<-receptor$zagl + receptor$run_time<-strftime(receptor$run_time,"%Y-%m-%d %H:%M",tz="UTC") + dat<-dat.all$particle + # STILT trajectories + dat.all.s<-readRDS(file=trajfile.s) + dat.s<-dat.all.s$particle + + #Check whether PNG file exists or not + if(pngfileTF){ + pngfile<-paste0(outputdir,"/",YYYYMMDDs[i],"/",YYYYMMDDs[i],".png") + if(!overwritepngfileTF){if(file.exists(pngfile)){print(paste("PNG file already exists!; skip:",pngfile));next}} + } #if(pngfileTF){ + + lats<-dat$lati;lons<-dat$long +if(!GoogleMapTF){ + dev.new();par(cex.axis=1.3,cex.lab=1.3,cex.main=1.1) + #image.plot(flon,flat,foot.log,xlab="Lon",ylab="Lat") + #image.plot(flon,flat,foot.log,xlab="Lon",ylab="Lat",col=terrain.colors(10)) + #image.plot(flon,flat,foot.log,xlab="Lon",ylab="Lat",col=snow.colors(10)) + #image.plot(flon,flat,foot.log,xlab="Lon",ylab="Lat",col=larry.colors()) + image.plot(flon,flat,foot.log,xlab="Lon",ylab="Lat",col=rev(gray(0:10/10))) + #add trajectories + # a) STILT-HYSPLIT trajectories + points(lons,lats,pch=16,col="darkblue",cex=0.4) + points(lons[dat$foot>0],lats[dat$foot>0],pch=16,col="lightblue",cex=0.2) + # b) STILT trajectories + points(dat.s$long,dat.s$lati,pch=16,col="darkgreen",cex=0.2) + points(dat.s$long[dat.s$foot>0],dat.s$lati[dat.s$foot>0],pch=16,col="lightgreen",cex=0.2) + #add receptor location + points(lon0,lat0,pch=17,col="red",cex=1.5) + map("state",add=TRUE) #;map("world",add=TRUE) + map.cities(minpop=100000,cex=1.2) + title(main=paste("Footprint (STILT model) at receptor (red triangle):\n",receptor$run_time,"UTC ", + round(lon0,digits=3),"LON ",round(lat0,digits=3),"LAT ",round(agl0,digits=3),"m-AGL"),cex.main=1.1) + # legend("topright",c("HYSPLIT-STILT","STILT"),lty=1,pch=16,lwd=2,col=c("black","darkgray")) + legend("topright",c("HYSPLIT-STILT","HYSPLIT-STILT(foot>0)","STILT","STILT(foot>0)"), + lty=1,pch=16,lwd=2,col=c("darkblue","lightblue","darkgreen","lightgreen"),text.col=c("darkblue","lightblue","darkgreen","lightgreen")) +}else{ + #map defined by bounding box + #bb <- qbbox(lat=lats, lon=lons) #use trajectories to define map + #bb <- qbbox(lat=flat, lon=flon) #use footprint to define map + bb <- qbbox(lat=YLIMS, lon=XLIMS) #use footprint to define map + MyMap <- GetMap.bbox(bb$lonR, bb$latR, maptype="terrain", frame=TRUE,GRAYSCALE=FALSE) + PlotOnStaticMap(MyMap,mar=c(2,1,3,1)) + + #Add trajectories + site.xy<- LatLon2XY.centered(MyMap, lat=as.numeric(lats), lon=as.numeric(lons)) + points(site.xy$newX,site.xy$newY,pch=16,col="darkgray",cex=0.2) + + #Add receptor location + site.xy<- LatLon2XY.centered(MyMap, lat=as.numeric(lat0), lon=as.numeric(lon0)) + points(site.xy$newX,site.xy$newY,pch=17,col="black",cex=1.5) + + title(main=paste("CO2-USA Site Footprint (STILT model) at receptor (black triangle):\n",receptor$run_time,"UTC ", + round(lon0,digits=3),"LON ",round(lat0,digits=3),"LAT ",round(agl0,digits=3),"m-AGL"),cex.main=1.1) + + #Add footprint + alpha<-0.45 #transparency + lats.m<-matrix(rep(flat,length(flon)),byrow=T,ncol=length(flat)) + lons.m<-matrix(rep(flon,length(flat)),ncol=length(flat)) + xys<- LatLon2XY.centered(MyMap, lat=as.vector(lats.m), lon=as.vector(lons.m)) #converting lat/lon coordinates to map coordinates + image.coords.x<-matrix(xys$newX,ncol=length(flat)) + image.coords.y<-matrix(xys$newY,ncol=length(flat)) + COLS<-c("white","darkgray","lightblue","RoyalBlue","darkgreen","yellow","orange","red") #colorscale + colsc<-designer.colors(n=64,alpha=alpha,col=COLS) + colsc.legend<-designer.colors(n=64,alpha=1.0,col=COLS) + + #generate image + image(x=image.coords.x[,1],y=image.coords.y[1,],z=foot.log, + col=colsc,ylim=c(MyMap$BBOX$ll[1],MyMap$BBOX$ur[1]),xlim=c(MyMap$BBOX$ll[2],MyMap$BBOX$ur[2]),add=T,zlim=ZLIMS) + + #Add legend + image.plot(x=image.coords.x[,1],y=image.coords.y[1,],z=foot.log, + col=colsc.legend,ylim=c(MyMap$BBOX$ll[1],MyMap$BBOX$ur[1]),xlim=c(MyMap$BBOX$ll[2],MyMap$BBOX$ur[2]),add=T, + legend.lab="log10(footprint)",legend.line=2,zlim=ZLIMS,legend.only=T,horizontal=T,legend.width=0.8) +} #if(!GoogleMapTF){ + + if(pngfileTF){ + dev.copy(png,filename=pngfile);dev.off() + print(paste(pngfile,"produced")) + } #if(pngfileTF){ + nc_close(footfile) + gc() +} #for(i in 1:length(YYYYMMDDs)){ + + diff --git a/convolve_foot_wellinfo.r b/convolve_foot_wellinfo.r new file mode 100644 index 0000000..f53b2ba --- /dev/null +++ b/convolve_foot_wellinfo.r @@ -0,0 +1,194 @@ +# Convolve footprints with well information to get footprint-weighted quantities +# V2(210406): generate MONTHLY foot*wellinfo MAPS, in addition to hourly, sum(foot*wellinfo) +# V3(210427): introduce 'mettype' to distinguish between met files driving STILT (to assess transport error)--e.g., "HRRR", "NARR", "NAM12", "WRF27" +# March 28th, 2021 by John C. Lin (John.Lin@utah.edu) + + +################ +YEARs <- 2015:2020 +#YEARs <- 2019:2020 +#SITE<-"HPL" +#SITE<-"ROO" +SITE<-"CSP" +if(SITE=="CSP")YEARs <- c(2016,2019,2020) +if(SITE=="ROO")YEARs <- 2015:2019 + +# meteorology driving STILT +mettype <- "HRRR" +#mettype <- "NAM12" +#mettype <- "WRF27" +if(mettype=="HRRR")STILTdir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_HYSPLIT-STILT_HRRR/out/" +if(mettype=="NARR")STILTdir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_HYSPLIT-STILT_NARR/out/" +if(mettype=="NAM12")STILTdir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_HYSPLIT-STILT_NAM12/out/" +if(mettype=="WRF27")STILTdir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_HYSPLIT-STILT_WRF27km/out/" +obsdir<-"/uufs/chpc.utah.edu/common/home/lin-group4/jcl/SimCity/" #where CH4 observations reside (to filter out times when SITE measured NAs) + +# subset of hours [UTC]; if NULL, then all hours +# UThrs <- NULL +UThrs <- c(20,21,22,23) #[UTC] only analyze afternoon hours, following Foster papers + +# output generated by 'grid_well_info_foot.r' +VARS <- c("Mcf.Gas","Bbls.Oil","Gas.Well","Oil.Well","Gas.Well_Producing","Oil.Well_Producing") +hourlyTF <- TRUE +################ + +outputdir<-paste0(STILTdir,"/",SITE,"/by-id/") + +require(ncdf4);require(ggplot2) +# Footprint grid settings from 'run_stilt_hysplit-stiltV1.r' +xmn <- -112.5; xmx <- -108 +ymn <- 38.5; ymx <- 41.5 +xres <- 0.01; yres <- xres +FLAT <- seq(ymn+yres/2,ymx-yres/2,yres) +FLON <- seq(xmn+xres/2,xmx-xres/2,xres) + +XFILES <- list.files(outputdir) + +############################################################################################################### +if(hourlyTF){ +for(yy in 1:length(YEARs)){ + YEAR <- YEARs[yy] + # Simulation timing, yyyy-mm-dd HH:MM:SS (UTC) + t_start <- paste0(YEAR,'-01-01 00:00:00') + t_end <- paste0(YEAR,'-12-31 23:00:00') + finalresultname <- paste0(SITE,'_foot.wellinfo_',YEAR,'_',mettype,'.rds') + run_times <- seq(from = as.POSIXct(t_start, tz = 'UTC'), to = as.POSIXct(t_end, tz = 'UTC'), by = 'hour') +if(!is.null(UThrs)){ + sel <- substring(run_times,12,13)%in%UThrs #only analyze afternoon hours, following Foster papers + run_times <- run_times[sel] +} # if(!is.null(UThrs)){ + + # object to store final results + RESULT <- matrix(NA,ncol=length(VARS),nrow=length(run_times),dimnames=list(NULL,VARS)) + RESULT <- data.frame(Time=run_times,RESULT) +for(vv in 1:length(VARS)){ +#for(vv in 4:length(VARS)){ + var <- VARS[vv] + print(paste("----------------- Processing: var=",var," (HOURLY) --------------------")) + well.all <- readRDS(paste0(var,"_gridded.rds")) +for(i in 1:nrow(RESULT)){ + print(paste("time: ",RESULT[i,"Time"])) + YYYYMMDDHH <- format(RESULT[i,"Time"],"%Y%m%d%H") + sel <- substring(XFILES,1,10)==YYYYMMDDHH + ncfilenm <- paste0(outputdir,"/",XFILES[sel],"/",XFILES[sel],"_foot.nc") + if(!file.exists(ncfilenm)){print(paste("no footprint output for:",YYYYMMDDHH));next} + + # retrieve footprint + footfile <- nc_open(ncfilenm) + foot.all <- ncvar_get(footfile,"foot") + foot <- apply(foot.all,c(1,2),sum) #sum over back times + # pull out correct month from well info + YYYYMM <- substring(YYYYMMDDHH,1,6) + well <- well.all[,,YYYYMM] + if((nrow(well)!=nrow(foot))|(ncol(well)!=ncol(foot))){print("dimensions of foot not match those of wellinfo; skip!");next} + # convolve footprint with well info + footvar <- sum(foot*well) + RESULT[i,var] <- footvar +} # for(i in 1:nrow(RESULT)){ + saveRDS(RESULT,finalresultname) + print(paste(finalresultname,"generated")) + gc() +} # for(vv in 1:length(VARS)){ + gc() +} #for(yy in 1:length(YEARs)){ +} # if(hourlyTF){ + + +############################################################################################################### +# Calculating MONTHLY foot*wellinfo MAPS, in addition to hourly, sum(foot*wellinfo) +for(vv in 1:length(VARS)){ + var <- VARS[vv] + print(paste("----------------- Processing: var=",var," (MONTHLY) --------------------")) + well.all <- readRDS(paste0(var,"_gridded.rds")) + # Months and Years to grid + MMs <- formatC(1:12,flag="0",width=2) + YYYYMMs <- paste0(rep(YEARs,each=length(MMs)),rep(MMs,length(YEARs))) + RESULT.MON <- array(data=0, dim=c(length(FLON),length(FLAT),length(YYYYMMs))) + dimnames(RESULT.MON) <- list(FLON,FLAT,YYYYMMs) + + finalresultname.gridded <- paste0(SITE,'_foot_',var,'_gridded_',mettype,'.rds') +for(yy in 1:length(YEARs)){ + YEAR <- YEARs[yy] + # Simulation timing, yyyy-mm-dd HH:MM:SS (UTC) + t_start <- paste0(YEAR,'-01-01 00:00:00') + t_end <- paste0(YEAR,'-12-31 23:00:00') + run_times <- seq(from = as.POSIXct(t_start, tz = 'UTC'), to = as.POSIXct(t_end, tz = 'UTC'), by = 'hour') + if(!is.null(UThrs)){ + sel <- substring(run_times,12,13)%in%UThrs #only analyze afternoon hours, following Foster papers + run_times <- run_times[sel] + } # if(!is.null(UThrs)){ + + # read in observations + objname <- paste0("SimCity_CH4_allsites_hrly_",YEARs[yy]) + print(paste("Reading in.....",objname)) + obs <- getr(objname,path=obsdir)[,c("Time",paste0("CH4_",SITE))] + + # skip footprint if CH4_SITE is NA + run_times.sub <- subset(run_times,run_times%in%obs$Time) + Times2rm <- obs[is.na(obs[,paste0("CH4_",SITE)]),'Time'] + sel <- run_times.sub%in%Times2rm + print(paste('Skipping',sum(sel),'hrs out of total of',length(run_times.sub),'hrs')) + run_times.sub <- run_times.sub[!sel] + +for(mm in 1:12){ + YYYYMM <- paste0(YEAR,formatC(mm,flag="0",width=2)) + run_times.tt <- run_times.sub[format(run_times.sub,"%Y%m")==YYYYMM] + footvar.sum <- 0 + N <- 0 #start counter +for(tt in 1:length(run_times.tt)){ + YYYYMMDDHH <- format(run_times.tt[tt],"%Y%m%d%H") + sel <- substring(XFILES,1,10)==YYYYMMDDHH + ncfilenm <- paste0(outputdir,"/",XFILES[sel],"/",XFILES[sel],"_foot.nc") + if(!file.exists(ncfilenm)){print(paste("no footprint output for:",YYYYMMDDHH));next} + # retrieve footprint + footfile <- nc_open(ncfilenm) + # print(paste(N,"Reading in:",ncfilenm)) + foot.all <- ncvar_get(footfile,"foot") + foot <- apply(foot.all,c(1,2),sum) #sum over back times + well <- well.all[,,YYYYMM] + if((nrow(well)!=nrow(foot))|(ncol(well)!=ncol(foot))){print("dimensions of foot not match those of wellinfo; skip!");next} + # convolve footprint with well info + footvar <- foot*well + footvar.sum <- footvar.sum + footvar + N <- N+1 + print(paste(N,"Processed",YYYYMMDDHH)) +} # for(tt in 1:length(run_times.tt)){ + RESULT.MON[,,YYYYMM] <- footvar.sum/N + saveRDS(RESULT.MON,finalresultname.gridded) + print(paste(finalresultname.gridded,"saved")) + gc() +} # for(mm in 1:12){ + gc() +} # for(yy in 1:length(YEARs)){ +} # for(vv in 1:length(VARS)){ + + +if(FALSE){ + require(fields) + + SITE <- "CSP" + if(SITE=="HPL"){LAT <- 40.1434; LON <- -109.4680; zagl <- 4.06} + if(SITE=="ROO"){LAT <- 40.2941; LON <- -110.0090; zagl <- 4.06} + if(SITE=="CSP"){LAT <- 40.0509; LON <- -110.0194; zagl <- 4.06} + + var <- c("Mcf.Gas","Bbls.Oil","Gas.Well","Oil.Well","Gas.Well_Producing","Oil.Well_Producing")[4] + + #YYYYMM <- "201604" + YYYYMM <- "202004" + resultnm <- paste0(SITE,"_foot_",var,"_gridded.rds") + dat.all <- readRDS(resultnm) + + dat <- dat.all[,,YYYYMM] + dat[dat==0] <- NA + dat.log10 <- log10(dat) + zlims <- NULL + if(var=="Mcf.Gas")zlims <- c(-12,-5) + if(var=="Bbls.Oil")zlims <- c(-12,-5) + flon <- as.numeric(rownames(dat)) + flat <- as.numeric(colnames(dat)) + dev.new();par(cex.lab=1.3,cex.axis=1.3,cex.main=1.3) + image.plot(flon,flat,dat.log10,main=paste("Site=",SITE," ",var,"\n",YYYYMM),zlim=zlims) + points(LON,LAT,pch=17,col="black") + + +} # if(FALSE){ diff --git a/dCH4_vs_foot_wellinfo.r b/dCH4_vs_foot_wellinfo.r new file mode 100644 index 0000000..8649cfd --- /dev/null +++ b/dCH4_vs_foot_wellinfo.r @@ -0,0 +1,524 @@ +# Plot dCH4 as function of footprint-convolved well # or footprint-convolved CH4 production +# y-axis is in [ppm], and x-axis is either +# a) [well # x ppm/(umole/m2/s)] or +# b) [umole CH4 produced x ppm/(umole/m2/s)] +# +# So the SLOPE from the linear regression has units of +# a) [(umole/m2s)/(well/m2/day)] or +# b) [(umole/m2/s)/(umole CH4 produced/m2/day)], which is a direct estimate of leak rate! +# +# This allows one to scale emissions to the Basin level using a) well # or b) CH4 produced +# b) is a direct measure of the overall leak rate! +# +# TO DO: see whether the slope from b) (leak rate) varies from year to year +# TO DO: see whether a strong linear relationship that emerges? +# +# V2(210403): add option to fill in gaps in FRU time series +# V3(210405): add calculation of leak from slope in each plot +# V4(210416): filter out days based on within 45o of observed wind direction (filterUV.TF), use Model-II regression, bootstrap to determine slope errs +# V5(210427): introduce 'mettype' to distinguish between met files driving STILT (to assess transport error)--e.g., "HRRR", "NARR", "NAM12", "WRF27" +# V6(210523): prepare figures for manuscript submission-e.g., Greek symbols, subscript, trim outliers, report N=? +# March 30th, 2021 by John C. Lin (John.Lin@utah.edu) + +################ +#SITE<-"HPL" +#SITE<-"ROO" +SITE<-"CSP" +# subset of months +MONsub<-4:9 +# subset of hours [UTC]; if NULL, then all hours +UThrs <- c(20,21,22,23) # [UTC] only analyze afternoon hours, following Foster papers + +# meteorology driving STILT +mettype <- "HRRR" +# mettype <- "NAM12" +# mettype <- "WRF27" + +#datname <- paste0(SITE,'_CH4_UVwind_foot.wellinfo_ALL.rds') # output from 'merge_CH4_UVwind_wellinfo.r' +datname <- paste0(SITE,'_CH4_UVwind_foot.wellinfo_ALL_',mettype,'.rds') # output from 'merge_CH4_UVwind_wellinfo.r' +fillFRU.TF <- TRUE # whether or not to fill in gaps in background (FRU) time series +filterUV.TF <- TRUE # filter times based on U/V (filter out times when HRRR is off--i.e., large transport errors) ? +CH4.VOLFRAC <- 0.89 # volume fraction of CH4 in natural gas of 0.89 [Karion et al. 2013] + +col.gas <- "#B91C1C" # red +col.oil <- "#1D4ED8" # blue +################ + +require(ggplot2); require(reshape2); require(lmodel2) + +if(SITE=="HPL"){LAT <- 40.1434; LON <- -109.4680; zagl <- 4.06; WINDsite <- "UBHSP"; WINDsite.2015 <- "A1633"} # UBHSP site not available in 2015, so use A1633 (RedWash) site instead in 2015 +if(SITE=="ROO"){LAT <- 40.2941; LON <- -110.0090; zagl <- 4.06; WINDsite <- "QRS"} +if(SITE=="CSP"){LAT <- 40.0509; LON <- -110.0194; zagl <- 4.06; WINDsite <- "UBCSP"} + +DAT.all <- readRDS(datname) + +# filter subset of dataset based on Months & UThrs +YYYYMMDDHH <- format(DAT.all$Time,'%Y%m%d%H') +DAT.all <- data.frame(YYYYMMDDHH, DAT.all) +sel <- substring(DAT.all$YYYYMMDDHH, 5, 6)%in%formatC(MONsub,width=2,flag='0') +sel <- sel&(substring(DAT.all$YYYYMMDDHH, 9, 10)%in%formatC(UThrs,width=2,flag='0')) +DAT.all <- DAT.all[sel,] + +if(fillFRU.TF){ + dat <- DAT.all + YYYYMM <- substring(dat$YYYYMMDDHH,1,6) + dat <- data.frame(YYYYMM,dat) + # calculate % of NAs in each month + N <- tapply(dat$YYYYMM,dat$YYYYMM,length) + sumNA <- function(x)return(sum(is.na(x))) + N.NA <- tapply(dat$CH4_FRU,dat$YYYYMM,sumNA) + print("% of data with NAs at FRU in each month:") + print(round(100*N.NA/N,2)) + + # calculate CV (coefficient of variation)--what % is variability in FRU compared to magnitude of CH4 enhancement? + xsigma <- sqrt(tapply(dat$CH4_FRU,dat$YYYYMM,var,na.rm=T)) + dCH4.ave <- tapply(dat[,paste0('CH4_',SITE)]-dat$CH4_FRU,dat$YYYYMM,mean,na.rm=T) + print('% stdev(FRU.hrly)/mean(SITE-FRU):') + print(round(100*xsigma/dCH4.ave,2)) + # calculate CV for DAILY-AVERAGED CH4 (for selected UThrs) + FRU.day <- tapply(dat$CH4_FRU,substring(dat$YYYYMMDDHH,1,8),mean,na.rm=T) + dCH4.day <- tapply(dat[,paste0('CH4_',SITE)]-dat$CH4_FRU,substring(dat$YYYYMMDDHH,1,8),mean,na.rm=T) + YYYYMMDD <- names(FRU.day) + xsigma <- sqrt(tapply(FRU.day,substring(YYYYMMDD,1,6),var,na.rm=T)) + dCH4.ave <- tapply(dCH4.day,substring(YYYYMMDD,1,6),mean,na.rm=T) + print('% stdev(FRU.day)/mean(SITE-FRU):') + print(round(100*xsigma/dCH4.ave,2)) + + # fill in gaps in FRU time series with monthly average + print("fill in gaps in FRU time series with monthly average") + f <- function(x){ + x.ave <- mean(x,na.rm=T) + x[is.na(x)] <- x.ave + return(x) + } # f <- function(x){ + FRU.gapfilled <- unlist(tapply(dat$CH4_FRU,dat$YYYYMM,f)) + if(FALSE){ + FRU.monave <- tapply(dat$CH4_FRU,dat$YYYYMM,mean,na.rm=T) + tt <- as.numeric(substring(names(FRU.monave),1,4))+as.numeric(substring(names(FRU.monave),5,6))/12 + dev.new(); plot(tt,FRU.monave,pch=16,type="o") + dev.new(); plot(dat[,c("Time","CH4_FRU")],pch=16) + isNA <- is.na(dat$CH4_FRU) + points(dat$Time[isNA],FRU.gapfilled[isNA],pch=16,col="orange") + } #if(FALSE){ + dat$CH4_FRU <- FRU.gapfilled + DAT.all <- dat +} # if(fillFRU.TF){ + +dCH4 <- DAT.all[,paste0("CH4_",SITE)] - DAT.all[,"CH4_FRU"] +DAT.all <- data.frame(dCH4,DAT.all) +rownames(DAT.all) <- DAT.all$YYYYMMDDHH + +if(filterUV.TF){ + dat <- DAT.all + Usim <- dat[,paste0("Usim.",WINDsite)]; Vsim <- dat[,paste0("Vsim.",WINDsite)] + Uobs <- dat[,paste0("Uobs.",WINDsite)]; Vobs <- dat[,paste0("Vobs.",WINDsite)] + if(SITE=="HPL"){ + # UBHSP site not available in 2015, so use A1633 (RedWash) site instead in 2015 + sel <- substring(as.character(dat$Time),1,4)=="2015" + Usim[sel] <- dat[sel,paste0("Usim.",WINDsite.2015)] + Vsim[sel] <- dat[sel,paste0("Vsim.",WINDsite.2015)] + Uobs[sel] <- dat[sel,paste0("Uobs.",WINDsite.2015)] + Vobs[sel] <- dat[sel,paste0("Vobs.",WINDsite.2015)] + } # if(SITE=="HPL"){ + isNA <- is.na(Usim)|is.na(Vsim)|is.na(Uobs)|is.na(Vobs) + # a) U/V filter: both simulated U & V have to have the same sign as the observed (i.e., wind vector has to be same QUADRANT as observed) + # SEL <- sign(Usim)!=sign(Uobs) + # SEL <- sel|(sign(Vsim)!=sign(Vobs)) + + # b) U/V filter: filter out days based on +/-45o of observed (better than quadrant method above, since not throw way data close to quadrant boundaries + WSPD.sim <- sqrt(Usim^2 + Vsim^2); WSPD.obs <- sqrt(Uobs^2 + Vobs^2) + print("Observed windspeeds (quantile):") + print(signif(quantile(WSPD.obs,c(0,0.05,0.1,0.25,0.5,0.75,1.0),na.rm=T),3)) + # see: http://tornado.sfsu.edu/geosciences/classes/m430/Wind/WindDirection.html + WDIR.sim <- (180/pi)*atan2(y=-1*Vsim,x=-1*Usim) # wind direction (-180o to +180o); direction wind is blowing FROM + WDIR.obs <- (180/pi)*atan2(y=-1*Vobs,x=-1*Uobs) # wind direction (-180o to +180o); direction wind is blowing FROM + dWDIR <- WDIR.sim - WDIR.obs + sel <- abs(dWDIR) > 180 + dWDIR[sel&!isNA] <- sign(dWDIR[sel&!isNA])*(360 - abs(dWDIR[sel&!isNA])) + SEL <- abs(dWDIR) > 45 # simulated wind direction has to be within 45o of observed + SEL[WSPD.obs < 1.0] <- FALSE # not remove data if windspeed is too low (wind anemometer not reliable) + + print(paste("Out of",sum(!isNA),"non-NA U/V times,",sum(SEL[!isNA]),"failed to pass U/V filter: ",signif(100*sum(SEL[!isNA])/sum(!isNA),3),"%")) + # assign NAs to dCH4 when failed to pass U/V filter + print(paste("Before filtering, dCH4 has",sum(!is.na(dat$dCH4)),"non-NA values")) + dat$dCH4[SEL] <- NA + print(paste(" After filtering, dCH4 has",sum(!is.na(dat$dCH4)),"non-NA values")) + DAT.all <- dat +} # if(filterUV.TF){ + +COLNMS <- c('Mcf.Gas','Bbls.Oil','Gas.Well','Oil.Well','Gas.Well_Producing','Oil.Well_Producing')[-c(3,4)] +DAT0 <- DAT.all[,c("dCH4",COLNMS)] + +slope2leak <- function(dydx,CH4.VOLFRAC=0.89) { + # calculate leak rate based on regression slope between dCH4 and Mcf.Gas + dydx <- dydx*(3600*24) # [umole*day]/[MCFnatgas*s] => [umole]/[MCFnatgas] + dydx <- dydx/CH4.VOLFRAC # [umole]/MCFnatgas] => [umole]/[MCF CH4] + dydx <- dydx/1000 # [umole]/[MCF CH4] => [umole]/[ft^3 CH4] + dydx <- dydx/0.0283 # [umole]/[ft^3 CH4] => [umole]/[m^3 CH4] + R.CH4<-8.3143*1000/16.043 #Ideal Gas constant of CH4 [J/K/kg]: (Rg/mCH4), where Rg is universal gas constant & mCH4 is molar mass of CH4 + rho.CH4<-(101300)/(R.CH4*288.7) #density of CH4 [kg/m3], using P=101.3kPa and T=288.7K [Karion et al. 2013] + dydx <- dydx/rho.CH4 # [umole]/[m^3 CH4] => [umole]/[kg CH4] + dydx <- dydx/1000 # [umole]/[kg CH4] => [umole]/[g CH4] + dydx <- dydx*1E-6 # [umole]/[g CH4] => [mole]/[g CH4] + dydx <- dydx*16 # [mole]/[g CH4] => [g CH4]/[g CH4] + dydx <- dydx*100 # [% CH4 emitted per CH4 in Natural Gas production] + leak <- dydx + return(leak) +} # slope2leak <- function(dydx) { + +bootlmodel2 <- function(x,y,N=1000,method="SMA") { + if(length(x)!=length(y))stop(paste("x and y need to be the same length!")) + # bootstrap to calculate errors in regression slope and intercept, based on Model-II regression + result <- matrix(NA,nrow=N,ncol=2); colnames(result) <- c("m","interc") + result <- data.frame(result) + for(i in 1:N){ + ind.s <- sample(1:length(x),replace=TRUE) + x2 <- x[ind.s] + y2 <- y[ind.s] + #xlm <- lmodel2(formula=y ~ x) + xlm <- lmodel2(formula=y2 ~ x2) + sel <- xlm$regression.results[,"Method"]==method #from Nick Murdoch's test code, appears that the SMA method (standard major axis regression) + result[i,"m"] <- xlm$regression.results[sel,"Slope"] + result[i,"interc"] <- xlm$regression.results[sel,"Intercept"] + } # for(i in 1:N){ + return(result) +} # bootlmodel2 <- (x,y,N=1000) { + +#----------------------------------------- All Selected Hours ---------------------------------------# +DAT <- DAT0 +wellvars <- COLNMS +# trim outliers in dCH4 +sel <- DAT$dCH4>quantile(DAT$dCH4,probs=0.99,na.rm=T) +DAT$dCH4[sel] <- NA +dCH4 <- DAT$dCH4 +data.table <- NULL +for(i in 1:length(wellvars)){ + Wellinfo <- wellvars[i] + # trim outliers in Wellinfo + sel <- DAT[,Wellinfo]>quantile(DAT[,Wellinfo],probs=0.99,na.rm=T) + DAT[sel&!is.na(DAT[,Wellinfo]),Wellinfo] <- NA + well <- DAT[,Wellinfo] + R <- cor(dCH4, well, use="na.or.complete") + # a) OLS regression + #xlm <- lm(dCH4 ~ well) + #m <- coef(xlm)[2]; interc <- coef(xlm)[1] + # b) Model-II regression + xlm <- lmodel2(formula=dCH4 ~ well) + sel <- xlm$regression.results[,"Method"]=="SMA" #from Nick Murdoch's test code, appears that the SMA method (standard major axis regression) + m <- xlm$regression.results[sel,"Slope"] + interc <- xlm$regression.results[sel,"Intercept"] + N <- xlm$n + tmp <- bootlmodel2(x=well,y=dCH4) + m.sd <- sqrt(var(tmp$m)) # error in regression slope from bootstrap + data.table <- rbind(data.table,data.frame(Wellinfo,R,m,m.sd,interc,N)) +} # for(i in 1:length(wellvars)){ +rownames(data.table) <- NULL +data.table$Wellinfo <- factor(data.table$Wellinfo, levels=c('Mcf.Gas', 'Gas.Well_Producing', 'Bbls.Oil', 'Oil.Well_Producing')) +levels(data.table$Wellinfo) <- c("Gas Production","Gas Well Density","Oil Production","Oil Well Density") + +dydx <- data.table$m[data.table$Wellinfo=="Gas Production"] # [ppm]/[MCF natgas day-1 m-2 * (ppm/umole/m2/s)] = [umole*day]/[MCFnatgas*s] +dydx <- slope2leak(dydx,CH4.VOLFRAC=CH4.VOLFRAC) # [umole*day]/[MCFnatgas*s] => [% CH4 emitted per CH4 in Natural Gas production] +dydx.sd <- data.table$m.sd[data.table$Wellinfo=="Gas Production"] # [ppm]/[MCF natgas day-1 m-2 * (ppm/umole/m2/s)] = [umole*day]/[MCFnatgas*s] +dydx.sd <- slope2leak(dydx.sd,CH4.VOLFRAC=CH4.VOLFRAC) # [umole*day]/[MCFnatgas*s] => [% CH4 emitted per CH4 in Natural Gas production] + +dat <- melt(DAT, id.vars='dCH4') +colnames(dat)[colnames(dat)=="variable"] <- "Wellinfo" +# adjust order of well info variable +dat$Wellinfo <- factor(dat$Wellinfo, levels=c('Mcf.Gas', 'Gas.Well_Producing', 'Bbls.Oil', 'Oil.Well_Producing')) +levels(dat$Wellinfo) <- c("Gas Production","Gas Well Density","Oil Production","Oil Well Density") + +YYYYMMDDHH <- rownames(DAT) +dev.new() +theme_set(theme_bw()) +g <- ggplot(dat, aes(x=value,y=dCH4)) + geom_point(,size=1.5) + + stat_smooth(method = "lm") + geom_abline() + facet_grid(. ~ Wellinfo) +g <- g + labs(title=paste0('SITE =',SITE,'; dCH4 vs sum(foot*wellinfo); mons: ',paste(unique(substring(YYYYMMDDHH,5,6)),collapse=","), + '\n UThrs=',paste(unique(substring(YYYYMMDDHH,9,10)),collapse=","), + '; Yrs=',paste(sort(unique(substring(YYYYMMDDHH,1,4))),collapse=",")), + caption=paste0('mettype=',mettype,'; fillFRU.TF=',fillFRU.TF,'; filterUV.TF=',filterUV.TF)) +g <- g + labs(subtitle=paste("CH4 leak rate from dCH4 vs GasProd slope:",round(dydx,2),"+/-",round(dydx.sd,2),"[%]")) +g <- g + theme(strip.text.x = element_text(size = 14, colour = "black"),plot.caption=element_text(size=12), + axis.title.y=element_text(size=16), axis.title.x=element_text(size=16), + axis.text.x=element_text(size=10), axis.text.y=element_text(size=14)) +ymax <- max(dat$dCH4,na.rm=T) +g <- g + ylim(c(0,ymax)) +g <- g + geom_text(data = data.table, hjust = 0, col='blue', + aes(x = 0.0, y = 0.8*ymax, label = paste0('R=', round(R, 2), '\n', + 'Slope=', signif(m, 3),'+/-',signif(m.sd,3), '\n','N=',N))) +g + facet_wrap(.~Wellinfo, nrow = 2, scales='free') +ggsave(paste0(SITE,'_dCH4_vs_foot_wellinfo_hrly_',mettype,'.png')) + + +#----------------------------------------- Average Selected Hours within a Day ---------------------------------------# +YYYYMMDD <- substring(DAT.all$YYYYMMDDHH,1,8) +dat <- DAT.all[,c('dCH4',COLNMS)] +dum <- NULL +for(cc in 1:ncol(dat)){ + tmp <- tapply(dat[,cc],YYYYMMDD,mean,na.rm=T) + dum <- cbind(dum,tmp) +} # for(cc in 1:ncol(dat)){ +colnames(dum) <- c('dCH4',COLNMS) +dum[is.nan(dum)] <- NA + +DAT <- data.frame(dum) +wellvars <- COLNMS +# trim outliers in dCH4 +sel <- DAT$dCH4>quantile(DAT$dCH4,probs=0.99,na.rm=T) +DAT$dCH4[sel] <- NA +dCH4 <- DAT$dCH4 +data.table <- NULL +for(i in 1:length(wellvars)){ + Wellinfo <- wellvars[i] + # trim outliers in Wellinfo + sel <- DAT[,Wellinfo]>quantile(DAT[,Wellinfo],probs=0.99,na.rm=T) + DAT[sel&!is.na(DAT[,Wellinfo]),Wellinfo] <- NA + well <- DAT[,Wellinfo] + R <- cor(dCH4, well, use="na.or.complete") # Pearson correlation coefficient + # Model-II regression + xlm <- lmodel2(formula=dCH4 ~ well) + sel <- xlm$regression.results[,"Method"]=="SMA" #from Nick Murdoch's test code, appears that the SMA method (standard major axis regression) + m <- xlm$regression.results[sel,"Slope"] + interc <- xlm$regression.results[sel,"Intercept"] + N <- xlm$n + tmp <- bootlmodel2(x=well,y=dCH4) + m.sd <- sqrt(var(tmp$m)) # error in regression slope from bootstrap + data.table <- rbind(data.table,data.frame(Wellinfo,R,m,m.sd,interc,N)) +} # for(i in 1:length(wellvars)){ +rownames(data.table) <- NULL +data.table$Wellinfo <- factor(data.table$Wellinfo, levels=c('Mcf.Gas', 'Gas.Well_Producing', 'Bbls.Oil', 'Oil.Well_Producing')) +levels(data.table$Wellinfo) <- c("Gas Production","Gas Well Density","Oil Production","Oil Well Density") + +dydx <- data.table$m[data.table$Wellinfo=="Gas Production"] # [ppm]/[MCF natgas day-1 m-2 * (ppm/umole/m2/s)] = [umole*day]/[MCFnatgas*s] +dydx <- slope2leak(dydx,CH4.VOLFRAC=CH4.VOLFRAC) # [umole*day]/[MCFnatgas*s] => [% CH4 emitted per CH4 in Natural Gas production] +dydx.sd <- data.table$m.sd[data.table$Wellinfo=="Gas Production"] # [ppm]/[MCF natgas day-1 m-2 * (ppm/umole/m2/s)] = [umole*day]/[MCFnatgas*s] +dydx.sd <- slope2leak(dydx.sd,CH4.VOLFRAC=CH4.VOLFRAC) # [umole*day]/[MCFnatgas*s] => [% CH4 emitted per CH4 in Natural Gas production] + +dat <- melt(DAT, id.vars='dCH4') +colnames(dat)[colnames(dat)=="variable"] <- "Wellinfo" +# adjust order of well info variable +dat$Wellinfo <- factor(dat$Wellinfo, levels=c('Mcf.Gas', 'Gas.Well_Producing', 'Bbls.Oil', 'Oil.Well_Producing')) +levels(dat$Wellinfo) <- c("Gas Production","Gas Well Density","Oil Production","Oil Well Density") +saveRDS(dat,file="dat.rds");print("dat.rds written out") +YYYYMMDD <- rownames(DAT); isNA <- is.na(DAT$dCH4) + +dev.new() +#options(scipen=-1) # tend to use Scientific Notation +theme_set(theme_bw()) +g <- ggplot(dat, aes(x=value,y=dCH4)) + geom_point(,size=1.2) + facet_wrap(. ~ Wellinfo) + + stat_smooth(method = "lm") + geom_abline() # + facet_wrap(. ~ Wellinfo) +#g <- g + labs(title=paste0('SITE =',SITE,'; dCH4 vs sum(foot*wellinfo); mons: ',paste(unique(substring(YYYYMMDD[!isNA],5,6)),collapse=","), +# '\n mean(UThrs=',paste(unique(substring(DAT.all$YYYYMMDDHH,9,10)),collapse=","), +# '); Yrs=',paste(sort(unique(substring(YYYYMMDD[!isNA],1,4))),collapse=","), +# '\n CH4 leak rate from dCH4 vs GasProd slope:',round(dydx,2),"+/-",round(dydx.sd,2),"[%]")) +MONuniq <- paste(unique(substring(YYYYMMDD[!isNA],5,6)),collapse=",") +HHuniq <- paste(unique(substring(DAT.all$YYYYMMDDHH,9,10)),collapse=",") +YYYYuniq <- paste(sort(unique(substring(YYYYMMDD[!isNA],1,4))),collapse=",") +g <- g + labs(title=substitute(paste(SITE," daily; ",Delta,"CH"[4]," vs GasProd slope(leak rate)=",dydx,"+/-",dydx.sd,"[%]"), + list(SITE=SITE,HHuniq=HHuniq,dydx=round(dydx,2),dydx.sd=round(dydx.sd,2)))) +g <- g + labs(caption=paste0('mettype=',mettype,'; fillFRU.TF=',fillFRU.TF,'; filterUV.TF=',filterUV.TF)) +g <- g + scale_x_continuous(labels = scales::scientific) # forces use of scientific notatiion on x-axis +g <- g + theme(strip.text.x = element_text(size = 14, colour = "black"), axis.title.y=element_text(size=16), axis.title.x=element_text(size=16), + axis.text.x=element_text(size=11.0), axis.text.y=element_text(size=14),plot.caption=element_text(size=12),plot.title=element_text(size=15.5)) +ymax <- max(dat$dCH4,na.rm=T) +g <- g + ylim(c(0,ymax)) +g <- g + geom_text(data = data.table, hjust = 0, col='blue', + aes(x = 0.0, y = 0.8*ymax, label = paste0('R=', round(R, 2), '\n', + 'Slope=', signif(m, 3),'+/-',signif(m.sd,3), '\n','N=',N))) +g <- g + labs(y=expression(paste(Delta,"CH4 [ppm]"))) +g + facet_wrap(.~Wellinfo, nrow = 2, scales='free') +ggsave(paste0(SITE,'_dCH4_vs_foot_wellinfo_daily_',mettype,'.png')) + + +# Multiple plot function +# Copied from http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/ +# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects) +# - cols: Number of columns in layout +# - layout: A matrix specifying the layout. If present, 'cols' is ignored. +# +# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE), +# then plot 1 will go in the upper left, 2 will go in the upper right, and +# 3 will go all the way across the bottom. +# +multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) { + library(grid) + # Make a list from the ... arguments and plotlist + plots <- c(list(...), plotlist) + numPlots = length(plots) + # If layout is NULL, then use 'cols' to determine layout + if (is.null(layout)) { + # Make the panel + # ncol: Number of columns of plots + # nrow: Number of rows needed, calculated from # of cols + layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), + ncol = cols, nrow = ceiling(numPlots/cols)) + } + if (numPlots==1) { + print(plots[[1]]) + } else { + # Set up the page + grid.newpage() + pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout)))) + + # Make each plot, in the correct location + for (i in 1:numPlots) { + # Get the i,j matrix positions of the regions that contain this subplot + matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE)) + print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, layout.pos.col = matchidx$col)) + } + } +} + + +# Plot each panel separately, in order to more closely control plotting and to enable different x-axis labels +theme_set(theme_bw()) +wellvar <- "Gas Production" +g <- ggplot(subset(dat,Wellinfo==wellvar), aes(x=value,y=dCH4)) + geom_point(,size=1.2) +SEL <- data.table$Wellinfo==wellvar +g <- g + geom_abline(slope=data.table$m[SEL],intercept=data.table$interc[SEL],colour="blue") +g <- g + labs(title=paste0(SITE,": ",wellvar)) +if(SITE=="HPL")g <- g + scale_x_continuous(labels = scales::scientific) # forces use of scientific notatiion on x-axis +g <- g + theme(strip.text.x = element_text(size = 14, colour = "black"), axis.title.y=element_text(size=14), axis.title.x=element_text(size=14), + axis.text.x=element_text(size=10.0), axis.text.y=element_text(size=14),plot.caption=element_text(size=12),plot.title=element_text(size=14.0)) +ymax <- max(dat$dCH4,na.rm=T) +g <- g + ylim(c(0,ymax)) +g <- g + geom_text(data = data.table[SEL,], hjust = 0, col='blue', aes(x = 0.0, y = 0.8*ymax, label = paste0('R=', round(R, 2), '\n', 'Slope=', signif(m, 3),'+/-',signif(m.sd,3), '\n','N=',N)),size=4) +xpos <- 1.2E-4 +if(SITE=="CSP")xpos <- 5.0E-5 +g <- g + geom_text( col=col.gas, aes(x = xpos, y = ymax, label = paste0("Leakage rate: ",round(dydx,2),"+/-",round(dydx.sd,2)," [%]")),size=4.5) +g1 <- g + labs(y=expression(paste(Delta,"CH"[4]," [ppm]")), x=expression(paste("[Mcf day"^-1," ppm/(",mu,"mole s"^-1,")]"))) + +theme_set(theme_bw()) +wellvar <- "Oil Production" +dat.sub <- dat[dat$Wellinfo==wellvar,] +g <- ggplot(subset(dat,Wellinfo==wellvar), aes(x=value,y=dCH4)) + geom_point(,size=1.2) +SEL <- data.table$Wellinfo==wellvar +g <- g + geom_abline(slope=data.table$m[SEL],intercept=data.table$interc[SEL],colour="blue") +g <- g + labs(title=paste0(SITE,": ",wellvar)) +g <- g + scale_x_continuous(labels = scales::scientific) # forces use of scientific notatiion on x-axis +g <- g + theme(strip.text.x = element_text(size = 14, colour = "black"), axis.title.y=element_text(size=14), axis.title.x=element_text(size=14), + axis.text.x=element_text(size=10.5), axis.text.y=element_text(size=14),plot.caption=element_text(size=12),plot.title=element_text(size=14.0)) +ymax <- max(dat$dCH4,na.rm=T) +g <- g + ylim(c(0,ymax)) +g <- g + geom_text(data = data.table[SEL,], hjust = 0, col='blue', aes(x = 0.0, y = 0.8*ymax, label = paste0('R=', round(R, 2), '\n', 'Slope=', signif(m, 3),'+/-',signif(m.sd,3), '\n','N=',N)),size=4) +g2 <- g + labs(y=expression(paste(Delta,"CH"[4]," [ppm]")), x=expression(paste("[Barrels day"^-1," ppm/(",mu,"mole s"^-1,")]"))) + +theme_set(theme_bw()) +wellvar <- "Gas Well Density" +dat.sub <- dat[dat$Wellinfo==wellvar,] +g <- ggplot(subset(dat,Wellinfo==wellvar), aes(x=value,y=dCH4)) + geom_point(,size=1.2) +SEL <- data.table$Wellinfo==wellvar +g <- g + geom_abline(slope=data.table$m[SEL],intercept=data.table$interc[SEL],colour="blue") +g <- g + labs(title=paste0(SITE,": ",wellvar)) +g <- g + scale_x_continuous(labels = scales::scientific) # forces use of scientific notatiion on x-axis +g <- g + theme(strip.text.x = element_text(size = 14, colour = "black"), axis.title.y=element_text(size=14), axis.title.x=element_text(size=14), + axis.text.x=element_text(size=10.5), axis.text.y=element_text(size=14),plot.caption=element_text(size=12),plot.title=element_text(size=14.0)) +ymax <- max(dat$dCH4,na.rm=T) +g <- g + ylim(c(0,ymax)) +g <- g + geom_text(data = data.table[SEL,], hjust = 0, col='blue', aes(x = 0.0, y = 0.8*ymax, label = paste0('R=', round(R, 2), '\n', 'Slope=', signif(m, 3),'+/-',signif(m.sd,3), '\n','N=',N)),size=4) +g3 <- g + labs(y=expression(paste(Delta,"CH"[4]," [ppm]")), x=expression(paste("[Well # day"^-1," ppm/(",mu,"mole s"^-1,")]"))) + +theme_set(theme_bw()) +wellvar <- "Oil Well Density" +dat.sub <- dat[dat$Wellinfo==wellvar,] +g <- ggplot(subset(dat,Wellinfo==wellvar), aes(x=value,y=dCH4)) + geom_point(,size=1.2) +SEL <- data.table$Wellinfo==wellvar +g <- g + geom_abline(slope=data.table$m[SEL],intercept=data.table$interc[SEL],colour="blue") +g <- g + labs(title=paste0(SITE,": ",wellvar)) +g <- g + scale_x_continuous(labels = scales::scientific) # forces use of scientific notatiion on x-axis +g <- g + theme(strip.text.x = element_text(size = 14, colour = "black"), axis.title.y=element_text(size=14), axis.title.x=element_text(size=14), + axis.text.x=element_text(size=10.5), axis.text.y=element_text(size=14),plot.caption=element_text(size=12),plot.title=element_text(size=14.0)) +ymax <- max(dat$dCH4,na.rm=T) +g <- g + ylim(c(0,ymax)) +g <- g + geom_text(data = data.table[SEL,], hjust = 0, col='blue', aes(x = 0.0, y = 0.8*ymax, label = paste0('R=', round(R, 2), '\n', 'Slope=', signif(m, 3),'+/-',signif(m.sd,3), '\n','N=',N)),size=4) +g4 <- g + labs(y=expression(paste(Delta,"CH"[4]," [ppm]")), x=expression(paste("[Well # day"^-1," ppm/(",mu,"mole s"^-1,")]"))) + +multiplot(g1,g2,g3,g4,cols=2) +#ggsave(paste0(SITE,'_dCH4_vs_foot_wellinfo_daily_',mettype,'_multiplot.png')) #ggsave only saves one of the panels, so use dev.copy instead... +#dev.copy(png,paste0(SITE,'_dCH4_vs_foot_wellinfo_daily_',mettype,'_multiplot.png'));dev.off() +dev.copy(pdf,paste0(SITE,'_dCH4_vs_foot_wellinfo_daily_',mettype,'_multiplot.pdf'));dev.off() + + +#----------------------------------------- Analyze Year-by-Year, and Average Selected Hours within a Day ---------------------------------------# +if(FALSE){ +YEARs <- unique(substring(DAT.all$YYYYMMDD,1,4)) + +for(yy in 1:length(YEARs)){ + +YEAR <- YEARs[yy] +sel <- substring(DAT.all$YYYYMMDD,1,4)==YEAR +YYYYMMDD <- substring(DAT.all$YYYYMMDDHH[sel],1,8) +dat <- DAT.all[sel,c('dCH4',COLNMS)] +dum <- NULL +for(cc in 1:ncol(dat)){ + tmp <- tapply(dat[,cc],YYYYMMDD,mean,na.rm=T) + dum <- cbind(dum,tmp) +} # for(cc in 1:ncol(dat)){ +colnames(dum) <- c('dCH4',COLNMS) +dum[is.nan(dum)] <- NA + +DAT <- data.frame(dum) +wellvars <- COLNMS +# trim outliers in dCH4 +sel <- DAT$dCH4>quantile(DAT$dCH4,probs=0.99,na.rm=T) +DAT$dCH4[sel] <- NA +dCH4 <- DAT$dCH4 +data.table <- NULL +for(i in 1:length(wellvars)){ + Wellinfo <- wellvars[i] + # trim outliers in Wellinfo + sel <- DAT[,Wellinfo]>quantile(DAT[,Wellinfo],probs=0.99,na.rm=T) + DAT[sel&!is.na(DAT[,Wellinfo]),Wellinfo] <- NA + well <- DAT[,Wellinfo] + R <- cor(dCH4, well, use="na.or.complete") + # a) OLS regression + #xlm <- lm(dCH4 ~ well) + #m <- coef(xlm)[2]; interc <- coef(xlm)[1] + # b) Model-II regression + xlm <- lmodel2(formula=dCH4 ~ well) + sel <- xlm$regression.results[,"Method"]=="SMA" #from Nick Murdoch's test code, appears that the SMA method (standard major axis regression) + m <- xlm$regression.results[sel,"Slope"] + interc <- xlm$regression.results[sel,"Intercept"] + N <- xlm$n + tmp <- bootlmodel2(x=well,y=dCH4) + m.sd <- sqrt(var(tmp$m)) # error in regression slope from bootstrap + data.table <- rbind(data.table,data.frame(Wellinfo,R,m,m.sd,interc,N)) +} # for(i in 1:length(wellvars)){ +rownames(data.table) <- NULL +data.table$Wellinfo <- factor(data.table$Wellinfo, levels=c('Mcf.Gas', 'Gas.Well_Producing', 'Bbls.Oil', 'Oil.Well_Producing')) +levels(data.table$Wellinfo) <- c("Gas Production","Gas Well Density","Oil Production","Oil Well Density") + +dydx <- data.table$m[data.table$Wellinfo=="Gas Production"] # [ppm]/[MCF natgas day-1 m-2 * (ppm/umole/m2/s)] = [umole*day]/[MCFnatgas*s] +dydx <- slope2leak(dydx,CH4.VOLFRAC=CH4.VOLFRAC) # [umole*day]/[MCFnatgas*s] => [% CH4 emitted per CH4 in Natural Gas production] +dydx.sd <- data.table$m.sd[data.table$Wellinfo=="Gas Production"] # [ppm]/[MCF natgas day-1 m-2 * (ppm/umole/m2/s)] = [umole*day]/[MCFnatgas*s] +dydx.sd <- slope2leak(dydx.sd,CH4.VOLFRAC=CH4.VOLFRAC) # [umole*day]/[MCFnatgas*s] => [% CH4 emitted per CH4 in Natural Gas production] + +dat <- melt(DAT, id.vars='dCH4') +colnames(dat)[colnames(dat)=="variable"] <- "Wellinfo" +# adjust order of well info variable +dat$Wellinfo <- factor(dat$Wellinfo, levels=c('Mcf.Gas', 'Gas.Well_Producing', 'Bbls.Oil', 'Oil.Well_Producing')) +levels(dat$Wellinfo) <- c("Gas Production","Gas Well Density","Oil Production","Oil Well Density") + +YYYYMMDD <- rownames(DAT); isNA <- is.na(DAT$dCH4) +dev.new() +theme_set(theme_bw()) +g <- ggplot(dat, aes(x=value,y=dCH4)) + geom_point(,size=1.5) + + stat_smooth(method = "lm") + geom_abline() + facet_grid(. ~ Wellinfo) +g <- g + labs(title=paste0('SITE =',SITE,'; dCH4 vs sum(foot*wellinfo); mons: ',paste(unique(substring(YYYYMMDD[!isNA],5,6)),collapse=","), + '\n mean(UThrs=',paste(unique(substring(DAT.all$YYYYMMDDHH,9,10)),collapse=","), + '); Yrs=',paste(sort(unique(substring(YYYYMMDD[!isNA],1,4))),collapse=",")), + caption=paste0('mettype=',mettype,'; fillFRU.TF=',fillFRU.TF,'; filterUV.TF=',filterUV.TF)) +g <- g + labs(subtitle=paste("CH4 leak rate from dCH4 vs GasProd slope:",round(dydx,2),"+/-",round(dydx.sd,2),"[%]")) +g <- g + theme(strip.text.x = element_text(size = 14, colour = "black"), axis.title.y=element_text(size=16), axis.title.x=element_text(size=16), + axis.text.x=element_text(size=10), axis.text.y=element_text(size=14),plot.caption=element_text(size=12)) +ymax <- max(dat$dCH4,na.rm=T) +g <- g + ylim(c(0,ymax)) +g <- g + geom_text(data = data.table, hjust = 0, col='blue', + aes(x = 0.0, y = 0.8*ymax, label = paste0('R=', round(R, 2), '\n', + 'Slope=', signif(m, 3),'+/-',signif(m.sd,3), '\n','N=',N))) +g + facet_wrap(.~Wellinfo, nrow = 2, scales='free') +ggsave(paste0(SITE,'_dCH4_vs_foot_wellinfo_daily_',YEAR,'_',mettype,'.png')) + +} # for(yy in 1:length(YEARs)){ + +} #if(FALSE){ + diff --git a/export_Fch4simple.r b/export_Fch4simple.r new file mode 100644 index 0000000..5956bbf --- /dev/null +++ b/export_Fch4simple.r @@ -0,0 +1,33 @@ +# Export merged dataset from 3 different meteorological datasets +# used to calculate CH4 emissions +# May 26th, 2021 by John C. Lin (John.Lin@utah.edu) + +################### +obsdir <- "/uufs/chpc.utah.edu/common/home/lin-group4/jcl/SimCity" +mettypes <- c("HRRR","NAM12","WRF27") +################### + +#---------- export observational datasets ----------# +sites <- c("FRU","HPL","CSP") +YEARs <- 2015:2020 +tracer <- "CH4" +for(yy in 1:length(YEARs)){ + YEAR <- YEARs[yy] + # generated by "/uufs/chpc.utah.edu/common/home/lin-group4/jcl/SimCity/extract_hrly_subsetsites.r" + outputname<-paste("UUCON_",tracer,"_",paste0(sites,collapse="_"),"_hrly_",YEAR,".csv",sep="") + file.copy(from=file.path(obsdir,outputname),to=file.path("./",outputname)) +} # for(yy in 1:length(YEARs)){ +system(paste0("rename UUCON_",tracer," UINTA_",tracer," UINTA_*.csv")) +system(paste0("tar cvf ~/public_html/Exchange/UINTA_",tracer,".tar UINTA_",tracer,"*.csv")) + +#---------- export output from Fch4_simple.r ----------# +SITE <- "HPL" +for(mm in 1:length(mettypes)){ + mettype <- mettypes[mm] + dat <- readRDS(paste0("Fch4_",SITE,"_",mettype,".rds")) + CSVfile <- paste0("Fch4_",SITE,"_",mettype,".csv") + write.csv(dat,file=CSVfile,row.names=FALSE) + print(paste(CSVfile,"written out")) +} # for(mm in 1:length(mettypes)){ +system(paste0("tar cvf ~/public_html/Exchange/Fch4.tar Fch4_",SITE,"_*.csv")) + diff --git a/export_Uintah_tseries.r b/export_Uintah_tseries.r new file mode 100644 index 0000000..4976636 --- /dev/null +++ b/export_Uintah_tseries.r @@ -0,0 +1,29 @@ +# Exports Uintah Basin CH4/CO2 time series +# December 14th, 2020 by John C. Lin (John.Lin@utah.edu) + +#################### +YEARs <- 2018:2020 +species <- c("CH4","CO2") +SITEs <- c("FRU","ROO","CSP","HPL") # sites to output +obsdir <- "/uufs/chpc.utah.edu/common/home/lin-group4/jcl/SimCity/" +#################### + +for(ss in species){ + specie <- ss +for(i in 1:length(YEARs)){ + objname <- paste0("SimCity_",specie,"_allsites_hrly_",YEARs[i]) + print(paste("Reading in.....",objname)) + tmp <- getr(objname,path=obsdir) + colnms <- colnames(tmp) + sel <- colnms %in% paste0(specie,"_",SITEs) + dat <- data.frame(Time_UTC=tmp$Time,tmp[,sel]) + + outputfilenm <- paste0("Uintah_",specie,"_",YEARs[i],".csv") + write.csv(dat,outputfilenm,row.names=FALSE) + print(paste(outputfilenm,"generated")) + gc() +} #for(i in 1:length(YEARs)){ +} #for(ss in species){ + + + diff --git a/grid_well_info_foot.r b/grid_well_info_foot.r new file mode 100644 index 0000000..1ccf562 --- /dev/null +++ b/grid_well_info_foot.r @@ -0,0 +1,149 @@ +# Creates gridded information of well numbers and well production that can be convolved with +# STILT footprints to quantitatively relate wells to receptors +# V2(210404): since fluxes are in [umole/m2/s], need to divide well quantities by area of footprint gridcell to cancel out [m-2] and number of days in a month to cancel out [s-1] +# March 22nd, 2021 by John C. Lin (John.Lin@utah.edu) + +###################### +datname <- "State_of_Utah_oil_gas_DATA/well_data_gas.oil.production_monthly.rds" # generated by 'State_of_Utah_oil_gas_DATA/merge_well_data_production.r' +# Months and Years to grid +MMs <- formatC(1:12,flag="0",width=2) +YYYYs <- 2015:2020 +YYYYMMs <- paste0(rep(YYYYs,each=length(MMs)),rep(MMs,length(YYYYs))) +# where footprints are stored +outputdir<-paste0("/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_HYSPLIT-STILT/out/HPL/by-id/") +plotTF <- FALSE +###################### + +require(ncdf4);require(fields) +require(maps);require(RgoogleMaps) +require(geosphere) + +DAT <- readRDS(datname) + +# Footprint grid settings from 'run_stilt_hysplit-stiltV1.r' +xmn <- -112.5; xmx <- -108 +ymn <- 38.5; ymx <- 41.5 +xres <- 0.01; yres <- xres + +# read in sample footprint to get footprint dimensions +YYYYMMDDs <- list.files(outputdir,pattern="201604")[1] +ncfilenm <- paste0(outputdir,"/",YYYYMMDDs[1],"/",YYYYMMDDs[1],"_foot.nc") +if(!file.exists(ncfilenm))stop(paste("no footprint output for:",ncfilenm)) +# Generate footprint +footfile <- nc_open(ncfilenm); foot.all <- ncvar_get(footfile,"foot") +footsum <- apply(foot.all,c(1,2),sum) #sum over back times +flat <- ncvar_get(footfile,"lat");flon <- ncvar_get(footfile,"lon") +# lat/lons refer to footprint grid CENTERS, so move them to lower left (southwest) corner to facilitate gridding +flat.ll <- flat - yres/2; flon.ll <- flon - xres/2 + +# determine footprint gridcell areas +flon.vec <- rep(flon,ncol(footsum)) +flat.vec <- rep(flat,each=nrow(footsum)) +dx <- distHaversine(p1=cbind(flon.vec,flat.vec),p2=cbind(flon.vec+xres,flat.vec)) +dy <- mean(distHaversine(p1=cbind(flon.vec,flat.vec),p2=cbind(flon.vec,flat.vec+yres))) # gridcell dimension in y-direction is almost a constant +AREA <- dx*dy # gridcell area [m^2] +AREA.mat <- matrix(AREA,nrow=length(flon)) #; image.plot(AREA.mat) + +# map well lat/lon to footprint grid index +foot.iy <- floor((DAT$Latitude - ymn)/yres) + 1 +foot.ix <- floor((DAT$Longitude - xmn)/yres) + 1 +DAT <- data.frame(DAT,foot.iy,foot.ix) + +# I. SUM total production in each footprint grid +vars <- c("Mcf.Gas","Bbls.Oil") +for(vv in 1:length(vars)){ + var <- vars[vv] + result <- array(data=0, dim=c(length(flon),length(flat),length(YYYYMMs))) + dimnames(result) <- list(flon, flat, YYYYMMs) + print(paste("Processing: ",var,"......")) +for(i in 1:length(YYYYMMs)){ + YYYYMM <- YYYYMMs[i] + sel <- substring(DAT$YYYYMMDD,1,6) == YYYYMM + dat <- DAT[sel,] + print(paste0(YYYYMM,"; Max # of days in month: ",max(dat$Days.Prod,na.rm=T))) + ind <- paste0(dat$foot.ix,",",dat$foot.iy) + xsum <- tapply(dat[,var],ind,sum) + tmp <- strsplit(names(xsum),",") + inds.mat <- matrix(as.numeric(unlist(tmp)),byrow=T,ncol=2) + result[,,i][inds.mat] <- xsum + result[,,i] <- result[,,i]/AREA.mat # normalize by gridcell area [m^2] + result[,,i] <- result[,,i]/max(dat$Days.Prod,na.rm=T) # normalize by number of days in month + gc() +} #for(i in 1:length(YYYYMMs)){ + resultname <- paste0(var,"_gridded.rds") + saveRDS(result,resultname) + print(paste(resultname,"generated")) + gc() +} #for(vv in 1:length(vars)){ + +# II. TALLY total well numbers of different types in each footprint grid +vars <- c("Gas Well","Oil Well","Unlabelled","Gas Well:Producing","Oil Well:Producing")[-3] +for(vv in 1:length(vars)){ + var <- vars[vv] + result <- array(data=0, dim=c(length(flon),length(flat),length(YYYYMMs))) + dimnames(result) <- list(flon, flat, YYYYMMs) + print(paste("Processing: ",var,"......")) +for(i in 1:length(YYYYMMs)){ + YYYYMM <- YYYYMMs[i] + sel <- substring(DAT$YYYYMMDD,1,6) == YYYYMM + if(sum(sel)==0){print(paste("No data:",YYYYMM,"; skip"));next} + dat <- DAT[sel,] + if(var%in%c("Gas Well","Oil Well","Unlabelled")){ + var2 <- var + if(var=="Unlabelled")var2 <- "" + sel <- dat$Well.Type.x == var2 + if(sum(sel)==0){print(paste("No data:",YYYYMM,"; ",var2,"skip"));next} + dat <- dat[sel,] + } else { + sel <- paste0(dat$Well.Type.x,":",dat$Well.Status.x)==var + dat <- dat[sel,] + } #if(){ + ind <- paste0(dat$foot.ix,",",dat$foot.iy) + N <- tapply(ind,ind,length) + tmp <- strsplit(names(N),",") + inds.mat <- matrix(as.numeric(unlist(tmp)),byrow=T,ncol=2) + result[,,i][inds.mat] <- N + result[,,i] <- result[,,i]/AREA.mat # normalize by gridcell area [m^2] + result[,,i] <- result[,,i]/max(dat$Days.Prod,na.rm=T) # normalize by number of days in month + print(paste(YYYYMM,var,": sum(N)=",sum(N),"; max # of days in month:",max(dat$Days.Prod,na.rm=T))) + gc() +} #for(i in 1:length(YYYYMMs)){ + if(length(grep(":",var))) var <- gsub(":","_",var) + resultname <- paste0(gsub(" ",".",var),"_gridded.rds") + saveRDS(result,resultname) + print(paste(resultname,"generated")) + gc() +} #for(vv in 1:length(vars)){ + + +# tally number of oil and gas producing wells +gas <- readRDS("Gas.Well_Producing_gridded.rds") +oil <- readRDS("Oil.Well_Producing_gridded.rds") +gas <- gas[,,"202009"]*AREA.mat +oil <- oil[,,"202009"]*AREA.mat +gas <- gas*30 # multiply by number of days in month +oil <- oil*30 # multiply by number of days in month +print(paste("Total of:",sum(gas),"natural gas producing wells\n",sum(oil),"oil producing wells")) + + + +if(plotTF){ + var <- "Mcf.Gas" + #var <- "Bbls.Oil" + resultname <- paste0(var,"_gridded.rds") + dat.all <- readRDS(resultname) + + # time series per year + dat <- apply(dat.all,c(3),sum) + dev.new() + YYYY <- as.numeric(substring(names(dat),1,4)) + MM <- as.numeric(substring(names(dat),5,6)) + plot(YYYY+(MM-1)/12,dat,pch=16,type="o",xlab="Year",ylab=var) + + i <- 65 + dev.new() + image.plot(x=flon,y=flat,z=dat.all[,,i],main=paste(var,"\n",dimnames(dat.all)[[3]][i])) + +} #if(plotTF){ + + diff --git a/merge_CH4_UVwind_wellinfo.r b/merge_CH4_UVwind_wellinfo.r new file mode 100644 index 0000000..6cfa04f --- /dev/null +++ b/merge_CH4_UVwind_wellinfo.r @@ -0,0 +1,79 @@ +# Merges CH4 observations and wind (observed + simulated UV) with results of convolution of footprints with well information +# March 30th, 2021 by John C. Lin (John.Lin@utah.edu) + +################# +YEARs <- 2015:2020 +#SITE<-"HPL" +#SITE<-"ROO" +SITE<-"CSP" +if(SITE=="CSP")YEARs <- c(2016,2019,2020) +if(SITE=="ROO")YEARs <- 2015:2019 + +# subset of hours [UTC]; if NULL, then all hours +# UThrs <- NULL +UThrs <- c(20,21,22,23) #[UTC] only analyze afternoon hours, following Foster papers + +# meteorology driving STILT +# mettype <- "NAM12" +# mettype <- "WRF27" +mettype <- "HRRR" + +obsdir<-"/uufs/chpc.utah.edu/common/home/lin-group4/jcl/SimCity/" #where CH4 observations reside +winddir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/Transporterr_stiltread/Output" #where wind obs and sim values are found +finalresultname <- paste0(SITE,'_CH4_UVwind_foot.wellinfo_ALL_',mettype,'.rds') +################ + +if(SITE=="HPL"){LAT <- 40.1434; LON <- -109.4680; zagl <- 4.06; WINDsite <- "UBHSP"; WINDsite.2015 <- "A1633"} # UBHSP site not available in 2015, so use A1633 (RedWash) site instead in 2015 +if(SITE=="ROO"){LAT <- 40.2941; LON <- -110.0090; zagl <- 4.06; WINDsite <- "QRS"} +if(SITE=="CSP"){LAT <- 40.0509; LON <- -110.0194; zagl <- 4.06; WINDsite <- "UBCSP"} + +obs.all <- NULL +wind.all <- NULL +dat.all <- NULL +for(yy in 1:length(YEARs)){ + # read in observations + objname <- paste0("SimCity_CH4_allsites_hrly_",YEARs[yy]) + print(paste("Reading in.....",objname)) + tmp <- getr(objname,path=obsdir)[,c("Time",paste0("CH4_",c(SITE,"FRU")))] + obs.all <- rbind(obs.all,tmp) + + # read in sum(foot*wellinfo)--i.e., output from 'convolve_foot_wellinfo.r' + #datname <- paste0(SITE,'_foot.wellinfo_',YEARs[yy],'.rds') + datname <- paste0(SITE,'_foot.wellinfo_',YEARs[yy],'_',mettype,'.rds') + print(paste("Reading in.....",datname)) + tmp.dat <- readRDS(datname) + dat.all <- rbind(dat.all,tmp.dat) + + # read in HRRR winds + observed values + windname <- paste0("HRRR_obs_",YEARs[yy],"01010000to",YEARs[yy],"12312300.RDS") + if(YEARs[yy]=="2018")windname <- paste0("HRRR_obs_",YEARs[yy],"01010000to",YEARs[yy],"12310000.RDS") # typo during 2018 + colnms <- c("Time",paste0(c("Usim.","Vsim.","Tsim.","Uobs.","Vobs.","Tobs."),WINDsite)) + # UBHSP site not available in 2015, so use A1633 (RedWash) site instead in 2015 + if(SITE=="HPL")colnms <- c("Time",paste0(c("Usim.","Vsim.","Tsim.","Uobs.","Vobs.","Tobs."),WINDsite), + paste0(c("Usim.","Vsim.","Tsim.","Uobs.","Vobs.","Tobs."),WINDsite.2015)) + if(file.exists(paste0(winddir,"/",windname))) { + print(paste("Reading in.....",windname)) + tmp0 <- readRDS(paste0(winddir,"/",windname))$dat + tmp <- matrix(NA,nrow=nrow(tmp0),ncol=length(colnms)) + colnames(tmp) <- colnms + tmp <- data.frame(tmp) + colnms.sel <- colnames(tmp0)%in%colnms + tmp[,colnames(tmp0)[colnms.sel]] <- tmp0[,colnms.sel] + } else { + print(paste("Doesn't exists:",windname,"; replaces with NAs")) + tmp <- matrix(NA,nrow=nrow(tmp.dat),ncol=length(colnms)) + colnames(tmp) <- colnms + tmp <- data.frame(tmp) + tmp$Time <- tmp.dat$Time + } # if(file.exists(paste0(winddir,"/",windname)){ + wind.all <- rbind(wind.all,tmp) + + gc() +} # for(yy in 1:length(YEARs)){ + +# merge observations with sum(foot*wellinfo) +result <- merge(x=obs.all,y=dat.all,by='Time',all.x=FALSE,all.y=TRUE) +result <- merge(x=wind.all,y=result,by='Time',all.x=FALSE,all.y=TRUE) +saveRDS(result,file=finalresultname) +print(paste(finalresultname,'written out')) + diff --git a/plot_CH4tseries_diurnal_seasonal.r b/plot_CH4tseries_diurnal_seasonal.r new file mode 100644 index 0000000..14e29a7 --- /dev/null +++ b/plot_CH4tseries_diurnal_seasonal.r @@ -0,0 +1,230 @@ +# Plots CH4 time series, diurnal, and seasonal cycle in Uintah Basin +# V2(210521): prepare figures for publicatoin: a) change colors b) use subscripts c) remove UThrs, mons labelling +# April 26th, 2021 by John C. Lin (John.Lin@utah.edu) + +#################### +UTsel <- 20:23 # [UTC] only analyze afternoon hours, following Foster papers +obsdir <- "/uufs/chpc.utah.edu/common/home/lin-group4/jcl/SimCity/" +col.gas <- "#B91C1C" # red +col.oil <- "#1D4ED8" # blue +#################### + +######################################################### +#I. Compare CH4 time series +########## +MONsel <- 4:9 # months to examine +YEARs <- 2015:2020 +SITEs <- c("FRU","CSP","HPL") +########## +# combine data from different years +dat.all<-NULL +colnms.sel <- paste0("CH4_",SITEs) +for(i in 1:length(YEARs)){ + objname<-paste0("SimCity_CH4_allsites_hrly_",YEARs[i]) + print(paste("Reading in.....",objname)) + tmp <- getr(objname,path=obsdir) + dat.add <- matrix(NA,nrow=nrow(tmp),ncol=length(colnms.sel)) + colnames(dat.add) <- colnms.sel + dat.add <- data.frame(dat.add) + + sel <- colnames(tmp)%in%colnms.sel + dat.add[,colnames(tmp)[sel]] <- tmp[,sel] + dat.add <- data.frame(Time=tmp$Time,dat.add) + dat.all<-rbind(dat.all,dat.add) + gc() +} #for(i in 1:length(YEARs)){ + +# average over UTsel +YYYYMMDDHH <- format(dat.all[,"Time"],format="%Y%m%d%H") +tmp <- dat.all[substring(YYYYMMDDHH,9,10)%in%formatC(UTsel,width=2,flag="0"),] +YYYYMMDD <- format(tmp[,"Time"],format="%Y%m%d") +CH4_FRU <- tapply(tmp$CH4_FRU,YYYYMMDD,mean,na.rm=T) +CH4_CSP <- tapply(tmp$CH4_CSP,YYYYMMDD,mean,na.rm=T) +CH4_HPL <- tapply(tmp$CH4_HPL,YYYYMMDD,mean,na.rm=T) +Time <- as.POSIXct(strptime(names(CH4_HPL),"%Y%m%d"),tz="GMT") + +#ylims <- 1000*range(c(CH4_FRU,CH4_CSP,CH4_HPL),na.rm=T) +#ylims <- c(1900,3000) +ylims <- c(1850,2500) +#ylims <- c(1850,2300) +#ylims <- NULL +dev.new(width=8,height=6);par(cex.axis=1.5,cex.main=1.5,cex.lab=1.5,mar=c(5,5,4,2)) +plot(Time,1000*CH4_HPL,pch=16,ylim=ylims,cex=0.6,col=col.gas,xlab="Year",ylab=expression("CH"[4]*" [ppb]")) +points(Time,1000*CH4_CSP,pch=16,cex=0.6,col=col.oil) +points(Time,1000*CH4_FRU,pch=16,cex=0.5,col="black") +title(main=expression("Observed Daily CH"[4]*" (afternoon-averaged)")) # \n UThrs:",paste(UTsel,collapse=","))) +# draw gray boxes for times not included in analyses +xleft <- as.POSIXct(strptime(x=paste0(c(YEARs),"-",formatC(max(MONsel)+1,flag="0",width=2),"-","01"), + format="%Y-%m-%d"),tz="GMT") +xright <- as.POSIXct(strptime(x=paste0(c(YEARs+1),"-",formatC(min(MONsel),flag="0",width=2),"-","01"), + format="%Y-%m-%d"),tz="GMT") +rect(xleft=xleft,ybottom=ylims[1]-100,xright=xright,ytop=ylims[2]+100,col=gray(level=0.3,alpha=0.5),border="darkgray") +legend("topleft",c("HPL","CSP","FRU"),col=c(col.gas,col.oil,"black"), + text.col=c(col.gas,col.oil,"black"),pch=16,bg="white",cex=1.3) +dev.copy(png,filename="CH4_tseries_allsites.png",width=480*8/6,height=480);dev.off() + + +######################################################### +#II. Compare DIURNAL CYCLES @ HPL in different years +########## +YEARs <- 2015:2020 +COLS <- c("darkblue","lightblue","darkgreen","lightgreen","orange","red") +names(COLS) <- YEARs +MONsel <- 4:9 # months to examine +#MONsel <- 4:5 # months to examine +#MONsel <- 7:8 # months to examine +SITE <- "HPL" +#ylims <- c(0,0.7) # [ppm] +#ylims <- c(0,0.7)*1000 # [ppb] +ylims <- c(0,0.5)*1000 # [ppb] +########## +MONsel.txt <- formatC(MONsel,width=2,flag="0") +# combine data from different years +dat.all<-NULL +for(i in 1:length(YEARs)){ + objname<-paste0("SimCity_CH4_allsites_hrly_",YEARs[i]) + print(paste("Reading in.....",objname)) + tmp <- getr(objname,path=obsdir) + sel <- paste0("CH4_",SITE)%in%colnames(tmp) + if(sum(sel)==length(SITE)){ + tmp2 <- tmp[,c("Time",paste0("CH4_",c(SITE,"FRU")))] + } else { + tmp2 <- data.frame(Time=tmp[,"Time"],tmp[,paste0("CH4_",c(SITE,"FRU"))]) + tmp2 <- data.frame(tmp2,NA) + colnames(tmp2)[ncol(tmp2)] <- paste0("CH4_",SITE[!sel]) + } # if(paste0("CH4_",SITE)%in%colnames(tmp)){ + dat.all<-rbind(dat.all,tmp2) + gc() +} #for(i in 1:length(YEARs)){ + +sitecolnm<-paste0("CH4_",SITE) +YYYYMMDDHH <- format(dat.all[,"Time"],format="%Y%m%d%H") +dat <- dat.all[substring(YYYYMMDDHH,5,6)%in%MONsel.txt,] +YYYYMMDDHH <- format(dat[,"Time"],format="%Y%m%d%H") + +dev.new();par(cex.axis=1.3,cex.main=1.3,cex.lab=1.3,mar=c(5,5,4,2)) +plot(0,0,type="n",xlab="",ylab="CH4 enhancement over FRU [ppb]",ylim=ylims,xlim=c(0,24),axes=FALSE) +title(main=paste("Average Diurnal Cycle:",SITE,"- FRU,\n Months:",paste(MONsel,collapse=","))) +Years <- unique(substring(YYYYMMDDHH,1,4)) +for(yy in 1:length(Years)){ + sel <- substring(YYYYMMDDHH,1,4)==Years[yy] + diur <- tapply(dat[sel,sitecolnm]-dat[sel,"CH4_FRU"],format(dat[sel,"Time"],format="%H"),mean,na.rm=T) + diur <- 1000*diur # [ppm] => [ppb] + UThrs <- as.numeric(names(diur)) + MST <- UThrs - 7; MST[MST<0] <- MST[MST<0] + 24 + # rearrange order according to MST + diur <- diur[order(MST)] + MST <- MST[order(MST)] + UThrs <- as.numeric(names(diur)) + lines(MST,diur,pch=16,type="l",col=COLS[Years[yy]],lwd=2) +} #for(yy in 1:length(Years)){ +box();axis(2) +axis(1,at=MST,labels=MST,line=0) +mtext("Hour of Day [MST]",side=1,line=2.3,cex=1.3) +#axis(1,at=MST,labels=UThrs,line=3.5);mtext("[UT]",side=1,line=4) +#abline(v=range(UTsel)-7,col="darkgray",lwd=2,lty=2) # selected hours [MST] +rect(xleft=-1,ybottom=ylims[1]-100,xright=UTsel[1]-7,ytop=ylims[2]+100,col=gray(level=0.3,alpha=0.5),border="darkgray") +rect(xleft=max(UTsel)+1-7,ybottom=ylims[1]-100,xright=24+1,ytop=ylims[2]+100,col=gray(level=0.3,alpha=0.5),border="darkgray") +legend("bottomleft",Years,bty="n",lwd=2,text.col=COLS,col=COLS,lty=1,bg="white") + + +######################################################### +#III. Compare SEASONAL CYCLES @ HPL in different years +########## +########## +dev.new();par(cex.axis=1.3,cex.main=1.3,cex.lab=1.3,mar=c(5,5,4,2)) +ylims <- c(0,0.1)*1000 # [ppb] +plot(0,0,type="n",xlab="Month",ylab="CH4 enhancement over FRU [ppb]",ylim=ylims,xlim=range(MONsel)) +title(main=paste("Average Seasonal Cycle:",SITE,"- FRU,\n UThrs:",paste(UTsel,collapse=","))) +Years <- unique(substring(YYYYMMDDHH,1,4)) +for(yy in 1:length(Years)){ + sel <- substring(YYYYMMDDHH,1,4)==Years[yy] + sel <- sel&substring(YYYYMMDDHH,9,10)%in%as.character(UTsel) + seas <- tapply(dat[sel,sitecolnm]-dat[sel,"CH4_FRU"],format(dat[sel,"Time"],format="%m"),mean,na.rm=T) + seas <- 1000*seas # [ppm] => [ppb] + Mons <- as.numeric(names(seas)) + lines(Mons,seas,pch=16,type="l",col=COLS[Years[yy]],lwd=2) +} #for(yy in 1:length(Years)){ +legend("topleft",Years,bty="n",lwd=2,text.col=COLS,col=COLS,lty=1) + + + + +######################################################### +#IV. Compare DIURNAL CYCLES @ CSP in different years +########## +YEARs <- (2015:2020)[c(2,6)] +COLS <- c("darkblue","lightblue","darkgreen","lightgreen","orange","red")[c(2,6)] +names(COLS) <- YEARs +#MONsel <- 4:9 # months to examine +MONsel <- 4:5 # months to examine +SITE <- "CSP" +ylims <- c(0,0.5)*1000 # [ppb] +########## +MONsel.txt <- formatC(MONsel,width=2,flag="0") +# combine data from different years +dat.all<-NULL +for(i in 1:length(YEARs)){ + objname<-paste0("SimCity_CH4_allsites_hrly_",YEARs[i]) + print(paste("Reading in.....",objname)) + tmp <- getr(objname,path=obsdir) + sel <- paste0("CH4_",SITE)%in%colnames(tmp) + if(sum(sel)==length(SITE)){ + tmp2 <- tmp[,c("Time",paste0("CH4_",c(SITE,"FRU")))] + } else { + tmp2 <- data.frame(Time=tmp[,"Time"],tmp[,paste0("CH4_",c(SITE,"FRU"))]) + tmp2 <- data.frame(tmp2,NA) + colnames(tmp2)[ncol(tmp2)] <- paste0("CH4_",SITE[!sel]) + } # if(paste0("CH4_",SITE)%in%colnames(tmp)){ + dat.all<-rbind(dat.all,tmp2) + gc() +} #for(i in 1:length(YEARs)){ +sitecolnm<-paste0("CH4_",SITE) +YYYYMMDDHH <- format(dat.all[,"Time"],format="%Y%m%d%H") +dat <- dat.all[substring(YYYYMMDDHH,5,6)%in%MONsel.txt,] +YYYYMMDDHH <- format(dat[,"Time"],format="%Y%m%d%H") + +dev.new();par(cex.axis=1.3,cex.main=1.3,cex.lab=1.3,mar=c(5,5,4,2)) +plot(0,0,type="n",xlab="",ylab="CH4 enhancement over FRU [ppb]",ylim=ylims,xlim=c(0,24),axes=FALSE) +title(main=paste("Average Diurnal Cycle:",SITE,"- FRU,\n Months:",paste(MONsel,collapse=","))) +Years <- unique(substring(YYYYMMDDHH,1,4)) +for(yy in 1:length(Years)){ + sel <- substring(YYYYMMDDHH,1,4)==Years[yy] + diur <- tapply(dat[sel,sitecolnm]-dat[sel,"CH4_FRU"],format(dat[sel,"Time"],format="%H"),mean,na.rm=T) + diur <- 1000*diur # [ppm] => [ppb] + UThrs <- as.numeric(names(diur)) + MST <- UThrs - 7; MST[MST<0] <- MST[MST<0] + 24 + # rearrange order according to MST + diur <- diur[order(MST)] + MST <- MST[order(MST)] + UThrs <- as.numeric(names(diur)) + lines(MST,diur,pch=16,type="l",col=COLS[Years[yy]],lwd=2) +} #for(yy in 1:length(Years)){ +box();axis(2) +axis(1,at=MST,labels=MST,line=0) +mtext("Hour of Day [MST]",side=1,line=2.3,cex=1.3) +#axis(1,at=MST,labels=UThrs,line=3.5);mtext("[UT]",side=1,line=4) +#abline(v=range(UTsel)-7,col="darkgray",lwd=2,lty=2) # selected hours [MST] +rect(xleft=-1,ybottom=ylims[1]-100,xright=UTsel[1]-7,ytop=ylims[2]+100,col=gray(level=0.3,alpha=0.5),border="darkgray") +rect(xleft=max(UTsel)+1-7,ybottom=ylims[1]-100,xright=24+1,ytop=ylims[2]+100,col=gray(level=0.3,alpha=0.5),border="darkgray") +legend("bottomleft",Years,bty="n",lwd=2,text.col=COLS,col=COLS,lty=1,bg="white") + + +######################################################### +#V. Compare SEASONAL CYCLES @ CSP in different years +########## +########## +dev.new();par(cex.axis=1.3,cex.main=1.3,cex.lab=1.3,mar=c(5,5,4,2)) +ylims <- c(0,0.1)*1000 # [ppb] +plot(0,0,type="n",xlab="Month",ylab="CH4 enhancement over FRU [ppb]",ylim=ylims,xlim=range(MONsel)) +title(main=paste("Average Seasonal Cycle:",SITE,"- FRU,\n UThrs:",paste(UTsel,collapse=","))) +Years <- unique(substring(YYYYMMDDHH,1,4)) +for(yy in 1:length(Years)){ + sel <- substring(YYYYMMDDHH,1,4)==Years[yy] + sel <- sel&substring(YYYYMMDDHH,9,10)%in%as.character(UTsel) + seas <- tapply(dat[sel,sitecolnm]-dat[sel,"CH4_FRU"],format(dat[sel,"Time"],format="%m"),mean,na.rm=T) + seas <- 1000*seas # [ppm] => [ppb] + Mons <- as.numeric(names(seas)) + lines(Mons,seas,pch=16,type="l",col=COLS[Years[yy]],lwd=2) +} #for(yy in 1:length(Years)){ +legend("topleft",Years,bty="n",lwd=2,text.col=COLS,col=COLS,lty=1) diff --git a/plot_foot_and_wells.r b/plot_foot_and_wells.r new file mode 100644 index 0000000..e04c065 --- /dev/null +++ b/plot_foot_and_wells.r @@ -0,0 +1,155 @@ +# Plots footprint of receptors and adds in well locations +# V2(210418): instead of calculating averaged footprints on the fly, reads in output from 'average_foot_monthly.r' +# V3(210420): add AVERAGE footprint map with well locations and site locations (meant for paper) +# March 21st, 2021 by John C. Lin (John.Lin@utah.edu) + +######################### +#SITE<-"HPL" +#SITE<-"ROO" +SITE<-"CSP" + +# Uintah Basin domain +#XLIMS <- c(-111.05,-108.9);YLIMS <- c(39.4,41.0) #set lat/lon limits for limits of map; set to NULL if dynamically determine +XLIMS <- c(-111.55,-108.9);YLIMS <- c(39.4,41.0) #set lat/lon limits for limits of map; set to NULL if dynamically determine +#XLIMS <- c(-111.15,-108.7);YLIMS <- c(39.4,41.0) #set lat/lon limits for limits of map; set to NULL if dynamically determine +#XLIMS <- c(-111.10,-108.8);YLIMS <- c(39.4,41.0) #set lat/lon limits for limits of map; set to NULL if dynamically determine +#range of footprint plot (log10 limits) + +HRs <- 20:23 #[UTC] only analyze afternoon hours, following Foster papers + +welldatname <- "State_of_Utah_oil_gas_DATA/well_data_gas.oil.production_monthly.rds" # generated by 'State_of_Utah_oil_gas_DATA/merge_well_data_production.r' + +# meteorology driving STILT +mettype <- "HRRR" +#mettype <- "NAM12" +#mettype <- "WRF27" + +#footmonavenm <- paste0(SITE,'_foot_monave.rds') # output from 'average_foot_monthly.r' +footmonavenm <- paste0(SITE,'_foot_monave_',mettype,'.rds') # output from 'average_foot_monthly.r' + +#YYYYMMs <- c("201604","201605","202004","202005","202008") # selected Years/Months +MONs <- 4:9; YEARs <- 2015:2020 +if(SITE=="CSP") YEARs <- c(2016,2019,2020) +YYYYMMs <- paste0(rep(YEARs,each=length(MONs)),rep(formatC(MONs,width=2,flag="0"),length(YEARs))) +YYYYMMs <- YYYYMMs[!YYYYMMs%in%c("201504","201505","201506")] +#log10TF <- TRUE; ZLIMS <- c(-7,-1) +#log10TF <- FALSE; ZLIMS <- c(0.00001,0.0001) +log10TF <- FALSE; ZLIMS <- c(0.000005,0.0001) +######################### + +require(ncdf4);require(fields) +require(maps);require(RgoogleMaps) +if(SITE=="HPL"){LAT <- 40.1434; LON <- -109.4680; zagl <- 4.06} +if(SITE=="ROO"){LAT <- 40.2941; LON <- -110.0090; zagl <- 4.06} +if(SITE=="CSP"){LAT <- 40.0509; LON <- -110.0194; zagl <- 4.06} + +SITES <- data.frame(matrix(NA,nrow=4,ncol=3,dimnames=list(c("HPL","ROO","CSP","FRU"),c("LAT","LON","zagl")))) +for(i in 1:nrow(SITES)){ + if(rownames(SITES)[i]=="HPL"){SITES$LAT[i] <- 40.1434; SITES$LON[i] <- -109.4680; SITES$zagl[i] <- 4.06} + if(rownames(SITES)[i]=="ROO"){SITES$LAT[i] <- 40.2941; SITES$LON[i] <- -110.0090; SITES$zagl[i] <- 4.06} + if(rownames(SITES)[i]=="CSP"){SITES$LAT[i] <- 40.0509; SITES$LON[i] <- -110.0194; SITES$zagl[i] <- 4.06} + if(rownames(SITES)[i]=="FRU"){SITES$LAT[i] <- 40.2087; SITES$LON[i] <- -110.8403; SITES$zagl[i] <- 4.06} +} #for(i in nrow(SITES)){ +# not show ROO in the overall plot due to contamination from nearby well... +SITES <- SITES[rownames(SITES)!="ROO",] + +DAT <- readRDS(welldatname) + +foot.ave.ALL <- readRDS(footmonavenm) + +#-------------------------------------------------# +# average of ALL of the selected Year/Months # +foot.AVE <- apply(foot.ave.ALL[,,YYYYMMs],c(1,2),mean,na.rm=T) + + +#-------------------------------------------------# +# Plot month-by-month # +#for(i in 0:length(YYYYMMs)){ +for(i in 0){ + if(i==0){ + foot <- foot.AVE + wdat <- DAT + # Well.Type.y and Well.Status.y is the type/status of the well at the time when data downloaded (2021)--March 2021 + # relevant files are "CH4_inversion_Uintah/State_of_Utah_oil_gas_DATA/well_data__all.csv" + oil.producing <- wdat[wdat$Well.Type.y=="Oil Well"&wdat$Well.Status.y=="Producing",] + gas.producing <- wdat[wdat$Well.Type.y=="Gas Well"&wdat$Well.Status.y=="Producing",] + YYYYMM <- paste("\nfootprint averaged over selected years/mons") + saveRDS(oil.producing,file="oil.producing.rds") + saveRDS(gas.producing,file="gas.producing.rds") + saveRDS(foot,file=paste0(SITE,"_footAVE.rds")) + system(paste0("tar cvf ~/public_html/Exchange/rds.tar oil.producing.rds gas.producing.rds ",SITE,"_footAVE.rds")) + } else { + YYYYMM <- YYYYMMs[i] + foot <- foot.ave.ALL[,,YYYYMM] + # prepare well data + sel <- substring(DAT$YYYYMMDD,1,6) == YYYYMM + wdat <- DAT[sel,] + oil.producing <- wdat[wdat$Well.Type.x=="Oil Well"&wdat$Well.Status.x=="Producing",] + gas.producing <- wdat[wdat$Well.Type.x=="Gas Well"&wdat$Well.Status.x=="Producing",] + } # } else { + + flon <- as.numeric(rownames(foot)); flat <- as.numeric(colnames(foot)) + foot[foot==0] <- NA + foot.log10 <- log10(foot) + + xmain<-paste0(SITE,": ",YYYYMM) + xsub<-paste0("UThrs: ",paste(HRs,collapse=","),"; mettype=",mettype) + #plot on Google Maps + alpha<-0.70 #transparency + bb <- qbbox(lat=YLIMS, lon=XLIMS, TYPE="quantile") + par(pty="s");MyMap <- GetMap.bbox(bb$lonR, bb$latR,urlBase = "http://mt1.google.com/vt/lyrs=m", tileDir= "./mapTiles/Google/") + PlotOnStaticMap(MyMap,lat=LAT,lon=LON,col="red",pch=17,FUN=points,TrueProj = TRUE,mar=c(2,2,3,5)) + PlotOnStaticMap(MyMap,lat=gas.producing$Latitude,lon=gas.producing$Longitude, + col="blue",pch=16,cex=0.4,FUN=points,TrueProj = TRUE,add=T) + PlotOnStaticMap(MyMap,lat=oil.producing$Latitude,lon=oil.producing$Longitude, + col="orange",pch=16,cex=0.4,FUN=points,TrueProj = TRUE,add=T) + title(main=xmain,cex.main=1.5) + mtext(xsub,side=1,cex=1.1) + + lats <- flat; lons <- flon + lons.mat<-matrix(rep(lons,times=length(lats)),ncol=length(lats)) + lats.mat<-matrix(rep(lats,times=length(lons)),byrow=T,nrow=length(lons)) + image.coords <- LatLon2XY.centered(MyMap, lat=as.vector(lats.mat), lon=as.vector(lons.mat)) + xnew.mat<-matrix(image.coords$newX,ncol=ncol(lons.mat),nrow=nrow(lons.mat)) + xnew<-xnew.mat[,1] + ynew.mat<-matrix(image.coords$newY,ncol=ncol(lats.mat),nrow=nrow(lats.mat)) + ynew<-ynew.mat[1,] + + zmat <- foot + #legend.lab <- "footprint" + legend.lab <- expression(paste("footprint [",mu,"mole ",m^-2," ",s^-1,"]",sep="")) + if(log10TF){zmat <- foot.log10;legend.lab <- "log10(footprint)"} + #figure out color scheme & zlims + COLS<-c("white","violet","lightblue","blue","yellow","orange","red") #JCL: just a subset of colors--not to conflict with trajectory color + # ZLIMS<-range(zmat,na.rm=T) + #colsc<-designer.colors(n=64,alpha=alpha,col=COLS); colsc.legend<-designer.colors(n=64,alpha=1.0,col=COLS) + #colsc <- tim.colors(n=64,alpha=alpha); colsc.legend <- tim.colors(n=64,alpha=1.0) + colsc <- gray(seq(1,0,-0.02),alpha=alpha); colsc.legend <- gray(seq(1,0,-0.02),alpha=1.0) + zmat2<-zmat + zmat2[zmat2>ZLIMS[2]]<-ZLIMS[2] #make sure everything gets plotted on colorscale + #zmat2[zmat20) + sel <- welldat$Well.Type.x == "Oil Well" + PlotOnStaticMap(MyMap,lat=welldat$Latitude[sel],lon=welldat$Longitude[sel],col="darkgray",FUN=points,TrueProj = TRUE,add=T,cex=0.4,pch=16) + sel <- welldat$Well.Type.x == "Gas Well" + PlotOnStaticMap(MyMap,lat=welldat$Latitude[sel],lon=welldat$Longitude[sel],col="black",FUN=points,TrueProj = TRUE,add=T,cex=0.4,pch=16) + } # if(plotwellsTF) + + lats <- FLAT; lons <- FLON + lons.mat<-matrix(rep(lons,times=length(lats)),ncol=length(lats)) + lats.mat<-matrix(rep(lats,times=length(lons)),byrow=T,nrow=length(lons)) + image.coords <- LatLon2XY.centered(MyMap, lat=as.vector(lats.mat), lon=as.vector(lons.mat)) + xnew.mat<-matrix(image.coords$newX,ncol=ncol(lons.mat),nrow=nrow(lons.mat)) + xnew<-xnew.mat[,1] + ynew.mat<-matrix(image.coords$newY,ncol=ncol(lats.mat),nrow=nrow(lats.mat)) + ynew<-ynew.mat[1,] + + dat <- dat.all[,,YYYYMM] + dat <- dat*AREA.mat # remove [m^2] from denominator + TOT <- sum(dat,na.rm=T) + dat[dat==0] <- NA + zmat <- dat + #figure out color scheme & zlims + #COLS<-c("white","violet","lightblue","blue","yellow","orange","red") #JCL: just a subset of colors--not to conflict with trajectory color + #colsc<-designer.colors(n=64,alpha=alpha,col=COLS);colsc.legend<-designer.colors(n=64,alpha=1.0,col=COLS) + #colsc <- plasma(n=20,alpha=alpha);colsc.legend <- plasma(n=20,alpha=1.0) + colsc<-tim.colors(n=60,alpha=alpha);colsc.legend<-tim.colors(n=60,alpha=1.0) + + zmat2 <- zmat + ZLIMS <- range(zmat2,na.rm=T) + if(VAR=='Mcf.Gas')ZLIMS<-c(0.05,0.5) + if(VAR=='Bbls.Oil')ZLIMS<-c(0.05,0.5) + zmat2[zmat2ZLIMS[2]]<-ZLIMS[2] #make sure everything gets plotted on colorscale + image(x=xnew,y=ynew,z=zmat2,add=TRUE,col=colsc,zlim=ZLIMS) + image.plot(x=xnew,y=ynew,z=zmat2,legend.only=TRUE,add=TRUE,col=colsc.legend,legend.lab=VAR,zlim=ZLIMS,horizontal=TRUE) + PlotOnStaticMap(MyMap,lat=LAT,lon=LON,col="black",pch=17,FUN=points,TrueProj = TRUE,add=T,cex=1.5) + PlotOnStaticMap(MyMap,lat=LAT+0.05,lon=LON,col="black",FUN=text,TrueProj = TRUE,add=TRUE,labels=SITE) + xmain <- paste0(SITE," ",VAR," ",YYYYMM,"\nUThrs: ",paste(UThrs,collapse=","),"; total=",signif(TOT,3),"; mettype=",mettype) + title(main=xmain,cex.main=1.3) + #mtext(xsub,side=1,cex=1.1) + + if(plotwellsTF)legend("topleft",c("Gas Well","Oil Well"),pch=16,col=c("black","darkgray")) + + figfilenm <- paste0(SITE,"_",VAR,"_foot_wellinfo_convolved_",YYYYMM,"_",mettype,".png") + dev.copy(png,filename=figfilenm);dev.off();print(paste(figfilenm,"generated")) + gc() +} # for(i in 1:length(YYYYMMs)){ + + + diff --git a/plot_foot_wellinfo_tseries_trend.r b/plot_foot_wellinfo_tseries_trend.r new file mode 100644 index 0000000..3468b86 --- /dev/null +++ b/plot_foot_wellinfo_tseries_trend.r @@ -0,0 +1,263 @@ +# Plot time series of footprint-convolved well # or footprint-convolved CH4 production and also determine trends +# +# V2(210417): filter out days based on within 45o of observed wind direction (filterUV.TF) +# April 5th, 2021 by John C. Lin (John.Lin@utah.edu) + +################ +SITE<-"HPL" +#SITE<-"ROO" +#SITE<-"CSP" +# subset of months +MONsub<-4:9 +# subset of hours [UTC]; if NULL, then all hours +UThrs <- c(20,21,22,23) # [UTC] only analyze afternoon hours, following Foster papers +#datname <- paste0(SITE,'_CH4_foot.wellinfo_ALL.rds') # output from 'merge_CH4_well_info.r' +datname <- paste0(SITE,'_CH4_UVwind_foot.wellinfo_ALL.rds') # output from 'merge_CH4_UVwind_wellinfo.r' +filterUV.TF <- TRUE # filter times based on U/V (filter out times when HRRR is off--i.e., large transport errors) ? +################ + +require(ggplot2); require(reshape2); require(geosphere) + +if(SITE=="HPL"){LAT <- 40.1434; LON <- -109.4680; zagl <- 4.06; WINDsite <- "UBHSP"; WINDsite.2015 <- "A1633"} # UBHSP site not available in 2015, so use A1633 (RedWash) site instead in 2015 +if(SITE=="ROO"){LAT <- 40.2941; LON <- -110.0090; zagl <- 4.06; WINDsite <- "QRS"} +if(SITE=="CSP"){LAT <- 40.0509; LON <- -110.0194; zagl <- 4.06; WINDsite <- "UBCSP"} + +DAT.all <- readRDS(datname) + +# filter subset of dataset based on Months & UThrs +YYYYMMDDHH <- format(DAT.all$Time,'%Y%m%d%H') +DAT.all <- data.frame(YYYYMMDDHH, DAT.all) +rownames(DAT.all) <- YYYYMMDDHH +sel <- substring(DAT.all$YYYYMMDDHH, 5, 6)%in%formatC(MONsub,width=2,flag='0') +sel <- sel&(substring(DAT.all$YYYYMMDDHH, 9, 10)%in%formatC(UThrs,width=2,flag='0')) +DAT.all <- DAT.all[sel,] +# filter out times with missing obs at SITE +sel <- !is.na(DAT.all[,paste0("CH4_",SITE)]) +DAT.all <- DAT.all[sel,] + +if(filterUV.TF){ + dat <- DAT.all + Usim <- dat[,paste0("Usim.",WINDsite)]; Vsim <- dat[,paste0("Vsim.",WINDsite)] + Uobs <- dat[,paste0("Uobs.",WINDsite)]; Vobs <- dat[,paste0("Vobs.",WINDsite)] + isNA <- is.na(Usim)|is.na(Vsim)|is.na(Uobs)|is.na(Vobs) + # a) U/V filter: both simulated U & V have to have the same sign as the observed (i.e., wind vector has to be same QUADRANT as observed) + # SEL <- sign(Usim)!=sign(Uobs) + # SEL <- sel|(sign(Vsim)!=sign(Vobs)) + + # b) U/V filter: filter out days based on +/-45o of observed (better than quadrant method above, since not throw way data close to quadrant boundaries + WSPD.sim <- sqrt(Usim^2 + Vsim^2); WSPD.obs <- sqrt(Uobs^2 + Vobs^2) + print("Observed windspeeds (quantile):") + print(signif(quantile(WSPD.obs,c(0,0.05,0.1,0.25,0.5,0.75,1.0),na.rm=T),3)) + # see: http://tornado.sfsu.edu/geosciences/classes/m430/Wind/WindDirection.html + WDIR.sim <- (180/pi)*atan2(y=-1*Vsim,x=-1*Usim) # wind direction (-180o to +180o); direction wind is blowing FROM + WDIR.obs <- (180/pi)*atan2(y=-1*Vobs,x=-1*Uobs) # wind direction (-180o to +180o); direction wind is blowing FROM + dWDIR <- WDIR.sim - WDIR.obs + sel <- abs(dWDIR) > 180 + dWDIR[sel&!isNA] <- sign(dWDIR[sel&!isNA])*(360 - abs(dWDIR[sel&!isNA])) + SEL <- abs(dWDIR) > 45 # simulated wind direction has to be within 45o of observed + SEL[WSPD.obs < 1.0] <- FALSE # not remove data if windspeed is too low (wind anemometer not reliable) + + print(paste("Out of",sum(!isNA),"non-NA U/V times,",sum(SEL[!isNA]),"failed to pass U/V filter: ",signif(100*sum(SEL[!isNA])/sum(!isNA),3),"%")) + # assign NAs to wellinfo when failed to pass U/V filter + print(paste("Before filtering, Mcf.Gas has",sum(!is.na(dat$Mcf.Gas)),"non-NA values")) + dat[SEL&!isNA,COLNMS] <- NA + print(paste(" After filtering, Mcf.Gas has",sum(!is.na(dat$Mcf.Gas)),"non-NA values")) + DAT.all <- dat +} # if(filterUV.TF){ + +# multiply by AREA of footprint grid to remove [1/m^2] factor and to make the slope values larger +# footprint grid settings from 'run_stilt_hysplit-stiltV1.r' +xmn <- -112.5; xmx <- -108 +ymn <- 38.5; ymx <- 41.5 +xres <- 0.01; yres <- xres +flat <- seq(ymn+yres/2,ymx-yres/2,yres) +flon <- seq(xmn+xres/2,xmx-xres/2,xres) +# lat/lons refer to footprint grid CENTERS, so move them to lower left (southwest) corner to facilitate gridding +flat.ll <- flat - yres/2; flon.ll <- flon - xres/2 +# determine footprint gridcell areas +flon.vec <- rep(flon,length(flat)) +flat.vec <- rep(flat,each=length(flon)) +dx <- distHaversine(p1=cbind(flon.vec,flat.vec),p2=cbind(flon.vec+xres,flat.vec)) +dy <- mean(distHaversine(p1=cbind(flon.vec,flat.vec),p2=cbind(flon.vec,flat.vec+yres))) # gridcell dimension in y-direction is almost a constant +AREA <- dx*dy # gridcell area [m^2] +AREA.mat <- matrix(AREA,nrow=length(flon)) #; image.plot(AREA.mat) +#!!!!! average footprint gridcell area is going to be used to remove [1/m^2] factor !!!!!# +AREA.ave <- mean(AREA.mat) + +COLNMS <- c('Mcf.Gas','Bbls.Oil','Gas.Well','Oil.Well','Gas.Well_Producing','Oil.Well_Producing') +DAT.all[,COLNMS] <- DAT.all[,COLNMS]*AREA.ave + + +DAT0 <- DAT.all[,c("Time",COLNMS)] + +#----------------------------------------- All Selected Hours ---------------------------------------# +DAT <- DAT0 +wellvars <- COLNMS +Time <- DAT$Time +data.table <- NULL +for(i in 1:length(wellvars)){ + Wellinfo <- wellvars[i] + well <- DAT[,Wellinfo] + xlm <- lm(well ~ Time) + m <- coef(xlm)[2] + p.value <- summary(xlm)$coefficients[,"Pr(>|t|)"]["Time"] + data.table <- rbind(data.table,data.frame(Wellinfo,m,p.value)) +} # for(i in 1:length(wellvars)){ +rownames(data.table) <- NULL +# slope ~ [1/s] due to Time being in seconds, so change to [1/year] +data.table$m <- data.table$m*(60*60*24*365) + +dev.new(); theme_set(theme_bw()) +dat <- melt(DAT, id.vars='Time') +colnames(dat)[colnames(dat)=="variable"] <- "Wellinfo" +g <- ggplot(dat, aes(x=Time,y=value)) + geom_point(,size=1.5) + + stat_smooth(method = "lm") + geom_abline() + facet_grid(. ~ Wellinfo) +g <- g + labs(title=paste0('SITE =',SITE,'; Time series of sum(foot*wellinfo)', + '\n UThrs=',paste(sort(unique(substring(rownames(DAT),9,10))),collapse=","),'; MONs: ',paste(MONsub,collapse=",")), + caption=paste0('filterUV.TF=',filterUV.TF)) +g <- g + theme(strip.text.x = element_text(size = 14, colour = "black"),plot.caption=element_text(size=12), + axis.title.y=element_text(size=16), axis.title.x=element_text(size=16), + axis.text.x=element_text(size=12), axis.text.y=element_text(size=14)) +g <- g + geom_text(data = data.table, hjust = 0, col='blue', + aes(x = dat$Time[1], y = Inf, label = paste0('\nSlope=',signif(m,3),'; p-val=',signif(p.value,2)))) +g + facet_wrap(.~Wellinfo, ncol=2, scales='free') + + +#----------------------------------------- Average Selected Hours within a Day ---------------------------------------# +YYYYMMDD <- substring(DAT.all$YYYYMMDDHH,1,8) +dat <- DAT.all[,c(COLNMS)] +dum <- NULL +for(cc in 1:ncol(dat)){ + tmp <- tapply(dat[,cc],YYYYMMDD,mean,na.rm=T) + dum <- cbind(dum,tmp) +} # for(cc in 1:ncol(dat)){ +colnames(dum) <- c(COLNMS) +dum[is.nan(dum)] <- NA + +Time <- as.POSIXct(rownames(dum),format="%Y%m%d",tz="GMT") +DAT <- data.frame(Time,dum) +wellvars <- COLNMS +data.table <- NULL +for(i in 1:length(wellvars)){ + Wellinfo <- wellvars[i] + well <- DAT[,Wellinfo] + xlm <- lm(well ~ Time) + m <- coef(xlm)[2] + p.value <- summary(xlm)$coefficients[,"Pr(>|t|)"]["Time"] + data.table <- rbind(data.table,data.frame(Wellinfo,m,p.value)) +} # for(i in 1:length(wellvars)){ +rownames(data.table) <- NULL +# slope ~ [1/s] due to Time being in seconds, so change to [1/year] +data.table$m <- data.table$m*(60*60*24*365) + +dev.new(); theme_set(theme_bw()) +dat <- melt(DAT, id.vars='Time') +colnames(dat)[colnames(dat)=="variable"] <- "Wellinfo" +g <- ggplot(dat, aes(x=Time,y=value)) + geom_point(,size=1.5) + + stat_smooth(method = "lm") + geom_abline() + facet_grid(. ~ Wellinfo) +g <- g + labs(title=paste0('SITE =',SITE,'; Time series of sum(foot*wellinfo)', + '\n mean(UThrs=',paste(UThrs,collapse=","),')','; MONs: ',paste(MONsub,collapse=",")), + caption=paste0('filterUV.TF=',filterUV.TF)) +g <- g + theme(strip.text.x = element_text(size = 14, colour = "black"),plot.caption=element_text(size=12), + axis.title.y=element_text(size=16), axis.title.x=element_text(size=16), + axis.text.x=element_text(size=12), axis.text.y=element_text(size=14)) +g <- g + geom_text(data = data.table, hjust = 0, col='blue', + aes(x = dat$Time[1], y = Inf, label = paste0('\nSlope=',signif(m,3),'; p-val=',signif(p.value,2)))) +g + facet_wrap(.~Wellinfo, ncol=2, scales='free') +figfilenm <- paste0(SITE,"_foot_wellinfo_tseries_daily.png") +ggsave(figfilenm);print(paste(figfilenm,"generated")) + +#----------------------------------------- Average over an entire MONTH ---------------------------------------# +YYYYMM <- substring(DAT.all$YYYYMMDDHH,1,6) +dat <- DAT.all[,c(COLNMS)] +dum <- NULL +for(cc in 1:ncol(dat)){ + tmp <- tapply(dat[,cc],YYYYMM,mean,na.rm=T) + dum <- cbind(dum,tmp) +} # for(cc in 1:ncol(dat)){ +colnames(dum) <- c(COLNMS) +dum[is.nan(dum)] <- NA + +YYYYMMDD <- paste0(rownames(dum),"15") # create this, since as.POSIXct doesn't work with just %Y%m, need DAYS too +Time <- as.POSIXct(YYYYMMDD,format="%Y%m%d",tz="GMT") +DAT <- data.frame(Time,dum) +wellvars <- COLNMS +data.table <- NULL +for(i in 1:length(wellvars)){ + Wellinfo <- wellvars[i] + well <- DAT[,Wellinfo] + xlm <- lm(well ~ Time) + m <- coef(xlm)[2] + p.value <- summary(xlm)$coefficients[,"Pr(>|t|)"]["Time"] + data.table <- rbind(data.table,data.frame(Wellinfo,m,p.value)) +} # for(i in 1:length(wellvars)){ +rownames(data.table) <- NULL +# slope ~ [1/s] due to Time being in seconds, so change to [1/year] +data.table$m <- data.table$m*(60*60*24*365) + +dev.new(); theme_set(theme_bw()) +dat <- melt(DAT, id.vars='Time') +colnames(dat)[colnames(dat)=="variable"] <- "Wellinfo" +g <- ggplot(dat, aes(x=Time,y=value)) + geom_point(,size=1.5) + + stat_smooth(method = "lm") + geom_abline() + facet_grid(. ~ Wellinfo) +g <- g + labs(title=paste0('SITE =',SITE,'; Time series of sum(foot*wellinfo)', + '\n MONTHLY mean; MONs: ',paste(MONsub,collapse=",")), + caption=paste0('filterUV.TF=',filterUV.TF)) +g <- g + theme(strip.text.x = element_text(size = 14, colour = "black"),plot.caption=element_text(size=12), + axis.title.y=element_text(size=16), axis.title.x=element_text(size=16), + axis.text.x=element_text(size=12), axis.text.y=element_text(size=14)) +g <- g + geom_text(data = data.table, hjust = 0, col='blue', + aes(x = dat$Time[1], y = Inf, label = paste0('\nSlope=',signif(m,3),'; p-val=',signif(p.value,2)))) +g + facet_wrap(.~Wellinfo, ncol=2, scales='free') +figfilenm <- paste0(SITE,"_foot_wellinfo_tseries_monthly.png") +ggsave(figfilenm);print(paste(figfilenm,"generated")) + +#----------------------------------------- Average over an entire YEAR ----------------------------------------# +YYYY <- substring(DAT.all$YYYYMMDDHH,1,4) +dat <- DAT.all[,c(COLNMS)] +dum <- NULL +for(cc in 1:ncol(dat)){ + tmp <- tapply(dat[,cc],YYYY,mean,na.rm=T) + dum <- cbind(dum,tmp) +} # for(cc in 1:ncol(dat)){ +colnames(dum) <- c(COLNMS) +dum[is.nan(dum)] <- NA + +YYYYMMDD <- paste0(rownames(dum),"0615") # create this, since as.POSIXct doesn't work with just %Y%m, need DAYS too +Time <- as.POSIXct(YYYYMMDD,format="%Y%m%d",tz="GMT") +DAT <- data.frame(Time,dum) +wellvars <- COLNMS +data.table <- NULL +for(i in 1:length(wellvars)){ + Wellinfo <- wellvars[i] + well <- DAT[,Wellinfo] + xlm <- lm(well ~ Time) + m <- coef(xlm)[2] + p.value <- summary(xlm)$coefficients[,"Pr(>|t|)"]["Time"] + data.table <- rbind(data.table,data.frame(Wellinfo,m,p.value)) +} # for(i in 1:length(wellvars)){ +rownames(data.table) <- NULL +# slope ~ [1/s] due to Time being in seconds, so change to [1/year] +data.table$m <- data.table$m*(60*60*24*365) + +dev.new(); theme_set(theme_bw()) +dat <- melt(DAT, id.vars='Time') +colnames(dat)[colnames(dat)=="variable"] <- "Wellinfo" +g <- ggplot(dat, aes(x=Time,y=value)) + geom_point(,size=1.5) + + stat_smooth(method = "lm") + geom_abline() + facet_grid(. ~ Wellinfo) +g <- g + labs(title=paste0('SITE =',SITE,'; Time series of sum(foot*wellinfo)', + '\n ANNUAL mean; MONs: ',paste(MONsub,collapse=",")), + caption=paste0('filterUV.TF=',filterUV.TF)) +g <- g + theme(strip.text.x = element_text(size = 14, colour = "black"),plot.caption=element_text(size=12), + axis.title.y=element_text(size=16), axis.title.x=element_text(size=16), + axis.text.x=element_text(size=12), axis.text.y=element_text(size=14)) +g <- g + geom_text(data = data.table, hjust = 0, col='blue', + aes(x = dat$Time[1], y = Inf, label = paste0('\nSlope=',signif(m,3),'; p-val=',signif(p.value,2)))) +g + facet_wrap(.~Wellinfo, ncol=2, scales='free') +figfilenm <- paste0(SITE,"_foot_wellinfo_tseries_yrly.png") +ggsave(figfilenm);print(paste(figfilenm,"generated")) + + + +#----------------------------------------- Separate out Seasons ---------------------------------------# + diff --git a/plot_gridded_wellinfoV1.r b/plot_gridded_wellinfoV1.r new file mode 100644 index 0000000..da18eee --- /dev/null +++ b/plot_gridded_wellinfoV1.r @@ -0,0 +1,115 @@ +# Plot the gridded well info onto Google Map +# March 28th, 2021 by John C. Lin (John.Lin@utah.edu) + +################ +# output generated by 'grid_well_info_foot.r' +VARS <- c("Mcf.Gas","Bbls.Oil","Gas.Well","Oil.Well","Gas.Well_Producing","Oil.Well_Producing")[c(2,5,6)] + +# year/month(s) that want to plot +YYYYMMs <- NULL # if NULL, then ALL of the YYYMMs +#YYYYMMs <- c("201604","202004") +#YYYYMMs <- c("201507","201707","201807","202007")[3] + +# Uintah Basin domain +XLIMS <- c(-111.05,-108.9);YLIMS <- c(39.4,41.0) #set lat/lon limits for limits of map; set to NULL if dynamically determine +################ + +require(ncdf4);require(fields) +require(maps);require(RgoogleMaps) +require(viridis);require(geosphere) + + +SITEs <- c("CSP","HPL","ROO") +LATs <- NULL; LONs <- NULL +for(SITE in SITEs){ + if(SITE=="HPL"){LAT <- 40.1434; LON <- -109.4680; zagl <- 4.06} + if(SITE=="ROO"){LAT <- 40.2941; LON <- -110.0090; zagl <- 4.06} + if(SITE=="CSP"){LAT <- 40.0509; LON <- -110.0194; zagl <- 4.06} + LATs <- c(LATs,LAT); LONs <- c(LONs,LON) +} #for(ss in 1:length(SITEs)){ + + +for(vv in 1:length(VARS)){ + var <- VARS[vv] + resultname <- paste0(var,"_gridded.rds") + dat.all <- readRDS(resultname) + + # determine footprint gridcell areas + FLON <- as.numeric(rownames(dat.all)); FLAT <- as.numeric(colnames(dat.all)) + xres <- unique(signif(diff(FLON),3)); yres <- unique(signif(diff(FLAT),3)) + flon.vec <- rep(FLON,ncol(dat.all)) + flat.vec <- rep(FLAT,each=nrow(dat.all)) + dx <- distHaversine(p1=cbind(flon.vec,flat.vec),p2=cbind(flon.vec+xres,flat.vec)) + dy <- mean(distHaversine(p1=cbind(flon.vec,flat.vec),p2=cbind(flon.vec,flat.vec+yres))) # gridcell dimension in y-direction is almost a constant + AREA <- dx*dy # gridcell area [m^2] + AREA.mat <- matrix(AREA,nrow=length(FLON)) #; image.plot(AREA.mat) + + YYYYMMsels <- YYYYMMs + if(is.null(YYYYMMsels))YYYYMMsels <- dimnames(dat.all)[[3]] +for(i in 1:length(YYYYMMsels)){ + YYYYMM <- YYYYMMsels[i] + dat <- dat.all[,,YYYYMM] + # remove [m^2] from denominator, resulting in units of [var/day] + dat <- dat*AREA.mat + + xmain<-paste0(var,": ",YYYYMM,"\ntotal=",signif(sum(dat,na.rm=T),3)," [per day]") + xsub<-"" + #plot on Google Maps + alpha<-0.50 #transparency + bb <- qbbox(lat=YLIMS, lon=XLIMS, TYPE="quantile") + par(pty="s");MyMap <- GetMap.bbox(bb$lonR, bb$latR,urlBase = "http://mt1.google.com/vt/lyrs=m", tileDir= "./mapTiles/Google/") + PlotOnStaticMap(MyMap,lat=LATs,lon=LONs,col="black",pch=17,FUN=points,TrueProj = TRUE,mar=c(2,2,3,2)) + #PlotOnStaticMap(MyMap,lat=LATs+0.03,lon=LONs,col="black",FUN=text,TrueProj = TRUE,add=TRUE,labels=SITEs) + title(main=xmain,cex.main=1.5) + mtext(xsub,side=1,cex=1.1) + + lats <- as.numeric(colnames(dat)); lons <- as.numeric(rownames(dat)) + lons.mat<-matrix(rep(lons,times=length(lats)),ncol=length(lats)) + lats.mat<-matrix(rep(lats,times=length(lons)),byrow=T,nrow=length(lons)) + image.coords <- LatLon2XY.centered(MyMap, lat=as.vector(lats.mat), lon=as.vector(lons.mat)) + xnew.mat<-matrix(image.coords$newX,ncol=ncol(lons.mat),nrow=nrow(lons.mat)) + xnew<-xnew.mat[,1] + ynew.mat<-matrix(image.coords$newY,ncol=ncol(lats.mat),nrow=nrow(lats.mat)) + ynew<-ynew.mat[1,] + + # set 0s to NA to increase transparency of plot + zmat <- dat + zmat[zmat==0] <- NA + ZLIMS<-range(zmat,na.rm=TRUE) + if(var=='Mcf.Gas')ZLIMS <- c(100,1E4) + if(var=='Bbls.Oil')ZLIMS <- c(0,150) + if(var=='Gas.Well'|var=='Gas.Well_Producing')ZLIMS <- c(0,1.5) # number of gas wells in each gridcell + if(var=='Oil.Well'|var=='Oil.Well_Producing')ZLIMS <- c(0,0.5) # number of oil wells in each gridcell + #figure out color scheme & zlims + # COLS<-c("white","violet","lightblue","blue","yellow","orange","red") #JCL: just a subset of colors--not to conflict with trajectory color + # colsc<-designer.colors(n=64,alpha=alpha,col=COLS) + # colsc.legend<-designer.colors(n=64,alpha=1.0,col=COLS) + #colsc <- tim.colors(n=64,alpha=alpha) + #colsc.legend <- tim.colors(n=64,alpha=1.0) + colsc <- plasma(n=20,alpha=alpha) + colsc.legend <- plasma(n=20,alpha=1.0) + + #make sure everything gets plotted on colorscale + # zmat[zmatZLIMS[2]] <- ZLIMS[2] + image(x=xnew,y=ynew,z=zmat,add=TRUE,col=colsc,zlim=ZLIMS) + image.plot(x=xnew,y=ynew,z=zmat,legend.only=TRUE,add=TRUE,col=colsc.legend,legend.lab=paste(var,"[per day]"),zlim=ZLIMS,horizontal=TRUE) + PlotOnStaticMap(MyMap,lat=LAT,lon=LON,col="black",pch=16,FUN=points,TrueProj = TRUE,add=T) + + #plot site locations/names + PlotOnStaticMap(MyMap,lat=LATs,lon=LONs,col="black",pch=17,FUN=points,TrueProj = TRUE,mar=c(2,2,3,2),add=TRUE) + PlotOnStaticMap(MyMap,lat=LATs+0.03,lon=LONs,col="black",FUN=text,TrueProj = TRUE,add=TRUE,labels=SITEs) + + #plot outliers for OIL production + #PlotOnStaticMap(MyMap,lat=c(40.255,40.255),lon=c(-110.315,-110.035),col="red",FUN=points,pch=16,TrueProj = TRUE,add=TRUE,cex=1.5) + + # figfilenm <- paste0(var,"_",YYYYMM,".jpg"); dev.copy(jpeg,filename=figfilenm) + figfilenm <- paste0(var,"_",YYYYMM,".png"); dev.copy(png,filename=figfilenm) + dev.off();print(paste(figfilenm,"generated")) + gc() +} # for(i in 1:length(YYYYMMs)){ + # create movie to facilitate visualization of changes in well activity + system(paste0("convert -dispose Background -loop 0 -delay 40 ",VARS[vv],"_??????.png ~/public_html/Exchange/",VARS[vv],".gif")) + gc() +} # for(vv in 1:length(VARS)){ + diff --git a/plot_tseries_and_COVID.r b/plot_tseries_and_COVID.r new file mode 100644 index 0000000..376ed22 --- /dev/null +++ b/plot_tseries_and_COVID.r @@ -0,0 +1,245 @@ +# Plots time series & diurnal cycle of CH4 in Uintah Basin and contrast 2020 (COVID-shutdown period) against previous years +# Based on "Fch4_simple.r" +# September 2nd, 2020 by John C. Lin (John.Lin@utah.edu) + +require(ncdf4); require(fields); require(geosphere) +#################### +YEARs<-2015:2020 +MONs<-1:12 +HRs<-20:23 #[UTC] only analyze afternoon hours, following Foster papers + +SITE<-"HPL" #site to focus on for calculating F_CH4 +#SITE<-"CSP" #site to focus on for calculating F_CH4 +obsdir<-"/uufs/chpc.utah.edu/common/home/lin-group4/jcl/SimCity/" +STILTdir<-"/uufs/chpc.utah.edu/common/home/lin-group7/jcl/CH4_inversion/CH4_inversion_Uintah/CH4_inversion_Uintah_STILT/out/" +resultname<-paste("Fch4.",SITE,".rds",sep="") +#XLIMS<-c(-111,-108.6) #longitude range of Uintah Basin +#YLIMS<-c(39.6,40.6) #latitude range of Uintah Basin +XLIMS<-c(-110.6,-109) #longitude range of Uintah Basin +YLIMS<-c(39.9,40.5) #latitude range of Uintah Basin +#################### + +######################################################### +#I. Calculate CH4 enhancements over background site (FRU) +YYYYMM.COVID <- c("202006","202007","202008")[2] # COVID period, in YYYYMM format +MMsel <- substring(YYYYMM.COVID,5,6) +dat.all<-NULL +for(i in 1:length(YEARs)){ + objname<-paste0("SimCity_CH4_allsites_hrly_",YEARs[i]) + print(paste("Reading in.....",objname)) + tmp<-getr(objname,path=obsdir) + if(paste0("CH4_",SITE)%in%colnames(tmp)){ + tmp2 <- tmp[,c("Time",paste0("CH4_",c(SITE,"FRU")))] + } else { + tmp2 <- data.frame(Time=tmp[,"Time"],NA,CH4_FRU=tmp[,"CH4_FRU"]) + colnames(tmp2)[2] <- paste0("CH4_",SITE) + } # if(paste0("CH4_",SITE)%in%colnames(tmp)){ + dat.all<-rbind(dat.all,tmp2) + gc() +} #for(i in 1:length(YEARs)){ + +dev.new();par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3) +plot(dat.all[,c("Time","CH4_FRU")],main=paste("FRU\nUThrs:",paste(HRs,collapse=" ")),pch=16, + xlab="Time") + +dCH4<-dat.all[,paste0("CH4_",SITE)]-dat.all[,"CH4_FRU"] +dCH4<-data.frame(Time=dat.all[,"Time"],dCH4) #[ppm] +YYYYMMDDHH<-format(dCH4[,"Time"],format="%Y%m%d%H") +sel1<-as.numeric(substring(YYYYMMDDHH,5,6))%in%MONs +sel2<-as.numeric(substring(YYYYMMDDHH,9,10))%in%HRs +dCH4<-data.frame(YYYYMMDDHH,dCH4,stringsAsFactors=FALSE) +dCH4<-dCH4[sel1&sel2,] + + +#calculate DAILY dCH4 +Time.day<-format(dCH4[,"Time"],format="%Y%m%d") +dCH4.ave<-tapply(dCH4[,"dCH4"],Time.day,mean,na.rm=T) +Time<-strptime(names(dCH4.ave),"%Y%m%d",tz="GMT") +dCH4.ave<-data.frame(YYYYMMDD=names(dCH4.ave),Time,dCH4.ave,stringsAsFactors=FALSE) + +dev.new();par(cex.axis=1.3,cex.lab=1.3,cex.main=1.3) +plot(dCH4.ave[,c("Time","dCH4.ave")],main=paste(SITE,"\nUThrs:",paste(HRs,collapse=" ")),pch=16, + xlab="Time",ylab="Daily Averge Enhancement over Background [ppm]") +# plot up months with COVID in 2020 (but also in previous years--want to remove periods with significant enhancements) +sel <- substring(dCH4.ave[,"YYYYMMDD"],5,6)%in%MMsel +points(dCH4.ave[sel,c("Time","dCH4.ave")],pch=16,col="red") +legend("topright",paste("months:",paste(MMsel,collapse=" ")),pch=16,col="red",bty="n",text.col="red",cex=1.3) + +dCH4.ave<-data.frame(dCH4.ave,footsum=NA,foot.basin=NA) +sel<-substring(dCH4.ave[,"YYYYMMDD"],1,6)%in%c("201501","201502","201503","201504","201505") #lack of obs--remove +dCH4.ave<-dCH4.ave[!sel,] + + +######################################################### +#II. Compare DIURNAL CYCLES @ CSP from 2020 COVID period versus those from previous years +YYYYMM.COVID <- c("202003","202004","202005")[-1] # COVID period, in YYYYMM format +MMsel<-substring(YYYYMM.COVID,5,6) +SITE <- "CSP" +dat.all<-NULL +for(i in 1:length(YEARs)){ + objname<-paste0("SimCity_CH4_allsites_hrly_",YEARs[i]) + print(paste("Reading in.....",objname)) + tmp<-getr(objname,path=obsdir) + if(paste0("CH4_",SITE)%in%colnames(tmp)){ + tmp2 <- tmp[,c("Time",paste0("CH4_",c(SITE,"FRU")))] + } else { + tmp2 <- data.frame(Time=tmp[,"Time"],NA,CH4_FRU=tmp[,"CH4_FRU"]) + colnames(tmp2)[2] <- paste0("CH4_",SITE) + } # if(paste0("CH4_",SITE)%in%colnames(tmp)){ + dat.all<-rbind(dat.all,tmp2) + gc() +} #for(i in 1:length(YEARs)){ +YYYYMMDDHH <- format(dat.all[,"Time"],format="%Y%m%d%H") +sel <- substring(YYYYMMDDHH,1,6)%in%c("201501","201502","201503","201504","201505") #lack of obs--remove +dat.all <- dat.all[!sel,] + +if(SITE=="CSP"){ + YYYYsel<-"2016" + YYYYMM.other <- paste0(YYYYsel,MMsel) +} #if(SITE=="CSP"){ + +sitecolnm<-paste0("CH4_",SITE) + +# check which months have most non-NA data +YYYYMMDDHH<-format(dat.all[,"Time"],format="%Y%m%d%H") +YYYYMM<-substring(YYYYMMDDHH,1,6) +sumNotNA<-function(x)return(sum(!is.na(x))) +dum <- tapply(dat.all[,sitecolnm],YYYYMM,sumNotNA) +print(dum) + + +sel<-YYYYMM%in%YYYYMM.COVID +Cdat<-dat.all[sel,] +sel<-YYYYMM%in%YYYYMM.other +Odat<-dat.all[sel,] + +# plot time series at SITE and FRU +YLIMs<-c(1.9,5) +dev.new();par(cex.axis=1.3,cex.main=1.3,cex.lab=1.3) +plot(Cdat[,c("Time",sitecolnm)],pch=16,col=c("orange"),ylim=YLIMs,ylab="CH4 [ppm]") +points(Cdat[,c("Time",paste0("CH4_FRU"))],pch=16,col=c("darkgray")) +legend("topright",pch=16,col=c("orange","darkgray"),legend=c(sitecolnm,"CH4_FRU")) +title(main=paste(YYYYMM.COVID,collapse=" ")) + +dev.new();par(cex.axis=1.3,cex.main=1.3,cex.lab=1.3) +plot(Odat[,c("Time",sitecolnm)],pch=16,col=c("orange"),ylim=YLIMs,ylab="CH4 [ppm]") +points(Odat[,c("Time",paste0("CH4_FRU"))],pch=16,col=c("darkgray")) +legend("topright",pch=16,col=c("orange","darkgray"),legend=c(sitecolnm,"CH4_FRU")) +title(main=paste(YYYYMM.other,collapse=" ")) + +# dev.new(); plot(Odat[1:48,c("Time","CH4_FRU")],pch=16,type="o",col=c("darkgray")); title(main=paste(YYYYMM.other,collapse=" ")) +# dev.new(); plot(Cdat[1:48,c("Time","CH4_FRU")],pch=16,type="o",col=c("darkgray"));title(main=paste(YYYYMM.COVID,collapse=" ")) + +# calculate diurnal cycle of CH4 +diur.C <- tapply(Cdat[,sitecolnm],format(Cdat[,"Time"],format="%H"),mean,na.rm=T) +diur.O <- tapply(Odat[,sitecolnm],format(Odat[,"Time"],format="%H"),mean,na.rm=T) +UThrs <- as.numeric(names(diur.C)) +dev.new();par(cex.axis=1.3,cex.main=1.3,cex.lab=1.3) +plot(UThrs,diur.C,ylim=range(c(diur.C,diur.O)),pch=16,type="o",col="red",xlab="Hour of Day [UT]",ylab="CH4 [ppm]") +lines(UThrs,diur.O,pch=16,type="o",col="black") +title(main=paste("Average Diurnal Cycle at",SITE)) +legend("topright",c(paste(YYYYMM.other,collapse=","),paste(YYYYMM.COVID,collapse=",")), + pch=16,lty=1,text.col=c("black","red"),col=c("black","red"),bty="n") + + +# calculate diurnal cycle of CH4 ENHANCEMENTS (CH4_CSP - CH4_FRU) +diurDCH4.C <- tapply(Cdat[,sitecolnm]-Cdat[,"CH4_FRU"],format(Cdat[,"Time"],format="%H"),mean,na.rm=T) +sigma.C <- tapply(Cdat[,sitecolnm]-Cdat[,"CH4_FRU"],format(Cdat[,"Time"],format="%H"),stdev,na.rm=T) +N.C <- tapply(Cdat[,sitecolnm]-Cdat[,"CH4_FRU"],format(Cdat[,"Time"],format="%H"),sumNotNA) +stderr.C <- sigma.C/sqrt(N.C) + +diurDCH4.O <- tapply(Odat[,sitecolnm]-Odat[,"CH4_FRU"],format(Odat[,"Time"],format="%H"),mean,na.rm=T) +sigma.O <- tapply(Odat[,sitecolnm]-Odat[,"CH4_FRU"],format(Odat[,"Time"],format="%H"),stdev,na.rm=T) +N.O <- tapply(Odat[,sitecolnm]-Odat[,"CH4_FRU"],format(Odat[,"Time"],format="%H"),sumNotNA) +stderr.O <- sigma.O/sqrt(N.O) + +dev.new();par(cex.axis=1.3,cex.main=1.3,cex.lab=1.3) +UThrs <- as.numeric(names(diurDCH4.C)) +plot(UThrs,diurDCH4.C,ylim=range(c(diurDCH4.C,diurDCH4.O)),pch=16, + type="o",col="red",xlab="Hour of Day [UT]",ylab="CH4 enhancement over FRU [ppm]") +segments(UThrs,diurDCH4.C-stderr.C,UThrs,diurDCH4.C+stderr.C,col="red") +lines(UThrs,diurDCH4.O,pch=16,type="o",col="black") +segments(UThrs,diurDCH4.O-stderr.O,UThrs,diurDCH4.O+stderr.O,col="black") +title(main=paste("Average Diurnal Cycle:",SITE,"- FRU")) +legend("topright",c(paste(YYYYMM.other,collapse=","),paste(YYYYMM.COVID,collapse=",")), + pch=16,lty=1,text.col=c("black","red"),col=c("black","red"),bty="n") +# number of non-NA data in different months--make sure there are no systematic differences between 2016 and 2020, since there is difference between March, April, and May +N.C <- tapply(Cdat[,sitecolnm]-Cdat[,"CH4_FRU"],format(Cdat[,"Time"],format="%m"),sumNotNA) +N.O <- tapply(Odat[,sitecolnm]-Odat[,"CH4_FRU"],format(Odat[,"Time"],format="%m"),sumNotNA) +Mons <- as.numeric(names(N.C)) +dev.new();par(cex.axis=1.3,cex.main=1.3,cex.lab=1.3) +plot(Mons,N.C,ylim=range(c(N.C,N.O)),pch=16,type="o",col="red",xlab="Month",ylab="# of non-NA data points") +lines(Mons,N.O,pch=16,type="o",col="black") +legend("topright",c(paste(YYYYMM.other,collapse=","),paste(YYYYMM.COVID,collapse=",")), + pch=16,lty=1,text.col=c("black","red"),col=c("black","red"),bty="n") +title(main=paste("# of non-NA CH4 obs during different months of:",SITE,"- FRU")) + + + +######################################################### +#III. Compare DIURNAL CYCLES @ HPL & WBB from 2020 COVID period versus those from previous years +########## +#YYYYMM.COVID <- c("202003","202004","202005")[-1] # COVID period, in YYYYMM format +YYYYMM.COVID <- c("202006","202007","202008")[2] # COVID period, in YYYYMM format +MMsel <- substring(YYYYMM.COVID,5,6) +SITEs <- c("WBB","HPL","ROO","FRU")[2] +########## + +YEARs <- 2015:2020 +dat.all<-NULL +for(i in 1:length(YEARs)){ + objname<-paste0("SimCity_CH4_allsites_hrly_",YEARs[i]) + print(paste("Reading in.....",objname)) + tmp <- getr(objname,path=obsdir) + sel <- paste0("CH4_",SITEs)%in%colnames(tmp) + if(sum(sel)==length(SITEs)){ + tmp2 <- tmp[,c("Time",paste0("CH4_",c(SITEs,"FRU")))] + } else { + tmp2 <- data.frame(Time=tmp[,"Time"],tmp[,paste0("CH4_",c(SITEs[sel],"FRU"))]) + tmp2 <- data.frame(tmp2,NA) + colnames(tmp2)[ncol(tmp2)] <- paste0("CH4_",SITE[!sel]) + } # if(paste0("CH4_",SITE)%in%colnames(tmp)){ + dat.all<-rbind(dat.all,tmp2) + gc() +} #for(i in 1:length(YEARs)){ + +for(ss in 1:length(SITEs)){ + SITE <- SITEs[ss] + sitecolnm<-paste0("CH4_",SITE) + YYYYMMDDHH <- format(dat.all[,"Time"],format="%Y%m%d%H") + dat <- dat.all[substring(YYYYMMDDHH,5,6)%in%MMsel,] + YYYYMMDDHH <- format(dat[,"Time"],format="%Y%m%d%H") + + dev.new();par(cex.axis=1.3,cex.main=1.3,cex.lab=1.3) + ylims <- c(1.95,2.7) + if(SITE=="WBB") ylims <- c(1.9,2.0) + if(SITE=="FRU") ylims <- c(1.9,2.1) + plot(0,0,type="n",xlab="Hour of Day [UT]",ylab="CH4 [ppm]",ylim=ylims,xlim=c(0,24)) + title(main=paste("Average Diurnal Cycle at",SITE,"\n Months:",paste(MMsel,collapse=","))) + Years <- unique(substring(YYYYMMDDHH,1,4)) + cols <- c("darkblue","lightblue","darkgreen","lightgreen","orange","red") +for(yy in 1:length(Years)){ + sel <- substring(YYYYMMDDHH,1,4)==Years[yy] + diur <- tapply(dat[sel,sitecolnm],format(dat[sel,"Time"],format="%H"),mean,na.rm=T) + UThrs <- as.numeric(names(diur)) + lines(UThrs,diur,pch=16,type="l",col=cols[yy],lwd=2) +} #for(yy in 1:length(Years)){ + legend("topright",Years,bty="n",lwd=2,text.col=cols,col=cols,lty=1) + + dev.new();par(cex.axis=1.3,cex.main=1.3,cex.lab=1.3) + ylims <- c(0,1.0) + if(SITE=="WBB") ylims <- c(0,0.1) + plot(0,0,type="n",xlab="Hour of Day [UT]",ylab="CH4 enhancement over FRU [ppm]",ylim=ylims,xlim=c(0,24)) + title(main=paste("Average Diurnal Cycle:",SITE,"- FRU,\n Months:",paste(MMsel,collapse=","))) + Years <- unique(substring(YYYYMMDDHH,1,4)) + cols <- c("darkblue","lightblue","darkgreen","lightgreen","orange","red") +for(yy in 1:length(Years)){ + sel <- substring(YYYYMMDDHH,1,4)==Years[yy] + diur <- tapply(dat[sel,sitecolnm]-dat[sel,"CH4_FRU"],format(dat[sel,"Time"],format="%H"),mean,na.rm=T) + UThrs <- as.numeric(names(diur)) + lines(UThrs,diur,pch=16,type="l",col=cols[yy],lwd=2) +} #for(yy in 1:length(Years)){ + legend("topright",Years,bty="n",lwd=2,text.col=cols,col=cols,lty=1) +} #for(ss in 1:length(SITEs)){ + +