Nothing
inst/compareOldTool/EvaluationGLPartitioning.R
In REddyProc: Data processing and plotting utilities of (half-)hourly eddy-covariance measurements
#Script for testing the new partitioning GL for REddyProc
library(REddyProc)
# alternative
.tmp.f <- function(){
library(Rcpp)
library(plyr)
source("R/GeoFunctions.R")
source("R/FileHandling.R")
source("R/DataFunctions.R")
source("R/PartitioningLasslop10.R")
pkg <- 'REddyProc'
loadDll <- function (dynFilenameLocal = file.path( system.file(package=pkg), "libs", "x64"
, paste0(pkg, .Platform$dynlib.ext)), pkg = pkg, isWarningDllNotFound = TRUE
){
if (file.exists(dynFilenameLocal)) {
dyn.load(dynFilenameLocal)
}
else if (isWarningDllNotFound) {
warning(paste("dynLib", dynFilenameLocal, "not found."))
}
}
loadDll(pkg=pkg)
source("R/RcppExports.R")
}
library(sirad)
library(scales) # for plotting (function alpha())
## TO DO
path <- "M:/work_3/REddyProcRelease/Eval_GL_Partitioning/"
flist <- list.files(paste0(path,"MR_GL_partitioning/"), pattern="*DataSetafterFluxpart.txt")[-31]
sites <- substr(flist,1,6)
latLongSites <- rbind(
data.frame(site="DE-Tha", lat=51, long=11, timeOffset=-1 )
,data.frame(site="IT-MBo", lat=45.0, long=1, timeOffset=0 )
)
tmp <- read.csv(file.path(path,"CommonAnc.csv"), colClasses=c(Site.ID="character", Latitude="numeric", Longitude="numeric", UTC="numeric"), na.strings="TBD")
latLongSites <- data.frame(site=tmp$Site.ID, lat=tmp$Latitude, long=tmp$Longitude, timeOffset=floor(-tmp$UTC))
## results arrays
# evaluation metrics
metrics <- c("N","pearson","MBE","RMBE","MAE","RMAE","RMSE","RRMSE","R2","slope",
"intercept","EF","SD","CRM","MPE","AC","ACu","ACs")
aggregation <- c("halfhourly","daily","monthly")
vars <- c("GPP","Reco")
NT_vs_DT_REddy <- DT_REddy_vs_pvwave <- array(NA,dim=c(length(sites),length(metrics),length(aggregation),length(vars)),
dimnames=list(sites,metrics,aggregation,vars))
check_quality <- TRUE # plot halfhourly NEE time series as a quality check?
###########################
#---- Reading data
# NEE_orig already ustar filtered!! as in Papale_2006 most conservative threshold,
# the seasons are computed in a different manner (not e.g. JFM)
# gapfilling according to Reichstein_2005
#s <- grep("CA-TP3",sites)[1]
for ( s in seq_along(sites)) {
fname <- flist[s]
year <- as.numeric(substr(fname, nchar(fname)-28, nchar(fname)-25))
fname.PVwave <- paste(sites[s],'.',year,'.','DataSetafterFluxpartGL2010.txt', sep="")
latLongSite <- unlist(subset(latLongSites, site==sites[s])[1,2:4])
#+++ Loading data from MR partitioning and data for running the partitioning
dfall <- fLoadTXTIntoDataframe(fname, paste0(path,"MR_GL_partitioning/"))
dfall.Lass.PVwave <- read.table(paste(path,"MR_GL_partitioning/",fname.PVwave,sep=""),skip=2)
title <- scan(paste(path,"MR_GL_partitioning/",fname.PVwave,sep=""), nlines = 1, sep = "", strip.white=TRUE,
what=list(rep('character',17)))
names(dfall.Lass.PVwave) <- title[[1]]
#+++ Add time stamp in POSIX time format
dfall$PotRad <- as.numeric(fCalcPotRadiation(dfall$julday,dfall$Hour,latLongSite["lat"],latLongSite["long"],latLongSite["timeOffset"]))
dfall$day <- (1 - dfall$night)*100
dfall_posix <- fConvertTimeToPosix(dfall, 'YMDH', Year.s = 'Year', Month.s='Month', Day.s = 'Day', Hour.s = 'Hr')
#+++ Initalize R5 reference class sEddyProc for processing of eddy data
#+++ with all variables needed for processing later
#EddyProc.C <- sEddyProc$new(sites[s], dfall_posix,
# c('NEE', 'NEE_f', 'NEE_fqc', 'Rg', 'Rg_f', 'Rg_fqc','Tair','Tair_fqc','Tsoil',
# 'VPD','VPD_f', 'VPD_fqc','Ustar', "night","day","PotRad"))
# EddyProc.C$sDATA$night
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# START - RUN THE REddyProc DT partitioning
df.REddy <- partitionNEEGL(dfall,NEEVar.s="NEE_f",QFNEEVar.s="NEE_fqc",QFNEEValue.n = 0,NEESdVar.s="NEE_fs_unc",
TempVar.s="Tair_f",QFTempVar.s="Tair_fqc",QFTempValue.n=0,VPDVar.s="VPD_f",QFVPDVar.s="VPD_fqc",
QFVPDValue.n=0,RadVar.s="Rg",PotRadVar.s="day",Suffix.s=""
,controlGLPart.l=partGLControl(nBootUncertainty=0L, isAssociateParmsToMeanOfValids=FALSE, isLasslopPriorsApplied=TRUE,
isBoundLowerNEEUncertainty=FALSE))
.tmp.debug <- function(){
df.REddy$DateTime <- dfall_posix$DateTime
iParRecs <- which( is.finite(df.REddy$FP_qc))
head(df.REddy[iParRecs,], 20)
plot( Reco_DT ~ DateTime, df.REddy)
#trace(partGLFitLRC, recover) #untrace(partGLFitLRC)
}
### add modelled NEE
df.REddy$NEE_DT <- -(df.REddy$GPP_DT - df.REddy$Reco_DT)
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Binding Data frame PVWave tool (all data frame and GL partitioning)
df.Pvwave <- cbind(dfall, dfall.Lass.PVwave[,c(6:17)])
# the columns "qcLE" and "qcH" exist only in DE-Tha. They are deleted here so that
# the data frame can be merged with the other sites (below)
if (sites[s] == "DE-Tha"){
df.Pvwave <- df.Pvwave[,-c(which(colnames(df.Pvwave) %in% c("qcLE","qcH")))]
}
## save data frames resulting from Pvwave and df.REddy
#save(df.REddy,file=paste0(path,"Results/",sites[s],"_df.REddy.RData")) # RData
write.table(df.Pvwave,file=paste0(path,"Results/",sites[s],"_df.Pvwave.txt"),row.names=F,col.names=T)
write.table(df.REddy,file=paste0(path,"Results/",sites[s],"_df.REddy.txt"),row.names=F,col.names=T)
#+++++++++++++++++++++
# Evaluation HH values
### 1) Comparison DT method Pvwave vs. REddyProc
DT_REddy_vs_pvwave[s,,"halfhourly","GPP"] <- c(unlist(modeval(df.Pvwave$GPP_HBLR, df.REddy$GPP_DT)))
DT_REddy_vs_pvwave[s,,"halfhourly","Reco"] <- c(unlist(modeval(df.Pvwave$Reco_HBLR, df.REddy$Reco_DT)))
### 2) nighttime vs. daytime in REddyProc
## nighttime REddyProc still missing!!!
# NT_vs_DT_REddy[s,,"GPP"] <- c(unlist(modeval(df.REddy$GPP_f, df.REddy$GPP_DT)))
# NT_vs_DT_REddy[s,,"Reco"] <- c(unlist(modeval(df.REddy$Reco_f, df.REddy$Reco_DT)))
#+++++++++++++++++++++++++++++
## add a few columns to the data frames for evaluation purposes
## the _agg columns are used for aggregation in aggregate() funciton below
df.Pvwave$julday_agg <- c(1,df.Pvwave$julday[1:(nrow(df.Pvwave)-1)])
df.Pvwave$Month_agg <- c(1,df.Pvwave$Month[1:(nrow(df.Pvwave)-1)])
df.Pvwave$Year_agg <- c(df.Pvwave$Year[1],df.Pvwave$Year[1:(nrow(df.Pvwave)-1)])
# Aggregation daily
if(nrow(df.REddy) != nrow(df.Pvwave)) stop("REddy vs Pvwave: row numbers do not match!")
df.Pvwave.dd <- aggregate(df.Pvwave,by=list(df.Pvwave$julday_agg),mean,na.rm=T)
df.REddy.dd <- aggregate(df.REddy,by=list(df.Pvwave$julday_agg),mean,na.rm=T)
DT_REddy_vs_pvwave[s,,"daily","GPP"] <- c(unlist(modeval(df.Pvwave.dd$GPP_HBLR, df.REddy.dd$GPP_DT)))
DT_REddy_vs_pvwave[s,,"daily","Reco"] <- c(unlist(modeval(df.Pvwave.dd$Reco_HBLR, df.REddy.dd$Reco_DT)))
#+++++++++++++++++++++++++++++
# Aggregation Monthly per site
df.Pvwave.mm <- aggregate(df.Pvwave,by=list(df.Pvwave$Month_agg),mean,na.rm=T)
df.REddy.mm <- aggregate(df.REddy,by=list(df.Pvwave$Month_agg),mean,na.rm=T)
df.Pvwave.mm <- cbind(sites[s],df.Pvwave.mm[,-1]); colnames(df.Pvwave.mm)[1] <- "Site"
df.REddy.mm <- cbind(sites[s],df.REddy.mm[,-1],df.Pvwave.mm$Month_agg); colnames(df.REddy.mm)[c(1,ncol(df.REddy.mm))] <- c("Site","Month_agg")
DT_REddy_vs_pvwave[s,,"monthly","GPP"] <- c(unlist(modeval(df.Pvwave.mm$GPP_HBLR, df.REddy.mm$GPP_DT)))
DT_REddy_vs_pvwave[s,,"monthly","Reco"] <- c(unlist(modeval(df.Pvwave.mm$Reco_HBLR, df.REddy.mm$Reco_DT)))
if (s == 1){
Pvwave.mm.all <- df.Pvwave.mm
REddy.mm.all <- df.REddy.mm
} else {
Pvwave.mm.all <- rbind(Pvwave.mm.all,df.Pvwave.mm)
REddy.mm.all <- rbind(REddy.mm.all,df.REddy.mm)
}
#+++++++++++++++++++++++++++++
# Aggregation Annual
df.Pvwave.yy <- aggregate(df.Pvwave,by=list(df.Pvwave$Year_agg),mean,na.rm=T)
df.REddy.yy <- aggregate(df.REddy,by=list(df.Pvwave$Year_agg),mean,na.rm=T)
df.Pvwave.yy <- cbind(sites[s],df.Pvwave.yy[,-1]); colnames(df.Pvwave.yy)[1] <- "Site"
df.REddy.yy <- cbind(sites[s],df.REddy.yy[,-1]); colnames(df.REddy.yy)[1] <- "Site"
if (s == 1){
Pvwave.yy.all <- df.Pvwave.yy
REddy.yy.all <- df.REddy.yy
} else {
Pvwave.yy.all <- rbind(Pvwave.yy.all,df.Pvwave.yy)
REddy.yy.all <- rbind(REddy.yy.all,df.REddy.yy)
}
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
### Plots
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## 1) data quality check
if (check_quality){
par(mfrow=c(2,2),oma=c(2,2,3,1),mar=c(4,4,1,1))
cex <- 0.9
plot(df.Pvwave$NEE_f,xlab="timestep",ylab="NEE_f (umol m-2 s-1)",las=1,pch=1,col="black",cex=cex)
points(df.Pvwave$NEEorig,col="blue")
legend("topleft",legend="NEEorig",col="blue",pch=1,bty="n",x.intersp=0.5,pt.lwd=2)
plot(df.Pvwave$Tair_f,xlab="timestep",ylab="Tair_f (degC)",las=1,pch=1,col="black",cex=cex)
plot(df.Pvwave$Rg_f,xlab="timestep",ylab="Rg_f (W m-2)",las=1,pch=1,col="black",cex=cex)
plot(df.Pvwave$VPD_f,xlab="timestep",ylab="VPD_f (hPa)",las=1,pch=1,col="black",cex=cex)
mtext(side=3,line=1,sites[s],cex=1.2,outer=T)
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_quality_check.pdf"),width=10,height=8,pointsize=11)
}
## 2) Scatterplots --> the good quality data are added on top of the others in a different color
col.pvwave <- "green3"
col.reddy <- "grey60"
col.reddy_good <- "black"
graphics.off()
par(mfrow=c(1,3),oma=c(0.5,0.5,3.5,0.5),mar=c(4,5,1,1))
cex.pt <- 1.2
transp <- 1 # transparency
pch <- 20
# halfhourly
plot(df.Pvwave$GPP_HBLR ~ df.REddy$GPP_DT,xlab="GPP_DT_REddyProc",ylab="GPP_DT_PVwave",las=1,
main="halfhourly",pch=pch,col=alpha(col.reddy,transp),cex=cex.pt)
points(df.Pvwave$GPP_HBLR[df.REddy$FP_qc < 1] ~ df.REddy$GPP_DT[df.REddy$FP_qc < 1],pch=pch,col=alpha(col.reddy_good,transp),cex=cex.pt)
legend("topleft",legend=paste0("R^2 = ",round(DT_REddy_vs_pvwave[s,"R2","halfhourly","GPP"],2)),bty="n",cex=0.9)
legend(y=0.1*max(df.Pvwave$GPP_HBLR,na.rm=T),x=0.35*max(df.REddy$GPP_DT,na.rm=T),legend=c("all","high quality"),
col=c(col.reddy,col.reddy_good),bty="n",pch=pch,x.intersp=0.5)
curve(1*x,from=-20,to=100,col="red",add=T)
# daily
plot(df.Pvwave.dd$GPP_HBLR ~ df.REddy.dd$GPP_DT,xlab="GPP_DT_REddyProc",ylab="GPP_DT_PVwave",las=1,
main="daily",pch=pch,col=alpha(col.reddy,transp),cex=cex.pt)
points(df.Pvwave.dd$GPP_HBLR[df.REddy.dd$FP_qc < 1] ~ df.REddy.dd$GPP_DT[df.REddy.dd$FP_qc < 1],pch=pch,col=alpha(col.reddy_good,transp),cex=cex.pt)
legend("topleft",legend=paste0("R^2 = ",round(DT_REddy_vs_pvwave[s,"R2","daily","GPP"],2)),bty="n",cex=0.9)
curve(1*x,from=-20,to=100,col="red",add=T)
mtext(side=3,line=2.2,sites[s],cex=1.1)
# monthly
plot(df.Pvwave.mm$GPP_HBLR ~ df.REddy.mm$GPP_DT,xlab="GPP_DT_REddyProc",ylab="GPP_DT_PVwave",las=1,
main="monthly",pch=1,col="black")
legend("topleft",legend=paste0("R^2 = ",round(DT_REddy_vs_pvwave[s,"R2","monthly","GPP"],2)),bty="n",cex=0.9)
curve(1*x,from=-20,to=100,col="red",add=T)
# write to file
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_GPP_DT_PVwave_REddyProc.pdf"),
width=9.5,height=5,pointsize=11)
#### Scatterplots of NEE
## a) comparison between Pvwave and REddyProc
graphics.off()
transp <- 1
par(mfrow=c(1,2),mar=c(5,5,1,0.5),oma=c(0.5,0.5,3,0.5))
daytime <- df.Pvwave$day > 50 # condition: daytime
good <- df.Pvwave$NEE_fqc < 1
# all data (daytime and nighttime)
plot(df.Pvwave$NEE_HBLR ~ df.REddy$NEE_DT,xlab="NEE_REddyProc",ylab="NEE_PVwave",las=1,
pch=pch,col=alpha(col.reddy,transp),cex=cex.pt,main="NEE (all data)",mgp=c(4,1,0))
points(df.Pvwave$NEE_HBLR[good] ~ df.REddy$NEE_DT[good],col=col.reddy_good,pch=pch)
curve(1*x,from=-100,to=100,col="red",add=T)
#mtext(side=3,line=2,text=sites[s])
title(main=sites[s],outer=T,cex=1.5)
# all data (only daytime)
plot(df.Pvwave$NEE_HBLR[daytime] ~ df.REddy$NEE_DT[daytime],xlab="NEE_REddyProc",ylab="NEE_PVwave",las=1,
pch=pch,col=alpha(col.reddy,transp),cex=cex.pt,main="NEE (daytime only)",mgp=c(4,1,0))
points(df.Pvwave$NEE_HBLR[daytime & good] ~ df.REddy$NEE_DT[daytime & good],col=col.reddy_good,pch=pch)
curve(1*x,from=-100,to=100,col="red",add=T)
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_NEE_PVwave_REddyProc.pdf"),
width=9.5,height=5,pointsize=11)
## b) comparison between the products and measurements
graphics.off()
transp <- 1
pch <- 1
par(mfrow=c(1,2),mar=c(5,5,1,0.5),oma=c(0.5,0.5,3,0.5))
xlim <- c(min(df.Pvwave$NEE_HBLR,df.REddy$NEE_DT,na.rm=T),max(df.Pvwave$NEE_HBLR,df.REddy$NEE_DT,na.rm=T))
ylim <- c(min(df.Pvwave$NEE_f,na.rm=T),max(df.Pvwave$NEE_f,na.rm=T))
# Pvwave
plot(df.Pvwave$NEE_f ~ df.Pvwave$NEE_HBLR,xlab="NEE_PVwave",ylab="NEE_f",las=1,xlim=xlim,ylim=ylim,
pch=pch,col=alpha("black",transp),cex=cex.pt)
curve(1*x,from=-100,to=100,col="red",add=T)
title(main=sites[s],outer=T,cex=1.5)
modPvwave <- lm(df.Pvwave$NEE_f ~ df.Pvwave$NEE_HBLR)
abline(modPvwave,col="blue")
mtext(side=3,line=-1.2,paste0("R^2 = ",round(summary(modPvwave)$r.squared,2)))
# REddyProc
plot(df.Pvwave$NEE_f ~ df.REddy$NEE_DT,xlab="NEE_REddyProc",ylab="NEE_f",las=1,xlim=xlim,ylim=ylim,
pch=pch,col=alpha("black",transp),cex=cex.pt)
curve(1*x,from=-100,to=100,col="red",add=T)
modREddy <- lm(df.Pvwave$NEE_f ~ df.REddy$NEE_DT)
abline(modREddy,col="blue")
mtext(side=3,line=-1.2,paste0("R^2 = ",round(summary(modREddy)$r.squared,2)))
#points(df.Pvwave$NEE_f[good] ~ df.REddy$NEE_DT[good],col="black",pch=20)
# Nighttime method
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_NEE_obs_mod.pdf"),
width=9.5,height=5,pointsize=11)
## 3) Timeseries of GPP and Reco
graphics.off()
par(mfrow=c(1,2),mar=c(5,5,1,0.5),oma=c(0.5,0.5,3,0.5))
GPP_good <- df.REddy[,"GPP_DT"]
Reco_good <- df.REddy[,"Reco_DT"]
GPP_good[df.REddy[,"FP_qc"] > 0.5] <- NA
Reco_good[df.REddy[,"FP_qc"] > 0.5] <- NA
# GPP
plot(df.Pvwave$GPP_HBLR,col=col.pvwave,xlab="Timestep",ylab=expression("GPP ("*mu*"mol m"^{-2}~"s"^{-1}*")"),las=1,
ylim=c(min(c(df.Pvwave$GPP_HBLR,df.REddy$GPP_DT),na.rm=T),max(df.Pvwave$GPP_HBLR,na.rm=T) + 0.15*max(df.Pvwave$GPP_HBLR,na.rm=T)))
points(df.REddy$GPP_DT,col=col.reddy)
points(GPP_good,col=col.reddy_good)
legend("topleft",legend=c("REddyProc (all)","REddyProc (high quality)","Pvwave"),col=c(col.reddy,col.reddy_good,col.pvwave),
bty="n",pch=1,x.intersp=0.5,y.intersp=0.8,pt.lwd=3)
mtext(side=3,line=1,sites[s],cex=1.2,outer=T)
# Reco
plot(df.Pvwave$Reco_HBLR,col=col.pvwave,xlab="Timestep",ylab=expression("Reco ("*mu*"mol m"^{-2}~"s"^{-1}*")"),las=1,
ylim=c(min(c(df.Pvwave$Reco_HBLR,df.REddy$Reco_DT),na.rm=T),max(df.Pvwave$Reco_HBLR,na.rm=T) + 0.15*max(df.Pvwave$Reco_HBLR,na.rm=T)))
points(df.REddy$Reco_DT,col=col.reddy)
points(Reco_good,col=col.reddy_good)
# write to file
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_GPP_Reco_DT_PVwave_REddyProc_timeseries.pdf"),
width=11,height=7,pointsize=11)
## 4) mean diurnal courses of GPP and Reco
## data selection --> only summer months .. ok in this case, as all sites are on the northern hemisphere
df.Pvwave.summer <- df.Pvwave[df.Pvwave[,"Month"] %in% c(5,6,7,8,9),]
df.REddy.summer <- df.REddy[df.Pvwave[,"Month"] %in% c(5,6,7,8,9),]
mean_diurnal_Pvwave <- aggregate(df.Pvwave.summer,by=list(df.Pvwave.summer[,9]),mean,na.rm=T) # 9 should be the second "Hour" column
mean_diurnal_REddyProc <- aggregate(df.REddy.summer,by=list(df.Pvwave.summer[,9]),mean,na.rm=T)
graphics.off()
par(mfrow=c(1,2),mar=c(5,5,1,0.5),oma=c(0.5,0.5,3,0.5))
plot(mean_diurnal_Pvwave$GPP_HBLR ~ mean_diurnal_Pvwave[,10],col=col.pvwave,xlab="Hour",ylab=expression("GPP ("*mu*"mol m"^{-2}~"s"^{-1}*")"),las=1,type="l",lwd=2,
ylim=c(min(c(mean_diurnal_Pvwave$GPP_HBLR,mean_diurnal_REddyProc$GPP_DT),na.rm=T),max(mean_diurnal_Pvwave$GPP_HBLR,na.rm=T) + 0.15*max(mean_diurnal_Pvwave$GPP_HBLR,na.rm=T)))
points(mean_diurnal_REddyProc$GPP_DT ~ mean_diurnal_Pvwave[,10],col=col.reddy_good,type="l",lwd=2)
legend("topleft",legend=c("REddyProc","Pvwave"),col=c(col.reddy_good,col.pvwave),
bty="n",lty=1,x.intersp=0.5,y.intersp=0.8,lwd=3)
mtext(side=3,line=1,sites[s],cex=1.2,outer=T)
plot(mean_diurnal_Pvwave$Reco_HBLR ~ mean_diurnal_Pvwave[,10],col=col.pvwave,xlab="Hour",ylab=expression("Reco ("*mu*"mol m"^{-2}~"s"^{-1}*")"),las=1,type="l",lwd=2,
ylim=c(min(c(mean_diurnal_Pvwave$Reco_HBLR,mean_diurnal_REddyProc$Reco_DT),na.rm=T),max(mean_diurnal_Pvwave$Reco_HBLR,na.rm=T) + 0.15*max(mean_diurnal_Pvwave$Reco_HBLR,na.rm=T)))
points(mean_diurnal_REddyProc$Reco_DT ~ mean_diurnal_Pvwave[,10],col=col.reddy_good,type="l",lwd=2)
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_GPP_Reco_DT_PVwave_REddyProc_mean_diurnal_courses_May_September.pdf"),
width=11,height=7,pointsize=11)
## 5) Timeseries of the parameters (Pvwave and REddy)
graphics.off()
par(mfrow=c(2,3),mar=c(5,5,1,1),oma=c(0.5,0.5,3,0.5))
pars_pvwave <- df.Pvwave[,c("Rrefopt_OrdE0_2_from","rb","beta","k","E0","alpha")]
pars_pvwave[pars_pvwave < -9000] <- NA
# str(df.REddy)
# R_ref
plot(df.REddy$FP_R_ref,col="black",xlab="timestep",ylab="R_ref",las=1,ylim=c(min(pars_pvwave$rb,df.REddy$FP_R_ref,na.rm=T),
max(pars_pvwave$rb,df.REddy$FP_R_ref,na.rm=T) + 0.2*max(pars_pvwave$rb,df.REddy$FP_R_ref,na.rm=T))) # leave some extra space for legend
points(pars_pvwave$rb,col=col.pvwave)
legend("topleft",legend=c("REddyProc","Pvwave"),col=c("black",col.pvwave),bty="n",pch=1,x.intersp=0.5)
# R_ref_night (for REddyProc only at the moment)
plot(df.REddy$FP_R_refNight,col="black",xlab="timestep",ylab="R_refNight",las=1,ylim=c(min(df.REddy$FP_R_ref,na.rm=T),
max(df.REddy$FP_R_ref,na.rm=T) + 0.2*max(df.REddy$FP_R_ref,na.rm=T)))
#points(pars_pvwave$Rrefopt_OrdE0_2_from,col=col.pvwave) # is it the right parameter?
# E0
plot(df.REddy$FP_E0,col="black",xlab="timestep",ylab="E0",las=1,ylim=c(min(pars_pvwave$E0,df.REddy$FP_E0,na.rm=T),
max(pars_pvwave$E0,df.REddy$FP_E0,na.rm=T)))
points(pars_pvwave$E0,col=col.pvwave)
# alpha
plot(df.REddy$FP_alpha,col="black",xlab="timestep",ylab="alpha",las=1,ylim=c(min(pars_pvwave$alpha,df.REddy$FP_alpha,na.rm=T),
max(pars_pvwave$alpha,df.REddy$FP_alpha,na.rm=T)))
points(pars_pvwave$alpha,col=col.pvwave)
# beta
plot(df.REddy$FP_beta,col="black",xlab="timestep",ylab="beta",las=1,ylim=c(min(pars_pvwave$beta,df.REddy$FP_beta,na.rm=T),
max(pars_pvwave$beta,df.REddy$FP_beta,na.rm=T)))
points(pars_pvwave$beta,col=col.pvwave)
# k
plot(df.REddy$FP_k,col="black",xlab="timestep",ylab="k",las=1,ylim=c(min(pars_pvwave$k,df.REddy$FP_k,na.rm=T),
max(pars_pvwave$k,df.REddy$FP_k,na.rm=T)))
points(pars_pvwave$k,col=col.pvwave)
mtext(side=3,line=1,sites[s],cex=1.2,outer=T)
# write to file
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_Parameters_Timeseries.pdf"),
width=9,height=6.5,pointsize=11)
## 6) Scattterplot of the parameters (Pvwave and REddy)
graphics.off()
par(mfrow=c(2,3),mar=c(5,5,1,1),oma=c(0.5,0.5,3,0.5),mgp=c(3,0.5,0))
pars_pvwave <- df.Pvwave[,c("rb","beta","k","E0","alpha")]
pars_pvwave[pars_pvwave < -9000] <- NA
# the parameter values of REddyProc and Pvwave are shifted one value towards each other
# temporary fix is to shift the Pvwave values one timestep forward (by deleting the first value):
pars_pvwave <- rbind(pars_pvwave[-1,],NA)
pars_pvwave <- rbind(pars_pvwave[-1,],NA)
# R_ref
plot(df.REddy$FP_R_ref ~ pars_pvwave$rb,xlab="Rb_Pvwave",ylab="Rb_REddy",las=1)
curve(1*x,from=-1000,to=1000,col="red",add=T)
# E0
plot(df.REddy$FP_E0 ~ pars_pvwave$E0,xlab="E0_Pvwave",ylab="E0_REddy",las=1)
curve(1*x,from=-1000,to=1000,col="red",add=T)
# alpha
plot(df.REddy$FP_alpha ~ pars_pvwave$alpha,xlab="alpha_Pvwave",ylab="alpha_REddy",las=1)
curve(1*x,from=-1000,to=1000,col="red",add=T)
# beta
plot(df.REddy$FP_beta ~ pars_pvwave$beta,xlab="beta_Pvwave",ylab="beta_REddy",las=1)
curve(1*x,from=-1000,to=1000,col="red",add=T)
# k
plot(df.REddy$FP_k ~ pars_pvwave$k,xlab="k_Pvwave",ylab="k_REddy",las=1)
curve(1*x,from=-1000,to=1000,col="red",add=T)
mtext(side=3,line=1,sites[s],cex=1.2,outer=T)
# write to file
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_Parameters_Scatterplot.pdf"),
width=9,height=6.5,pointsize=11)
## 7) Comparison to nighttime approach
graphics.off()
par(mfrow=c(2,2),oma=c(1,1,2.5,0.5))
col.nt <- "orange"
lim.reco <- c(min(df.Pvwave$Reco,df.Pvwave$Reco_HBLR,df.REddy$Reco_DT,na.rm=T),
max(df.Pvwave$Reco,df.Pvwave$Reco_HBLR,df.REddy$Reco_DT,na.rm=T))
lim.gpp <- c(min(df.Pvwave$GPP_f,df.Pvwave$GPP_HBLR,df.REddy$GPP_DT,na.rm=T),
max(df.Pvwave$GPP_f,df.Pvwave$GPP_HBLR,df.REddy$GPP_DT,na.rm=T))
# plot(df.Pvwave$Reco,col=col.nt,ylim=ylim)
# points(df.Pvwave$Reco_HBLR,col=col.pvwave)
# points(df.REddy$Reco_DT,col=col.reddy_good)
### Reco
plot(df.Pvwave$Reco ~ df.Pvwave$Reco_HBLR,ylim=lim.reco,xlim=lim.reco,xlab="Daytime_Pvwave",ylab="Nighttime_MR",main="Reco")
curve(1*x,from=-1000,to=1000,col="red",add=T)
mod.pvwave <- lm(df.Pvwave$Reco ~ df.Pvwave$Reco_HBLR)
mtext(side=3,line=-1.8,paste0("R^2 = ",round(summary(mod.pvwave)$r.squared,2)))
title(main=sites[s],outer=T,cex=2)
plot(df.Pvwave$Reco ~ df.REddy$Reco_DT,ylim=lim.reco,xlim=lim.reco,xlab="Daytime_REddyProc",ylab="Nighttime_MR",main="Reco")
curve(1*x,from=-1000,to=1000,col="red",add=T)
mod.reddy <- lm(df.Pvwave$Reco ~ df.REddy$Reco_DT)
mtext(side=3,line=-1.8,paste0("R^2 = ",round(summary(mod.reddy)$r.squared,2)))
### GPP
plot(df.Pvwave$GPP_f ~ df.Pvwave$GPP_HBLR,ylim=lim.gpp,xlim=lim.gpp,xlab="Daytime_Pvwave",ylab="Nighttime_MR",main="GPP")
curve(1*x,from=-1000,to=1000,col="red",add=T)
mod.pvwave <- lm(df.Pvwave$GPP_f ~ df.Pvwave$GPP_HBLR)
mtext(side=3,line=-1.8,paste0("R^2 = ",round(summary(mod.pvwave)$r.squared,2)))
title(main=sites[s],outer=T,cex=1.5)
plot(df.Pvwave$GPP_f ~ df.REddy$GPP_DT,ylim=lim.gpp,xlim=lim.gpp,xlab="Daytime_REddyProc",ylab="Nighttime_MR",main="GPP")
curve(1*x,from=-1000,to=1000,col="red",add=T)
mod.reddy <- lm(df.Pvwave$GPP_f ~ df.REddy$GPP_DT)
mtext(side=3,line=-1.8,paste0("R^2 = ",round(summary(mod.reddy)$r.squared,2)))
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_Reco_GPP_DT_NT_scatterplots.pdf"),
width=9.5,height=9.5,pointsize=11)
# Reco and GPP timeseries
graphics.off()
par(mfrow=c(1,2),oma=c(1,1,2.5,0.5))
plot(df.Pvwave$Reco,col=col.nt,ylim=lim.reco,type="l",xlab="timestep",ylab="Reco",las=1)
points(df.Pvwave$Reco_HBLR,col=col.pvwave,type="l")
points(df.REddy$Reco_DT,col=col.reddy_good,type="l")
title(main=sites[s],outer=T,cex=1.5)
legend("topleft",legend=c("Nighttime","DT_Pvwave","DT_REddyProc"),col=c(col.nt,col.pvwave,col.reddy_good),
lty=1,bty="n",x.intersp=0.5,y.intersp=0.5,cex=0.7,seg.len=0.7)
plot(df.Pvwave$GPP_f,col=col.nt,ylim=lim.gpp,type="l",xlab="timestep",ylab="GPP",las=1)
points(df.Pvwave$GPP_HBLR,col=col.pvwave,type="l")
points(df.REddy$GPP_DT,col=col.reddy_good,type="l")
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_Reco_GPP_DT_NT_timeseries.pdf"),
width=9.5,height=5,pointsize=11)
} # end site loop
#+++++++++++++++++++++++++++++
# Evaluation Monthly All
graphics.off()
par(mfrow=c(2,1),mar=c(4,4,1,1),mgp=c(1.8,0.4,0))
transp <- 0.5
r2_monthly <- round(summary(lm(Pvwave.mm.all$GPP_HBLR ~ REddy.mm.all$GPP_DT))$r.squared,3)
plot(Pvwave.mm.all$GPP_HBLR ~ REddy.mm.all$GPP_DT,xlab="GPP_REddyProc",ylab="GPP_Pvwave",
pch=20,col=alpha("black",transp),las=1,tcl=-0.2,xlim=c(0,20),ylim=c(0,20),cex=1.3)
curve(1*x,from=-10,to=100,col="red",add=T)
legend("bottomright","monthly",bty="n")
legend("topleft",paste0("R^2 = ",r2_monthly),bty="n")
#++++++++++++++++++++++++++++
# Evaluation Annual
r2_yearly <- round(summary(lm(Pvwave.yy.all$GPP_HBLR ~ REddy.yy.all$GPP_DT))$r.squared,3)
plot(Pvwave.yy.all$GPP_HBLR ~ REddy.yy.all$GPP_DT,xlab="GPP_REddyProc",ylab="GPP_Pvwave",
pch=20,col=alpha("black",transp),las=1,tcl=-0.2,xlim=c(0,11),ylim=c(0,11),cex=1.3)
curve(1*x,from=-10,to=100,col="red",add=T)
legend("bottomright","yearly",bty="n")
legend("topleft",paste0("R^2 = ",r2_yearly),bty="n")
dev.copy2pdf(file=paste0(path,"Plots/all_sites_monthly_yearly.pdf"),
width=7.5,height=8,pointsize=11)
## save as RData:
save(NT_vs_DT_REddy,DT_REddy_vs_pvwave,file=paste0(path,"Results/eval_metrics.RData")) # 1) evaluation metrics
save(Pvwave.mm.all,REddy.mm.all,file=paste0(path,"Results/all_sites_monthly.RData")) # 2) monthly aggregated results for all sites
save(Pvwave.yy.all,REddy.yy.all,file=paste0(path,"Results/all_sites_annual.RData")) # 3) annual aggregated results for all sites
Try the REddyProc package in your browser
Any scripts or data that you put into this service are public.
REddyProc documentation built on May 2, 2019, 5:19 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#Script for testing the new partitioning GL for REddyProc
library(REddyProc)
# alternative
.tmp.f <- function(){
library(Rcpp)
library(plyr)
source("R/GeoFunctions.R")
source("R/FileHandling.R")
source("R/DataFunctions.R")
source("R/PartitioningLasslop10.R")
pkg <- 'REddyProc'
loadDll <- function (dynFilenameLocal = file.path( system.file(package=pkg), "libs", "x64"
, paste0(pkg, .Platform$dynlib.ext)), pkg = pkg, isWarningDllNotFound = TRUE
){
if (file.exists(dynFilenameLocal)) {
dyn.load(dynFilenameLocal)
}
else if (isWarningDllNotFound) {
warning(paste("dynLib", dynFilenameLocal, "not found."))
}
}
loadDll(pkg=pkg)
source("R/RcppExports.R")
}
library(sirad)
library(scales) # for plotting (function alpha())
## TO DO
path <- "M:/work_3/REddyProcRelease/Eval_GL_Partitioning/"
flist <- list.files(paste0(path,"MR_GL_partitioning/"), pattern="*DataSetafterFluxpart.txt")[-31]
sites <- substr(flist,1,6)
latLongSites <- rbind(
data.frame(site="DE-Tha", lat=51, long=11, timeOffset=-1 )
,data.frame(site="IT-MBo", lat=45.0, long=1, timeOffset=0 )
)
tmp <- read.csv(file.path(path,"CommonAnc.csv"), colClasses=c(Site.ID="character", Latitude="numeric", Longitude="numeric", UTC="numeric"), na.strings="TBD")
latLongSites <- data.frame(site=tmp$Site.ID, lat=tmp$Latitude, long=tmp$Longitude, timeOffset=floor(-tmp$UTC))
## results arrays
# evaluation metrics
metrics <- c("N","pearson","MBE","RMBE","MAE","RMAE","RMSE","RRMSE","R2","slope",
"intercept","EF","SD","CRM","MPE","AC","ACu","ACs")
aggregation <- c("halfhourly","daily","monthly")
vars <- c("GPP","Reco")
NT_vs_DT_REddy <- DT_REddy_vs_pvwave <- array(NA,dim=c(length(sites),length(metrics),length(aggregation),length(vars)),
dimnames=list(sites,metrics,aggregation,vars))
check_quality <- TRUE # plot halfhourly NEE time series as a quality check?
###########################
#---- Reading data
# NEE_orig already ustar filtered!! as in Papale_2006 most conservative threshold,
# the seasons are computed in a different manner (not e.g. JFM)
# gapfilling according to Reichstein_2005
#s <- grep("CA-TP3",sites)[1]
for ( s in seq_along(sites)) {
fname <- flist[s]
year <- as.numeric(substr(fname, nchar(fname)-28, nchar(fname)-25))
fname.PVwave <- paste(sites[s],'.',year,'.','DataSetafterFluxpartGL2010.txt', sep="")
latLongSite <- unlist(subset(latLongSites, site==sites[s])[1,2:4])
#+++ Loading data from MR partitioning and data for running the partitioning
dfall <- fLoadTXTIntoDataframe(fname, paste0(path,"MR_GL_partitioning/"))
dfall.Lass.PVwave <- read.table(paste(path,"MR_GL_partitioning/",fname.PVwave,sep=""),skip=2)
title <- scan(paste(path,"MR_GL_partitioning/",fname.PVwave,sep=""), nlines = 1, sep = "", strip.white=TRUE,
what=list(rep('character',17)))
names(dfall.Lass.PVwave) <- title[[1]]
#+++ Add time stamp in POSIX time format
dfall$PotRad <- as.numeric(fCalcPotRadiation(dfall$julday,dfall$Hour,latLongSite["lat"],latLongSite["long"],latLongSite["timeOffset"]))
dfall$day <- (1 - dfall$night)*100
dfall_posix <- fConvertTimeToPosix(dfall, 'YMDH', Year.s = 'Year', Month.s='Month', Day.s = 'Day', Hour.s = 'Hr')
#+++ Initalize R5 reference class sEddyProc for processing of eddy data
#+++ with all variables needed for processing later
#EddyProc.C <- sEddyProc$new(sites[s], dfall_posix,
# c('NEE', 'NEE_f', 'NEE_fqc', 'Rg', 'Rg_f', 'Rg_fqc','Tair','Tair_fqc','Tsoil',
# 'VPD','VPD_f', 'VPD_fqc','Ustar', "night","day","PotRad"))
# EddyProc.C$sDATA$night
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# START - RUN THE REddyProc DT partitioning
df.REddy <- partitionNEEGL(dfall,NEEVar.s="NEE_f",QFNEEVar.s="NEE_fqc",QFNEEValue.n = 0,NEESdVar.s="NEE_fs_unc",
TempVar.s="Tair_f",QFTempVar.s="Tair_fqc",QFTempValue.n=0,VPDVar.s="VPD_f",QFVPDVar.s="VPD_fqc",
QFVPDValue.n=0,RadVar.s="Rg",PotRadVar.s="day",Suffix.s=""
,controlGLPart.l=partGLControl(nBootUncertainty=0L, isAssociateParmsToMeanOfValids=FALSE, isLasslopPriorsApplied=TRUE,
isBoundLowerNEEUncertainty=FALSE))
.tmp.debug <- function(){
df.REddy$DateTime <- dfall_posix$DateTime
iParRecs <- which( is.finite(df.REddy$FP_qc))
head(df.REddy[iParRecs,], 20)
plot( Reco_DT ~ DateTime, df.REddy)
#trace(partGLFitLRC, recover) #untrace(partGLFitLRC)
}
### add modelled NEE
df.REddy$NEE_DT <- -(df.REddy$GPP_DT - df.REddy$Reco_DT)
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Binding Data frame PVWave tool (all data frame and GL partitioning)
df.Pvwave <- cbind(dfall, dfall.Lass.PVwave[,c(6:17)])
# the columns "qcLE" and "qcH" exist only in DE-Tha. They are deleted here so that
# the data frame can be merged with the other sites (below)
if (sites[s] == "DE-Tha"){
df.Pvwave <- df.Pvwave[,-c(which(colnames(df.Pvwave) %in% c("qcLE","qcH")))]
}
## save data frames resulting from Pvwave and df.REddy
#save(df.REddy,file=paste0(path,"Results/",sites[s],"_df.REddy.RData")) # RData
write.table(df.Pvwave,file=paste0(path,"Results/",sites[s],"_df.Pvwave.txt"),row.names=F,col.names=T)
write.table(df.REddy,file=paste0(path,"Results/",sites[s],"_df.REddy.txt"),row.names=F,col.names=T)
#+++++++++++++++++++++
# Evaluation HH values
### 1) Comparison DT method Pvwave vs. REddyProc
DT_REddy_vs_pvwave[s,,"halfhourly","GPP"] <- c(unlist(modeval(df.Pvwave$GPP_HBLR, df.REddy$GPP_DT)))
DT_REddy_vs_pvwave[s,,"halfhourly","Reco"] <- c(unlist(modeval(df.Pvwave$Reco_HBLR, df.REddy$Reco_DT)))
### 2) nighttime vs. daytime in REddyProc
## nighttime REddyProc still missing!!!
# NT_vs_DT_REddy[s,,"GPP"] <- c(unlist(modeval(df.REddy$GPP_f, df.REddy$GPP_DT)))
# NT_vs_DT_REddy[s,,"Reco"] <- c(unlist(modeval(df.REddy$Reco_f, df.REddy$Reco_DT)))
#+++++++++++++++++++++++++++++
## add a few columns to the data frames for evaluation purposes
## the _agg columns are used for aggregation in aggregate() funciton below
df.Pvwave$julday_agg <- c(1,df.Pvwave$julday[1:(nrow(df.Pvwave)-1)])
df.Pvwave$Month_agg <- c(1,df.Pvwave$Month[1:(nrow(df.Pvwave)-1)])
df.Pvwave$Year_agg <- c(df.Pvwave$Year[1],df.Pvwave$Year[1:(nrow(df.Pvwave)-1)])
# Aggregation daily
if(nrow(df.REddy) != nrow(df.Pvwave)) stop("REddy vs Pvwave: row numbers do not match!")
df.Pvwave.dd <- aggregate(df.Pvwave,by=list(df.Pvwave$julday_agg),mean,na.rm=T)
df.REddy.dd <- aggregate(df.REddy,by=list(df.Pvwave$julday_agg),mean,na.rm=T)
DT_REddy_vs_pvwave[s,,"daily","GPP"] <- c(unlist(modeval(df.Pvwave.dd$GPP_HBLR, df.REddy.dd$GPP_DT)))
DT_REddy_vs_pvwave[s,,"daily","Reco"] <- c(unlist(modeval(df.Pvwave.dd$Reco_HBLR, df.REddy.dd$Reco_DT)))
#+++++++++++++++++++++++++++++
# Aggregation Monthly per site
df.Pvwave.mm <- aggregate(df.Pvwave,by=list(df.Pvwave$Month_agg),mean,na.rm=T)
df.REddy.mm <- aggregate(df.REddy,by=list(df.Pvwave$Month_agg),mean,na.rm=T)
df.Pvwave.mm <- cbind(sites[s],df.Pvwave.mm[,-1]); colnames(df.Pvwave.mm)[1] <- "Site"
df.REddy.mm <- cbind(sites[s],df.REddy.mm[,-1],df.Pvwave.mm$Month_agg); colnames(df.REddy.mm)[c(1,ncol(df.REddy.mm))] <- c("Site","Month_agg")
DT_REddy_vs_pvwave[s,,"monthly","GPP"] <- c(unlist(modeval(df.Pvwave.mm$GPP_HBLR, df.REddy.mm$GPP_DT)))
DT_REddy_vs_pvwave[s,,"monthly","Reco"] <- c(unlist(modeval(df.Pvwave.mm$Reco_HBLR, df.REddy.mm$Reco_DT)))
if (s == 1){
Pvwave.mm.all <- df.Pvwave.mm
REddy.mm.all <- df.REddy.mm
} else {
Pvwave.mm.all <- rbind(Pvwave.mm.all,df.Pvwave.mm)
REddy.mm.all <- rbind(REddy.mm.all,df.REddy.mm)
}
#+++++++++++++++++++++++++++++
# Aggregation Annual
df.Pvwave.yy <- aggregate(df.Pvwave,by=list(df.Pvwave$Year_agg),mean,na.rm=T)
df.REddy.yy <- aggregate(df.REddy,by=list(df.Pvwave$Year_agg),mean,na.rm=T)
df.Pvwave.yy <- cbind(sites[s],df.Pvwave.yy[,-1]); colnames(df.Pvwave.yy)[1] <- "Site"
df.REddy.yy <- cbind(sites[s],df.REddy.yy[,-1]); colnames(df.REddy.yy)[1] <- "Site"
if (s == 1){
Pvwave.yy.all <- df.Pvwave.yy
REddy.yy.all <- df.REddy.yy
} else {
Pvwave.yy.all <- rbind(Pvwave.yy.all,df.Pvwave.yy)
REddy.yy.all <- rbind(REddy.yy.all,df.REddy.yy)
}
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
### Plots
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## 1) data quality check
if (check_quality){
par(mfrow=c(2,2),oma=c(2,2,3,1),mar=c(4,4,1,1))
cex <- 0.9
plot(df.Pvwave$NEE_f,xlab="timestep",ylab="NEE_f (umol m-2 s-1)",las=1,pch=1,col="black",cex=cex)
points(df.Pvwave$NEEorig,col="blue")
legend("topleft",legend="NEEorig",col="blue",pch=1,bty="n",x.intersp=0.5,pt.lwd=2)
plot(df.Pvwave$Tair_f,xlab="timestep",ylab="Tair_f (degC)",las=1,pch=1,col="black",cex=cex)
plot(df.Pvwave$Rg_f,xlab="timestep",ylab="Rg_f (W m-2)",las=1,pch=1,col="black",cex=cex)
plot(df.Pvwave$VPD_f,xlab="timestep",ylab="VPD_f (hPa)",las=1,pch=1,col="black",cex=cex)
mtext(side=3,line=1,sites[s],cex=1.2,outer=T)
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_quality_check.pdf"),width=10,height=8,pointsize=11)
}
## 2) Scatterplots --> the good quality data are added on top of the others in a different color
col.pvwave <- "green3"
col.reddy <- "grey60"
col.reddy_good <- "black"
graphics.off()
par(mfrow=c(1,3),oma=c(0.5,0.5,3.5,0.5),mar=c(4,5,1,1))
cex.pt <- 1.2
transp <- 1 # transparency
pch <- 20
# halfhourly
plot(df.Pvwave$GPP_HBLR ~ df.REddy$GPP_DT,xlab="GPP_DT_REddyProc",ylab="GPP_DT_PVwave",las=1,
main="halfhourly",pch=pch,col=alpha(col.reddy,transp),cex=cex.pt)
points(df.Pvwave$GPP_HBLR[df.REddy$FP_qc < 1] ~ df.REddy$GPP_DT[df.REddy$FP_qc < 1],pch=pch,col=alpha(col.reddy_good,transp),cex=cex.pt)
legend("topleft",legend=paste0("R^2 = ",round(DT_REddy_vs_pvwave[s,"R2","halfhourly","GPP"],2)),bty="n",cex=0.9)
legend(y=0.1*max(df.Pvwave$GPP_HBLR,na.rm=T),x=0.35*max(df.REddy$GPP_DT,na.rm=T),legend=c("all","high quality"),
col=c(col.reddy,col.reddy_good),bty="n",pch=pch,x.intersp=0.5)
curve(1*x,from=-20,to=100,col="red",add=T)
# daily
plot(df.Pvwave.dd$GPP_HBLR ~ df.REddy.dd$GPP_DT,xlab="GPP_DT_REddyProc",ylab="GPP_DT_PVwave",las=1,
main="daily",pch=pch,col=alpha(col.reddy,transp),cex=cex.pt)
points(df.Pvwave.dd$GPP_HBLR[df.REddy.dd$FP_qc < 1] ~ df.REddy.dd$GPP_DT[df.REddy.dd$FP_qc < 1],pch=pch,col=alpha(col.reddy_good,transp),cex=cex.pt)
legend("topleft",legend=paste0("R^2 = ",round(DT_REddy_vs_pvwave[s,"R2","daily","GPP"],2)),bty="n",cex=0.9)
curve(1*x,from=-20,to=100,col="red",add=T)
mtext(side=3,line=2.2,sites[s],cex=1.1)
# monthly
plot(df.Pvwave.mm$GPP_HBLR ~ df.REddy.mm$GPP_DT,xlab="GPP_DT_REddyProc",ylab="GPP_DT_PVwave",las=1,
main="monthly",pch=1,col="black")
legend("topleft",legend=paste0("R^2 = ",round(DT_REddy_vs_pvwave[s,"R2","monthly","GPP"],2)),bty="n",cex=0.9)
curve(1*x,from=-20,to=100,col="red",add=T)
# write to file
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_GPP_DT_PVwave_REddyProc.pdf"),
width=9.5,height=5,pointsize=11)
#### Scatterplots of NEE
## a) comparison between Pvwave and REddyProc
graphics.off()
transp <- 1
par(mfrow=c(1,2),mar=c(5,5,1,0.5),oma=c(0.5,0.5,3,0.5))
daytime <- df.Pvwave$day > 50 # condition: daytime
good <- df.Pvwave$NEE_fqc < 1
# all data (daytime and nighttime)
plot(df.Pvwave$NEE_HBLR ~ df.REddy$NEE_DT,xlab="NEE_REddyProc",ylab="NEE_PVwave",las=1,
pch=pch,col=alpha(col.reddy,transp),cex=cex.pt,main="NEE (all data)",mgp=c(4,1,0))
points(df.Pvwave$NEE_HBLR[good] ~ df.REddy$NEE_DT[good],col=col.reddy_good,pch=pch)
curve(1*x,from=-100,to=100,col="red",add=T)
#mtext(side=3,line=2,text=sites[s])
title(main=sites[s],outer=T,cex=1.5)
# all data (only daytime)
plot(df.Pvwave$NEE_HBLR[daytime] ~ df.REddy$NEE_DT[daytime],xlab="NEE_REddyProc",ylab="NEE_PVwave",las=1,
pch=pch,col=alpha(col.reddy,transp),cex=cex.pt,main="NEE (daytime only)",mgp=c(4,1,0))
points(df.Pvwave$NEE_HBLR[daytime & good] ~ df.REddy$NEE_DT[daytime & good],col=col.reddy_good,pch=pch)
curve(1*x,from=-100,to=100,col="red",add=T)
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_NEE_PVwave_REddyProc.pdf"),
width=9.5,height=5,pointsize=11)
## b) comparison between the products and measurements
graphics.off()
transp <- 1
pch <- 1
par(mfrow=c(1,2),mar=c(5,5,1,0.5),oma=c(0.5,0.5,3,0.5))
xlim <- c(min(df.Pvwave$NEE_HBLR,df.REddy$NEE_DT,na.rm=T),max(df.Pvwave$NEE_HBLR,df.REddy$NEE_DT,na.rm=T))
ylim <- c(min(df.Pvwave$NEE_f,na.rm=T),max(df.Pvwave$NEE_f,na.rm=T))
# Pvwave
plot(df.Pvwave$NEE_f ~ df.Pvwave$NEE_HBLR,xlab="NEE_PVwave",ylab="NEE_f",las=1,xlim=xlim,ylim=ylim,
pch=pch,col=alpha("black",transp),cex=cex.pt)
curve(1*x,from=-100,to=100,col="red",add=T)
title(main=sites[s],outer=T,cex=1.5)
modPvwave <- lm(df.Pvwave$NEE_f ~ df.Pvwave$NEE_HBLR)
abline(modPvwave,col="blue")
mtext(side=3,line=-1.2,paste0("R^2 = ",round(summary(modPvwave)$r.squared,2)))
# REddyProc
plot(df.Pvwave$NEE_f ~ df.REddy$NEE_DT,xlab="NEE_REddyProc",ylab="NEE_f",las=1,xlim=xlim,ylim=ylim,
pch=pch,col=alpha("black",transp),cex=cex.pt)
curve(1*x,from=-100,to=100,col="red",add=T)
modREddy <- lm(df.Pvwave$NEE_f ~ df.REddy$NEE_DT)
abline(modREddy,col="blue")
mtext(side=3,line=-1.2,paste0("R^2 = ",round(summary(modREddy)$r.squared,2)))
#points(df.Pvwave$NEE_f[good] ~ df.REddy$NEE_DT[good],col="black",pch=20)
# Nighttime method
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_NEE_obs_mod.pdf"),
width=9.5,height=5,pointsize=11)
## 3) Timeseries of GPP and Reco
graphics.off()
par(mfrow=c(1,2),mar=c(5,5,1,0.5),oma=c(0.5,0.5,3,0.5))
GPP_good <- df.REddy[,"GPP_DT"]
Reco_good <- df.REddy[,"Reco_DT"]
GPP_good[df.REddy[,"FP_qc"] > 0.5] <- NA
Reco_good[df.REddy[,"FP_qc"] > 0.5] <- NA
# GPP
plot(df.Pvwave$GPP_HBLR,col=col.pvwave,xlab="Timestep",ylab=expression("GPP ("*mu*"mol m"^{-2}~"s"^{-1}*")"),las=1,
ylim=c(min(c(df.Pvwave$GPP_HBLR,df.REddy$GPP_DT),na.rm=T),max(df.Pvwave$GPP_HBLR,na.rm=T) + 0.15*max(df.Pvwave$GPP_HBLR,na.rm=T)))
points(df.REddy$GPP_DT,col=col.reddy)
points(GPP_good,col=col.reddy_good)
legend("topleft",legend=c("REddyProc (all)","REddyProc (high quality)","Pvwave"),col=c(col.reddy,col.reddy_good,col.pvwave),
bty="n",pch=1,x.intersp=0.5,y.intersp=0.8,pt.lwd=3)
mtext(side=3,line=1,sites[s],cex=1.2,outer=T)
# Reco
plot(df.Pvwave$Reco_HBLR,col=col.pvwave,xlab="Timestep",ylab=expression("Reco ("*mu*"mol m"^{-2}~"s"^{-1}*")"),las=1,
ylim=c(min(c(df.Pvwave$Reco_HBLR,df.REddy$Reco_DT),na.rm=T),max(df.Pvwave$Reco_HBLR,na.rm=T) + 0.15*max(df.Pvwave$Reco_HBLR,na.rm=T)))
points(df.REddy$Reco_DT,col=col.reddy)
points(Reco_good,col=col.reddy_good)
# write to file
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_GPP_Reco_DT_PVwave_REddyProc_timeseries.pdf"),
width=11,height=7,pointsize=11)
## 4) mean diurnal courses of GPP and Reco
## data selection --> only summer months .. ok in this case, as all sites are on the northern hemisphere
df.Pvwave.summer <- df.Pvwave[df.Pvwave[,"Month"] %in% c(5,6,7,8,9),]
df.REddy.summer <- df.REddy[df.Pvwave[,"Month"] %in% c(5,6,7,8,9),]
mean_diurnal_Pvwave <- aggregate(df.Pvwave.summer,by=list(df.Pvwave.summer[,9]),mean,na.rm=T) # 9 should be the second "Hour" column
mean_diurnal_REddyProc <- aggregate(df.REddy.summer,by=list(df.Pvwave.summer[,9]),mean,na.rm=T)
graphics.off()
par(mfrow=c(1,2),mar=c(5,5,1,0.5),oma=c(0.5,0.5,3,0.5))
plot(mean_diurnal_Pvwave$GPP_HBLR ~ mean_diurnal_Pvwave[,10],col=col.pvwave,xlab="Hour",ylab=expression("GPP ("*mu*"mol m"^{-2}~"s"^{-1}*")"),las=1,type="l",lwd=2,
ylim=c(min(c(mean_diurnal_Pvwave$GPP_HBLR,mean_diurnal_REddyProc$GPP_DT),na.rm=T),max(mean_diurnal_Pvwave$GPP_HBLR,na.rm=T) + 0.15*max(mean_diurnal_Pvwave$GPP_HBLR,na.rm=T)))
points(mean_diurnal_REddyProc$GPP_DT ~ mean_diurnal_Pvwave[,10],col=col.reddy_good,type="l",lwd=2)
legend("topleft",legend=c("REddyProc","Pvwave"),col=c(col.reddy_good,col.pvwave),
bty="n",lty=1,x.intersp=0.5,y.intersp=0.8,lwd=3)
mtext(side=3,line=1,sites[s],cex=1.2,outer=T)
plot(mean_diurnal_Pvwave$Reco_HBLR ~ mean_diurnal_Pvwave[,10],col=col.pvwave,xlab="Hour",ylab=expression("Reco ("*mu*"mol m"^{-2}~"s"^{-1}*")"),las=1,type="l",lwd=2,
ylim=c(min(c(mean_diurnal_Pvwave$Reco_HBLR,mean_diurnal_REddyProc$Reco_DT),na.rm=T),max(mean_diurnal_Pvwave$Reco_HBLR,na.rm=T) + 0.15*max(mean_diurnal_Pvwave$Reco_HBLR,na.rm=T)))
points(mean_diurnal_REddyProc$Reco_DT ~ mean_diurnal_Pvwave[,10],col=col.reddy_good,type="l",lwd=2)
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_GPP_Reco_DT_PVwave_REddyProc_mean_diurnal_courses_May_September.pdf"),
width=11,height=7,pointsize=11)
## 5) Timeseries of the parameters (Pvwave and REddy)
graphics.off()
par(mfrow=c(2,3),mar=c(5,5,1,1),oma=c(0.5,0.5,3,0.5))
pars_pvwave <- df.Pvwave[,c("Rrefopt_OrdE0_2_from","rb","beta","k","E0","alpha")]
pars_pvwave[pars_pvwave < -9000] <- NA
# str(df.REddy)
# R_ref
plot(df.REddy$FP_R_ref,col="black",xlab="timestep",ylab="R_ref",las=1,ylim=c(min(pars_pvwave$rb,df.REddy$FP_R_ref,na.rm=T),
max(pars_pvwave$rb,df.REddy$FP_R_ref,na.rm=T) + 0.2*max(pars_pvwave$rb,df.REddy$FP_R_ref,na.rm=T))) # leave some extra space for legend
points(pars_pvwave$rb,col=col.pvwave)
legend("topleft",legend=c("REddyProc","Pvwave"),col=c("black",col.pvwave),bty="n",pch=1,x.intersp=0.5)
# R_ref_night (for REddyProc only at the moment)
plot(df.REddy$FP_R_refNight,col="black",xlab="timestep",ylab="R_refNight",las=1,ylim=c(min(df.REddy$FP_R_ref,na.rm=T),
max(df.REddy$FP_R_ref,na.rm=T) + 0.2*max(df.REddy$FP_R_ref,na.rm=T)))
#points(pars_pvwave$Rrefopt_OrdE0_2_from,col=col.pvwave) # is it the right parameter?
# E0
plot(df.REddy$FP_E0,col="black",xlab="timestep",ylab="E0",las=1,ylim=c(min(pars_pvwave$E0,df.REddy$FP_E0,na.rm=T),
max(pars_pvwave$E0,df.REddy$FP_E0,na.rm=T)))
points(pars_pvwave$E0,col=col.pvwave)
# alpha
plot(df.REddy$FP_alpha,col="black",xlab="timestep",ylab="alpha",las=1,ylim=c(min(pars_pvwave$alpha,df.REddy$FP_alpha,na.rm=T),
max(pars_pvwave$alpha,df.REddy$FP_alpha,na.rm=T)))
points(pars_pvwave$alpha,col=col.pvwave)
# beta
plot(df.REddy$FP_beta,col="black",xlab="timestep",ylab="beta",las=1,ylim=c(min(pars_pvwave$beta,df.REddy$FP_beta,na.rm=T),
max(pars_pvwave$beta,df.REddy$FP_beta,na.rm=T)))
points(pars_pvwave$beta,col=col.pvwave)
# k
plot(df.REddy$FP_k,col="black",xlab="timestep",ylab="k",las=1,ylim=c(min(pars_pvwave$k,df.REddy$FP_k,na.rm=T),
max(pars_pvwave$k,df.REddy$FP_k,na.rm=T)))
points(pars_pvwave$k,col=col.pvwave)
mtext(side=3,line=1,sites[s],cex=1.2,outer=T)
# write to file
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_Parameters_Timeseries.pdf"),
width=9,height=6.5,pointsize=11)
## 6) Scattterplot of the parameters (Pvwave and REddy)
graphics.off()
par(mfrow=c(2,3),mar=c(5,5,1,1),oma=c(0.5,0.5,3,0.5),mgp=c(3,0.5,0))
pars_pvwave <- df.Pvwave[,c("rb","beta","k","E0","alpha")]
pars_pvwave[pars_pvwave < -9000] <- NA
# the parameter values of REddyProc and Pvwave are shifted one value towards each other
# temporary fix is to shift the Pvwave values one timestep forward (by deleting the first value):
pars_pvwave <- rbind(pars_pvwave[-1,],NA)
pars_pvwave <- rbind(pars_pvwave[-1,],NA)
# R_ref
plot(df.REddy$FP_R_ref ~ pars_pvwave$rb,xlab="Rb_Pvwave",ylab="Rb_REddy",las=1)
curve(1*x,from=-1000,to=1000,col="red",add=T)
# E0
plot(df.REddy$FP_E0 ~ pars_pvwave$E0,xlab="E0_Pvwave",ylab="E0_REddy",las=1)
curve(1*x,from=-1000,to=1000,col="red",add=T)
# alpha
plot(df.REddy$FP_alpha ~ pars_pvwave$alpha,xlab="alpha_Pvwave",ylab="alpha_REddy",las=1)
curve(1*x,from=-1000,to=1000,col="red",add=T)
# beta
plot(df.REddy$FP_beta ~ pars_pvwave$beta,xlab="beta_Pvwave",ylab="beta_REddy",las=1)
curve(1*x,from=-1000,to=1000,col="red",add=T)
# k
plot(df.REddy$FP_k ~ pars_pvwave$k,xlab="k_Pvwave",ylab="k_REddy",las=1)
curve(1*x,from=-1000,to=1000,col="red",add=T)
mtext(side=3,line=1,sites[s],cex=1.2,outer=T)
# write to file
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_Parameters_Scatterplot.pdf"),
width=9,height=6.5,pointsize=11)
## 7) Comparison to nighttime approach
graphics.off()
par(mfrow=c(2,2),oma=c(1,1,2.5,0.5))
col.nt <- "orange"
lim.reco <- c(min(df.Pvwave$Reco,df.Pvwave$Reco_HBLR,df.REddy$Reco_DT,na.rm=T),
max(df.Pvwave$Reco,df.Pvwave$Reco_HBLR,df.REddy$Reco_DT,na.rm=T))
lim.gpp <- c(min(df.Pvwave$GPP_f,df.Pvwave$GPP_HBLR,df.REddy$GPP_DT,na.rm=T),
max(df.Pvwave$GPP_f,df.Pvwave$GPP_HBLR,df.REddy$GPP_DT,na.rm=T))
# plot(df.Pvwave$Reco,col=col.nt,ylim=ylim)
# points(df.Pvwave$Reco_HBLR,col=col.pvwave)
# points(df.REddy$Reco_DT,col=col.reddy_good)
### Reco
plot(df.Pvwave$Reco ~ df.Pvwave$Reco_HBLR,ylim=lim.reco,xlim=lim.reco,xlab="Daytime_Pvwave",ylab="Nighttime_MR",main="Reco")
curve(1*x,from=-1000,to=1000,col="red",add=T)
mod.pvwave <- lm(df.Pvwave$Reco ~ df.Pvwave$Reco_HBLR)
mtext(side=3,line=-1.8,paste0("R^2 = ",round(summary(mod.pvwave)$r.squared,2)))
title(main=sites[s],outer=T,cex=2)
plot(df.Pvwave$Reco ~ df.REddy$Reco_DT,ylim=lim.reco,xlim=lim.reco,xlab="Daytime_REddyProc",ylab="Nighttime_MR",main="Reco")
curve(1*x,from=-1000,to=1000,col="red",add=T)
mod.reddy <- lm(df.Pvwave$Reco ~ df.REddy$Reco_DT)
mtext(side=3,line=-1.8,paste0("R^2 = ",round(summary(mod.reddy)$r.squared,2)))
### GPP
plot(df.Pvwave$GPP_f ~ df.Pvwave$GPP_HBLR,ylim=lim.gpp,xlim=lim.gpp,xlab="Daytime_Pvwave",ylab="Nighttime_MR",main="GPP")
curve(1*x,from=-1000,to=1000,col="red",add=T)
mod.pvwave <- lm(df.Pvwave$GPP_f ~ df.Pvwave$GPP_HBLR)
mtext(side=3,line=-1.8,paste0("R^2 = ",round(summary(mod.pvwave)$r.squared,2)))
title(main=sites[s],outer=T,cex=1.5)
plot(df.Pvwave$GPP_f ~ df.REddy$GPP_DT,ylim=lim.gpp,xlim=lim.gpp,xlab="Daytime_REddyProc",ylab="Nighttime_MR",main="GPP")
curve(1*x,from=-1000,to=1000,col="red",add=T)
mod.reddy <- lm(df.Pvwave$GPP_f ~ df.REddy$GPP_DT)
mtext(side=3,line=-1.8,paste0("R^2 = ",round(summary(mod.reddy)$r.squared,2)))
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_Reco_GPP_DT_NT_scatterplots.pdf"),
width=9.5,height=9.5,pointsize=11)
# Reco and GPP timeseries
graphics.off()
par(mfrow=c(1,2),oma=c(1,1,2.5,0.5))
plot(df.Pvwave$Reco,col=col.nt,ylim=lim.reco,type="l",xlab="timestep",ylab="Reco",las=1)
points(df.Pvwave$Reco_HBLR,col=col.pvwave,type="l")
points(df.REddy$Reco_DT,col=col.reddy_good,type="l")
title(main=sites[s],outer=T,cex=1.5)
legend("topleft",legend=c("Nighttime","DT_Pvwave","DT_REddyProc"),col=c(col.nt,col.pvwave,col.reddy_good),
lty=1,bty="n",x.intersp=0.5,y.intersp=0.5,cex=0.7,seg.len=0.7)
plot(df.Pvwave$GPP_f,col=col.nt,ylim=lim.gpp,type="l",xlab="timestep",ylab="GPP",las=1)
points(df.Pvwave$GPP_HBLR,col=col.pvwave,type="l")
points(df.REddy$GPP_DT,col=col.reddy_good,type="l")
dev.copy2pdf(file=paste0(path,"Plots/",sites[s],"_Reco_GPP_DT_NT_timeseries.pdf"),
width=9.5,height=5,pointsize=11)
} # end site loop
#+++++++++++++++++++++++++++++
# Evaluation Monthly All
graphics.off()
par(mfrow=c(2,1),mar=c(4,4,1,1),mgp=c(1.8,0.4,0))
transp <- 0.5
r2_monthly <- round(summary(lm(Pvwave.mm.all$GPP_HBLR ~ REddy.mm.all$GPP_DT))$r.squared,3)
plot(Pvwave.mm.all$GPP_HBLR ~ REddy.mm.all$GPP_DT,xlab="GPP_REddyProc",ylab="GPP_Pvwave",
pch=20,col=alpha("black",transp),las=1,tcl=-0.2,xlim=c(0,20),ylim=c(0,20),cex=1.3)
curve(1*x,from=-10,to=100,col="red",add=T)
legend("bottomright","monthly",bty="n")
legend("topleft",paste0("R^2 = ",r2_monthly),bty="n")
#++++++++++++++++++++++++++++
# Evaluation Annual
r2_yearly <- round(summary(lm(Pvwave.yy.all$GPP_HBLR ~ REddy.yy.all$GPP_DT))$r.squared,3)
plot(Pvwave.yy.all$GPP_HBLR ~ REddy.yy.all$GPP_DT,xlab="GPP_REddyProc",ylab="GPP_Pvwave",
pch=20,col=alpha("black",transp),las=1,tcl=-0.2,xlim=c(0,11),ylim=c(0,11),cex=1.3)
curve(1*x,from=-10,to=100,col="red",add=T)
legend("bottomright","yearly",bty="n")
legend("topleft",paste0("R^2 = ",r2_yearly),bty="n")
dev.copy2pdf(file=paste0(path,"Plots/all_sites_monthly_yearly.pdf"),
width=7.5,height=8,pointsize=11)
## save as RData:
save(NT_vs_DT_REddy,DT_REddy_vs_pvwave,file=paste0(path,"Results/eval_metrics.RData")) # 1) evaluation metrics
save(Pvwave.mm.all,REddy.mm.all,file=paste0(path,"Results/all_sites_monthly.RData")) # 2) monthly aggregated results for all sites
save(Pvwave.yy.all,REddy.yy.all,file=paste0(path,"Results/all_sites_annual.RData")) # 3) annual aggregated results for all sites
Try the REddyProc package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.