Nothing
partitionNEEGL=function(
### Partitioning NEE fluxes into GP and Reco after daytime method of Lasslop et al. (2010)
ds ##<< dataset with all the specified input columns and full days in equidistant times
,NEEVar.s=paste0('NEE',SuffixDash.s,'_f') ##<< Variable of Net Ecosystem Exchange flux
,QFNEEVar.s=paste0('NEE',SuffixDash.s,'_fqc') ##<< Quality flag of variable
,QFNEEValue.n=0 ##<< Value of quality flag for _good_ (original) data
,NEESdVar.s=paste0('NEE',SuffixDash.s,'_fsd') ##<< Variable of standard deviation of net ecosystem fluxes
,TempVar.s=paste0('Tair_f') ##<< Filled air or soil temperature variable (degC)
,QFTempVar.s=paste0('Tair_fqc') ##<< Quality flag of filled temperature variable
,QFTempValue.n=0 ##<< Value of temperature quality flag for _good_ (original) data
,VPDVar.s=paste0('VPD_f') ##<< Filled Vapor Pressure Deficit - VPD - (hPa)
,QFVPDVar.s=paste0('VPD_fqc') ##<< Quality flag of filled VPD variable
,QFVPDValue.n=0 ##<< Value of VPD quality flag for _good_ (original) data
,RadVar.s='Rg' ##<< Unfilled (original) radiation variable
,PotRadVar.s="PotRad_NEW" ##<< Variable name of potential radiation (W/m2)
,Suffix.s = "" ##<< string inserted into column names before identifier (see \code{\link{sMDSGapFillAfterUstar}}).
,controlGLPart=partGLControl() ##<< further default parameters, see \code{\link{partGLControl}}
,isVerbose=TRUE ##<< set to FALSE to suppress output messages
,nRecInDay=48L ##<< number of records within one day (for half-hourly data its 48)
)
##description<<
## daytime-based partitioning of measured net ecosystem fluxes into gross primary production (GPP) and ecosystem respiration (Reco)
##author<<
## MM, TW
##references<<
## Lasslop G, Reichstein M, Papale D, et al. (2010) Separation of net ecosystem exchange into assimilation and respiration using
## a light response curve approach: critical issues and global evaluation. Global Change Biology, Volume 16, Issue 1, Pages 187-208
{
'Partitioning of measured net ecosystem fluxes into gross primary production (GPP) and ecosystem respiration (Reco) using Lasslop et al., 2010'
# Check if specified columns exist in sDATA or sTEMP and if numeric and plausible. Then apply quality flag
SuffixDash.s <- paste( (if(fCheckValString(Suffix.s)) "_" else ""), Suffix.s, sep="")
'Partitioning of measured net ecosystem fluxes into gross primary production (GPP) and ecosystem respiration (Reco) using Lasslop et al., 2010'
# Check if specified columns exist in sDATA or sTEMP and if numeric and plausible. Then apply quality flag
fCheckColNames(ds, c(NEEVar.s, QFNEEVar.s, TempVar.s, QFTempVar.s, RadVar.s, PotRadVar.s, NEESdVar.s), 'sGLFluxPartition')
fCheckColNum(ds, c(NEEVar.s, QFNEEVar.s, TempVar.s, QFTempVar.s, RadVar.s, PotRadVar.s, NEESdVar.s), 'sGLFluxPartition')
fCheckColPlausibility(ds, c(NEEVar.s, QFNEEVar.s, TempVar.s, QFTempVar.s, RadVar.s, PotRadVar.s), 'sGLFluxPartition')
Var.V.n <- fSetQF(ds, NEEVar.s, QFNEEVar.s, QFNEEValue.n, 'sGLFluxPartition')
if( isVerbose ) message('Start daytime flux partitioning for variable ', NEEVar.s, ' with temperature ', TempVar.s, '.')
##value<< data.frame with columns
## \item{Reco_DT_<suffix>}{predicted ecosystem respiraiton: mumol CO2/m2/second}
## \item{GPP_DT_<suffix>}{predicted gross primary production mumol CO2/m2/second}
RecoDTVar.s <- paste0('Reco_DT',SuffixDash.s)
GPPDTVar.s <- paste0('GPP_DT',SuffixDash.s)
RecoDTSdVar.s <- paste0(RecoDTVar.s,"_SD")
GPPDTSdVar.s <- paste0(GPPDTVar.s,"_SD")
dsAns <- data.frame(
FP_VARnight=rep(NA_real_,nrow(ds)) ##<< NEE filtered for nighttime records (others NA)
,FP_VARday=NA_real_ ##<< NEE filtered for daytime recores (others NA)
,NEW_FP_Temp=NA_real_ ##<< temperature after filtering for quality flag degree Celsius
,NEW_FP_VPD=NA_real_ ##<< vapour pressure deficit after filtering for quality flag, hPa
,FP_R_refNight=NA_real_ ##<< basal respiration estimated from LRC of daytime window (W/m2)
,FP_R_ref=NA_real_ ##<< basal respiration estimated from LRC of daytime window (W/m2)
,FP_E0=NA_real_ ##<< temperature sensitivity estimated from nighttime NEE window in Kelvin (degK)
,FP_alpha=NA_real_ ##<< canopy light utilization efficiency and represents the initial slope of the light-response curve (mumol s-1 W)
,FP_beta=NA_real_ ##<< maximum CO2 uptake rate of the canopy at light saturation (mumol m-2 s-1)
,FP_k=NA_real_ ##<< parameter controlling the VPD limitation of GPP
,FP_qc=NA_integer_ ##<< quality flag: 0: good parameter fit, 1: some parameters out of range, required refit, 2: next parameter estimate is more than two weeks away
,FP_dRecPar=NA_integer_ ##<< records until or after closest record that has a parameter estimate associated
)
## \item{}{Light response curve parameters \code{FP_X} are estimated for windows, and are reported with the first record of the window}
##end<<
# Filter night time values only
#! Note: Rg <= 4 congruent with Lasslop et al., 2010 to define Night for the calculation of E_0.n
# Should be unfilled (original) radiation variable, therefore dataframe set to sDATA only
isNight <- (ds[,RadVar.s] <= 4 & ds[[PotRadVar.s]] == 0)
dsAns$FP_VARnight <- ifelse(isNight, Var.V.n, NA)
attr(dsAns$FP_VARnight, 'varnames') <- paste(attr(Var.V.n, 'varnames'), '_night', sep='')
attr(dsAns$FP_VARnight, 'units') <- attr(Var.V.n, 'units')
# Filter day time values only
#! Note: Rg > 4 congruent with Lasslop et al., 2010 to define Day for the calculation of paremeters of Light Response Curve
# Should be unfilled (original) radiation variable, therefore dataframe set to sDATA only
isDay=(ds[,RadVar.s] > 4 & ds[[PotRadVar.s]] != 0)
dsAns$FP_VARday <- ifelse(isDay, Var.V.n, NA)
attr(dsAns$FP_VARday, 'varnames') <- paste(attr(Var.V.n, 'varnames'), '_day', sep='')
attr(dsAns$FP_VARday, 'units') <- attr(Var.V.n, 'units')
#! New code: Slightly different subset than PV-Wave due to time zone correction (avoids timezone offset between Rg and PotRad)
# Apply quality flag for temperature and VPD
# TODO: docu meteo filter, standard FALSE
dsAns$NEW_FP_Temp <- if( isTRUE(controlGLPart$isFilterMeteoQualityFlag) ) fSetQF(ds, TempVar.s, QFTempVar.s, QFTempValue.n, 'partitionNEEGL') else ds[[TempVar.s]]
dsAns$NEW_FP_VPD <- if( isTRUE(controlGLPart$isFilterMeteoQualityFlag) ) fSetQF(ds, VPDVar.s, QFVPDVar.s, QFVPDValue.n, 'partitionNEEGL') else ds[[VPDVar.s]]
#Estimate Parameters of light response curve: R_ref, alpha, beta and k according to Table A1 (Lasslop et al., 2010)
# save(ds, file="tmp/dsTestPartitioningLasslop10.RData")
# extract the relevant columns in df with defined names (instead of passing many variables)
dsR <- data.frame(
sDateTime=ds[[1]] # not used, but usually first column is a dateTime is kept for aiding debug
,NEE=Var.V.n
,sdNEE=ds[[NEESdVar.s]]
, Temp=dsAns$NEW_FP_Temp
, VPD=dsAns$NEW_FP_VPD
, Rg=ds[[RadVar.s]]
, isDay=isDay
, isNight=isNight
)
##seealso<< \code{\link{partGLFitLRCWindows}}
resParms <- partGLFitLRCWindows( dsR
, nRecInDay=nRecInDay
, controlGLPart=controlGLPart
)
if( FALSE ){
resLRCOld <- tmp <- .depr.partGLFitLRCWindows( dsR
, nRecInDay=nRecInDay
, controlGLPart.l=controlGLPart
)
resParms <- resLRCOld
}
# append parameter fits to the central record of day window
#iGood <- which(resLRC$summary$parms_out_range == 0L)
# resLRC provides parameters only for a subset of rows. For each row in the original data.frame
# The parameter estimates are associated for interpolation of predicted fluxes
# either to the central record of the window,
# or the record with the time corresponding to the mean of all valid records in the window
# default is isAssociateParmsToMeanOfValids=TRUE (double check partGLControl argument)
colNameAssoc <- if( isTRUE(controlGLPart$isAssociateParmsToMeanOfValids) ) "iMeanRec" else "iCentralRec"
# for the output, always report at central record
dsAns[resParms$summary$iCentralRec,c("FP_R_refNight","FP_E0","FP_R_ref","FP_alpha","FP_beta","FP_k","FP_qc")] <-
resParms$summary[,c("R_ref_night","E_0","R_ref","a","b","k","parms_out_range")]
matchFP_qc <- NA_integer_; matchFP_qc[resParms$summary[[colNameAssoc]] ] <- resParms$summary$parms_out_range # here maybe indexed by meanRec
#
##seealso<< \code{\link{partGLInterpolateFluxes}}
dsAnsFluxes <- partGLInterpolateFluxes( ds[,RadVar.s]
#, dsAns$NEW_FP_VPD, dsAns$NEW_FP_Temp
, ds[[VPDVar.s]], ds[[TempVar.s]] # do prediction also using gap-Filled values
, resParms
, controlGLPart=controlGLPart
)
# set quality flag to 2 where next parameter estimate is more than 14 days away
dsAns$FP_dRecPar <- dsAnsFluxes$dRecNextEstimate
dDaysPar <- round(dsAns$FP_dRecPar / nRecInDay)
dsAns$FP_qc <- matchFP_qc[ 1:nrow(dsAns) + dsAns$FP_dRecPar ] # associate quality flag of parameter estimate to each record
dsAns$FP_qc[ dDaysPar > 14] <- 2L # set quality flag to 2 for records where next estimate is more than 14 days away
#dsAns[is.finite(dsAns$FP_beta),]
#
dsAns[[RecoDTVar.s]] <- dsAnsFluxes$Reco
attr(dsAns[[RecoDTVar.s]], 'varnames') <- RecoDTVar.s
attr(dsAns[[RecoDTVar.s]], 'units') <- attr(Var.V.n, 'units')
dsAns[[GPPDTVar.s]] <- dsAnsFluxes$GPP
attr(dsAns[[GPPDTVar.s]], 'varnames') <- GPPDTVar.s
attr(dsAns[[GPPDTVar.s]], 'units') <- attr(Var.V.n, 'units')
if( controlGLPart$isSdPredComputed ){
dsAns[[RecoDTSdVar.s]] <- dsAnsFluxes$sdReco
attr(dsAns[[RecoDTSdVar.s]], 'varnames') <- RecoDTSdVar.s
attr(dsAns[[RecoDTSdVar.s]], 'units') <- attr(Var.V.n, 'units')
dsAns[[GPPDTSdVar.s]] <- dsAnsFluxes$sdGPP
attr(dsAns[[GPPDTSdVar.s]], 'varnames') <- GPPDTSdVar.s
attr(dsAns[[GPPDTSdVar.s]], 'units') <- attr(Var.V.n, 'units')
}
#sTEMP$GPP_DT_fqc <<- cbind(sDATA,sTEMP)[,QFFluxVar.s]
#! New code: MDS gap filling information are not copied from NEE_fmet and NEE_fwin to GPP_fmet and GPP_fwin
# (since not known within this pure partitioning function)
return(dsAns)
}
partGLControl <- function(
### Default list of parameters for Lasslop 2010 daytime flux partitioning
LRCFitConvergenceTolerance=1e-3 ##<< convergence criterion for LRC fit.
## If relative improvement of reducing residual sum of squares between predictions and
## observations is less than this criterion, assume convergence.
## Decrease to get more precise parameter estimates, Increase for speedup.
,nBootUncertainty=30L ##<< number of bootstrap samples for estimating uncertainty.
## Set to zero to derive uncertainty from curvature of a single fit
,minNRecInDayWindow = 10L ##<< Minimum number of data points for regression
,isAssociateParmsToMeanOfValids=TRUE ##<< set to FALSE to associate parameters to
## the first record of the window for interpolation
## instead of mean across valid records inside a window
,isLasslopPriorsApplied=FALSE ##<< set to TRUE to apply strong fixed priors on LRC fitting.
## Returned parameter estimates claimed valid for some case where not enough data was available
,isSdPredComputed=TRUE ##<< set to FALSE to avoid computing standard errors
## of Reco and GPP for small performance increase
,isFilterMeteoQualityFlag=FALSE ##<< set to TRUE to use only records where quality flag
## of meteo drivers (Radation, Temperatrue, VPD) is zero, i.e. non-gapfilled for parameter estimation.
## For prediction, the gap-filled value is used always, to produce predictions also for gaps.
,isBoundLowerNEEUncertainty=TRUE ##<< set to FALSE to avoid adjustment of very low uncertainties before
## day-Time fitting that avoids the high leverage those records with unreasonable low uncertainty.
,smoothTempSensEstimateAcrossTime=TRUE ##<< set to FALSE to use independent estimates of temperature
## sensitivity on each windows instead of a vector of E0 that is smoothed over time
){
##author<< TW
##seealso<< \code{\link{partitionNEEGL}}
##description<<
## For highest compatibility to Lasslop10 use
## \code{nBootUncertainty=0L, isAssociateParmsToMeanOfValids=FALSE
## , isLasslopPriorsApplied=TRUE, isBoundLowerNEEUncertainty=FALSE
## , smoothTempSensEstimateAcrossTime=FALSE
## }
ctrl <- list(
LRCFitConvergenceTolerance=LRCFitConvergenceTolerance
,nBootUncertainty=nBootUncertainty
,minNRecInDayWindow=minNRecInDayWindow
,isAssociateParmsToMeanOfValids=isAssociateParmsToMeanOfValids
,isLasslopPriorsApplied=isLasslopPriorsApplied
,isSdPredComputed=isSdPredComputed
,isFilterMeteoQualityFlag=isFilterMeteoQualityFlag
,isBoundLowerNEEUncertainty=isBoundLowerNEEUncertainty
,smoothTempSensEstimateAcrossTime=smoothTempSensEstimateAcrossTime
)
#display warning message for the following variables that we advise not to be changed
#if (corrCheck != 0.5) warning("WARNING: parameter corrCheck set to non default value!")
##value<< list with entries of given arguments.
ctrl
}
attr(partGLControl,"ex") <- function(){
partGLControl(nBootUncertainty=40L)
}
partGLFitLRCWindows=function(
### Estimate temperature sensitivity and parameters of Rectangular Hyperbolic Light Response Curve function (a,b,R_ref, k) for successive periods
ds ##<< data.frame with numeric columns NEE, sdNEE, Temp (degC), VPD, Rg, and logical columns isNight and isDay
,WinSizeDays.i=4L ##<< Window size in days for daytime fits
,WinSizeNight.i=3L*WinSizeDays.i ##<< Window size in days for nighttime fits
,DayStep.i=floor(WinSizeDays.i / 2L)##<< step in days for shifting the windows
,isVerbose=TRUE ##<< set to FALSE to suppress messages
,nRecInDay=48L ##<< number of records within one day (for half-hourly data its 48)
,winExtendSizes = WinSizeNight.i*c(2L,4L) ##<< successively increased nighttime windows, to obtain a night-time fit
,controlGLPart=partGLControl() ##<< list of further default parameters
){
##seealso<< \code{\link{partGLFitNightTempSensOneWindow}}
if( isVerbose ) message(" Estimating temperature sensitivity from night time NEE ", appendLF = FALSE)
resNight <- applyWindows(ds, partGLFitNightTempSensOneWindow, prevRes=list(prevE0=NA)
, winSizeRefInDays = WinSizeDays.i
, winSizeInDays=WinSizeNight.i
,isVerbose=isVerbose
,nRecInDay=nRecInDay
#,isNAAllowed=TRUE
#
,controlGLPart=controlGLPart
)
iNoFit <- which( is.na(resNight$summary$E0) )
iExtend <- 1
while( length(iNoFit) && (iExtend <= length(winExtendSizes)) ){
if( isVerbose ) message(" increase window size to ",winExtendSizes[iExtend], appendLF = FALSE)
resNightExtend <- applyWindows(ds, partGLFitNightTempSensOneWindow, prevRes=list(prevE0=NA)
, winSizeRefInDays = WinSizeDays.i
,winSizeInDays=winExtendSizes[iExtend]
,isVerbose=isVerbose
,nRecInDay=nRecInDay
#,isNAAllowed= (iExtend < length(winExtendSizes))
#
,controlGLPart=controlGLPart
)
resNight$resOptList[iNoFit] <- resNightExtend$resOptList[iNoFit]
resNight$summary[iNoFit,] <- resNightExtend$summary[iNoFit,]
#all(resNight$summary$iCentralRec == resNight24$summary$iCentralRec)
iNoFit <- which( is.na(resNight$summary$E0) )
iExtend <- iExtend + 1L
}
# there might still some windows, where RRefNight has not been determined, but its needed as starting value for LRC estimate
resNight$summary$RRefFit <- resNight$summary$RRef # store original estimate
resNight$summary$RRef <- fillNAForward(resNight$summary$RRef, firstValue=
resNight$summary$RRef[which(is.finite(resNight$summary$RRef))[1] ])
#
##seealso<< \code{\link{partGLSmoothTempSens}}
# remember E0 and sdE0 before overidden by smoothing
resNight$summary$E0Fit <- resNight$summary$E0
resNight$summary$sdE0Fit <- resNight$summary$sdE0
E0Smooth <- if( isTRUE(controlGLPart$smoothTempSensEstimateAcrossTime) ){
if( isVerbose ) message(" Smoothing temperature sensitivity estimates")
partGLSmoothTempSens( resNight$summary )
}else {
E0Smooth <- resNight$summary
iNonFiniteE0 <- which(!is.finite(E0Smooth$E0))
E0Smooth$sdE0[iNonFiniteE0] <- quantile(E0Smooth$sdE0, 0.9, na.rm=TRUE) # set uncertainty to the 90% quantile of the distribution of uncertainties
E0Smooth$E0 <- fillNAForward(E0Smooth$E0) # fill NA with value from previous window
E0Smooth
}
# now all E0 and sdE0 are defined
#
##seealso<< \code{\link{partGLFitNightRespRefOneWindow}}
if( isVerbose ) message(" Estimating respiration at reference temperature for smoothed temperature sensitivity from night time NEE ", appendLF = FALSE)
resRef15 <- applyWindows(ds, partGLFitNightRespRefOneWindow
, winSizeRefInDays = WinSizeDays.i
, winSizeInDays=WinSizeDays.i
,isVerbose=isVerbose
,nRecInDay=nRecInDay
#
,E0Win = E0Smooth
,controlGLPart=controlGLPart
)
E0Smooth$RRef <- fillNAForward( resRef15$summary$RRef ) # may contain NA if not enough night-time records
# special case NA on first record fill backward
iFirstFinite <- min(which(is.finite(E0Smooth$RRef)))
if(iFirstFinite > 1) E0Smooth$RRef[1:(iFirstFinite-1)] <- E0Smooth$RRef[iFirstFinite]
#
##seealso<< \code{\link{partGLFitLRCOneWindow}}
if( isVerbose ) message(" Estimating light response curve parameters from day time NEE ", appendLF = FALSE)
resLRC <- applyWindows(ds, partGLFitLRCOneWindow, prevRes<-list(lastGoodParameters=rep(NA_real_, 5L))
, winSizeInDays=4L
,isVerbose=isVerbose
,nRecInDay=nRecInDay
#
,E0Win = E0Smooth
,controlGLPart=controlGLPart
)
resParms <- resLRC
resParms$summary$E_0 <- E0Smooth$E0
resParms$summary$E_0_sd <- E0Smooth$sdE0
resParms$summary$R_ref_night <- E0Smooth$RRef
# summary$iMeanRec yet based on window instead of entire time, need to add beginning of window
resParms$summary$iMeanRec <- resParms$summary$iRecStart-1L + resParms$summary$iMeanRec
# omit records where NULL was returned
iWinNoFit <- which( is.na(resParms$summary$parms_out_range) )
if( length(iWinNoFit) ){
resParms$summary <- resParms$summary[-iWinNoFit, ]
resParms$resOptList <- resParms$resOptList[-iWinNoFit]
}
resParms
}
.tmp.f <- function(){
plot( E0 ~ dayStart, E0Smooth)
}
fillNAForward <- function(
### replace NA by value of previous record
x ##<< numeric vector to fill NAs
, firstValue=median(x,na.rm=TRUE) ##<< value to be used for NA at the beginning of x
){
iMissing <- which(!is.finite(x))
if( length(iMissing) && (iMissing[1] == 1L)){
# set first vluae
x[1L] <- firstValue
iMissing <- iMissing[-1]
}
if(length(iMissing)){
for( i in iMissing ){
# set to value from previous window
x[i] <- x[i-1L]
}
}
return(x)
}
applyWindows <- function(
### apply a function to several windows of a data.frame
ds ##<< dataframe to iterate
,FUN ##<< function to apply to subsets of the data.frame
## taking a subset of the data.frame as first argument
## the second: a one-row data.frame with window information (iWindow, dayStart, dayEnd, iRecStart, iRecEnd, iCentralRec)
## the third a list that can transport values from previous fits, i.e. the third item of its return value
## with first item a list (objects that cannot be rowbound)
## second item a single row data.frame (summary that can be row-bound)
## third item a list that is will be provided to the call of this function on the next window as the second argument
,prevRes=list() ##<< initial values for the list that is carried between windows
,winSizeInDays=winSizeRefInDays ##<< Window size in days
,winSizeRefInDays=4L ##<< Window size in days for reference window (e.g. day-Window for night time)
,strideInDays=floor(winSizeRefInDays/2L) ##<< step in days for shifting the window, for alligning usually a factor of winSizeRef
,isVerbose=TRUE ##<< set to FALSE to suppress messages
,nRecInDay=48L ##<< number of records within one day (for half-hourly data its 48)
,... ##<< further arguments to FUN
){
##details<<
## Assumes equidistant rows with nRecInDay records forming one day and reporting full days
## i.e. all of the nRecInDay records are in the first day.
##details<<
## In order to have a common reference winSizeRefInDays is given so that results by a different window
## size correspond to each window of shifting a window of winSizeRefInDays
## Each window is anchord so that the center equals the center of the reference window.
## This becomes important when selecting records at the edges of the series.
nRec <- nrow(ds)
nDay <- ceiling( nRec / nRecInDay)
nDayLastWindow <- nDay - (winSizeRefInDays/2) # center of the reference window still in records
#iDayOfRec <- ((c(1:nRec)-1L) %/% nRecInDay)+1L # specifying the day for each record assuming equidistand records
startDaysRef <- seq(1, nDayLastWindow, strideInDays) # starting days for each reference window
iCentralRec <- ((startDaysRef-1L)+winSizeRefInDays/2)*nRecInDay+1L # assuming equidistant records
nWindow <- length(startDaysRef)
# precomputing the starting and end records for all periods in vectorized way
dayStart0 <- as.integer(startDaysRef + winSizeRefInDays/2 - winSizeInDays/2) # may become negative, negative needed for computation of dayEnd
dayStart <- pmax(1L, dayStart0)
dayEnd <- pmin(nDay, dayStart0-1L+winSizeInDays)
iRecStart0 <- iCentralRec -winSizeInDays/2*nRecInDay # may become negative, negative needed for computation of iRecEnd
iRecStart <- pmax(1L, iRecStart0 )
iRecEnd <- pmin(nRec, iCentralRec-1L +winSizeInDays/2*nRecInDay )
# central record in each window, only last record, refers to reference window of size winSizeRefInDays
dsRec <- data.frame(
iWindow=1:nWindow
,dayStart = dayStart
,dayEnd=dayEnd
,iRecStart=iRecStart
,iRecEnd=iRecEnd
,iCentralRec=iCentralRec
)
res1 <- vector("list",nWindow)
res2List <- vector("list",nWindow) # each will hold a data.frame to be row-bound afterwards (dont know the columns yet)
for( iWindow in 1:nWindow){
if( isVerbose ) message(",",startDaysRef[iWindow], appendLF = FALSE)
startRec <- iRecStart[iWindow]
endRec <- iRecEnd[iWindow]
# range( which(iDayOfRec >= startDay & iDayOfRec <= endDay)) # slower but check for startRec and endRec
dsWin <- ds[startRec:endRec,]
resFun <- FUN(dsWin, dsRec[iWindow,], prevRes, ...)
if( length(resFun) ){
res1[[iWindow]] <- resFun[[1]]
res2List[[iWindow]] <- resFun[[2]]
prevRes <- resFun[[3]]
}
}
if( isVerbose ) message("") # LineFeed
#tmp <- res2List
#res2List <- tmp
# rbind res2List, but keep NULL results, therefore add dummy dataframe, and delete dummy column afterward
res2List[ sapply(res2List, is.null) ] <- list(data.frame(..DUMMY=1))
res2 <- rbind.fill(res2List) # maybe later upgrade to dplyr
res2$..DUMMY <- NULL
dsSummary <- cbind(dsRec, res2)
##value<< a list with components that correspond to the three result components of FUN
return(list(
resOptList=res1 ##<< a list with restults of length nrow(ds)
,summary=dsSummary ##<< a data.frame with nrow(ds) rows, also including window information (iWindow, dayStart, dayEnd, iRecStart, iRecEnd, iCentralRec)
, prevRes=prevRes) ##<< last result of the values carried over between windows
)
}
partGLFitNightTempSensOneWindow=function(
### Estimate parameters of the Rectangular Hyperbolic Light Response Curve function (a,b,R_ref, k) for successive periods
dss ##<< data.frame with numeric columns NEE, sdNEE, Temp (degC), VPD, Rg, and logical columns isNight and isDay
,winInfo ##<< one-row data.frame with window information, including iWindow
,prevRes ##<< component prevRes from previous result, here with item prevE0
,isVerbose=TRUE ##<< set to FALSE to suppress messages
,nRecInDay=48L ##<< number of records within one day (for half-hourly data its 48)
,controlGLPart=partGLControl() ##<< list of further default parameters
){
##author<<
## TW
##seealso<< \code{\link{partitionNEEGL}}
##description<<
## Estimation of respiration at reference temperature (R_Ref) and temperature (E_0) for one window.
isValid <- isValidNightRecord(dss)
# check that there are enough night and enough day-values for fitting, else continue with next window
if( sum(isValid) < controlGLPart$minNRecInDayWindow ) return(NULL)
dssNight <- dss[isValid,]
# tmp <- dss[!is.na(dss$isNight) & dss$isNight & !is.na(dss$NEE), ]; plot(NEE ~ Temp, tmp)
# points(NEE ~ Temp, dssNight, col="blue" )
##seealso<< \code{\link{partGLEstimateTempSensInBoundsE0Only}}
#TODO here only E0 is of interest, do not need to reestimate RRef
#resNightFitOld <- partGLEstimateTempSensInBounds(dssNight$NEE, fConvertCtoK(dssNight$Temp)
# , prevE0=prevRes$prevE0)
resNightFit <- partGLEstimateTempSensInBoundsE0Only(dssNight$NEE, fConvertCtoK(dssNight$Temp)
, prevE0=prevRes$prevE0)
return(list(
resNightFit
,data.frame(E0=resNightFit$E0, sdE0=resNightFit$sdE0, TRefFit=resNightFit$TRefFit, RRefFit=resNightFit$RRefFit)
,list(prevE0=if(is.finite(resNightFit$E0)) resNightFit$E0 else prevRes$prevE0)
))
}
.tmp.f <- function(){
plot( dss$NEE ~ dss$sDateTime)
plot( dssNight$NEE ~ dssNight$sDateTime)
}
isValidNightRecord <- function(
### compute logical vector of each rows in ds is its a valid night record
ds ##<< data.frame with columns isNight, NEE, Temp (degC)
){
isValid <- !is.na(ds$isNight) & ds$isNight & !is.na(ds$NEE) & is.finite(ds$Temp)
##details<<
## For robustness, data is trimmed to conditions at temperature > 1 degC
## but only timmed if there are more at least 12 records left
isFreezing <- ds$Temp[isValid] <= -1
if( sum(is.finite(isFreezing) & !isFreezing) >= 12L ) isValid[isValid][isFreezing|is.na(isFreezing)] <- FALSE
##value<< a logical vector of length nrow(ds)
return(isValid)
}
partGLEstimateTempSensInBoundsE0Only <- function(
### Estimate temperature sensitivity E_0 and R_ref of ecosystem respiration, and apply bounds or previous estimate
REco ##<< numeric vector: night time NEE, i.e. ecosytem respiration
,temperatureKelvin ##<< numeric vector: temperature in Kelvin of same length as REco
,prevE0 = NA ##<< numeric scalar: the previous guess of Temperature Sensitivity
){
##author<< MM, TW
##seealso<< \code{\link{partGLFitLRCWindows}}
#twutz: using nls to avoid additional package dependency
#resFitLM <- NLS.L <- nlsLM(formula=R_eco ~ fLloydTaylor(R_ref, E_0, Temp, T_ref.n=273.15+15), algorithm='default', trace=FALSE,
# data=as.data.frame(cbind(R_eco=REco.V.n,Temp=temperatureKelvin.V.n)), start=list(R_ref=mean(REco.V.n,na.rm=TRUE),E_0=100)
# ,control=nls.lm.control(maxiter = 20))
##details<<
## Basal respiration is reported for temperature of 15 degree Celsius. However during the fit
## a reference temperature of the median of the dataset is used. This is done to avoid
## strong correlations between estimated parameters E0 and R_ref, that occure if reference temperature
## is outside the center of the data.
TRefFit <- median(temperatureKelvin, na.rm=TRUE) # formerly 273.15+15
resFit <- try(
nls(formula=R_eco ~ fLloydTaylor(R_ref, E_0, Temp, T_ref.n=TRefFit), algorithm='default', trace=FALSE,
data=as.data.frame(cbind(R_eco=REco,Temp=temperatureKelvin))
, start=list(R_ref= mean(REco,na.rm=TRUE)
,E_0=as.vector({if(is.finite(prevE0)) prevE0 else 100}))
,control=nls.control(maxiter = 20L)
)
, silent=TRUE)
#plot( REco.V.n ~ I(temperatureKelvin.V.n-273.15) )
#plot( REcoCorrected ~ I(temperatureKelvin.V.n-273.15)[isNotFreezing] )
if( inherits(resFit, "try-error")){
#stop("debug partGLEstimateTempSensInBounds")
#plot( REco.V.n ~ temperatureKelvin.V.n )
E0 <- NA
sdE0 <- NA
RRefFit <- NA
} else {
E0 <- coef(resFit)['E_0']
sdE0 <- coef(summary(resFit))['E_0',2]
RRefFit <- coef(resFit)['R_ref']
}
# resFit$convInfo$isConv
##details<<
## If E_0 is out of bounds [50,400] then report E0 as NA
if( is.na(E0) || (E0 < 50) || (E0 > 400)){
E0 <- NA
sdE0 <- NA
RRefFit <- NA
}
##value<< list with entries
return(list(
E0=E0 ##<< numeric scalar of estimated temperature sensitivty E0 bounded to [50,400]
,sdE0=sdE0 ##<< numeric scalar of standard deviation of E0
,TRefFit=TRefFit ##<< numeric scalar reference temperature used in the E0 fit
,RRefFit=RRefFit ##<< numeric scalar respiration at TRefFit
))
#
# refit R_Ref with bounded E0 for 15 degC, instead of calling fLoydAndTaylor do a simple regression
# starting value from forward model
# RRef15 <- fLloydTaylor( RRefFit, E_0Bounded.V.n, TRef15, T_ref.n=TRefFit)
# resFit15 <- nls(formula=R_eco ~ fLloydTaylor(R_ref, E_0, Temp, T_ref.n=TRef15), algorithm='default', trace=FALSE,
# data=as.data.frame(cbind(R_eco=REcoFitting,Temp=temperatureKelvin.V.n[isNotFreezing], E_0=E_0Bounded.V.n))
# , start=list(R_ref=RRef15)
# ,control=nls.control(maxiter = 20L)
# )
# (R_ref_ <- coef(resFit15)[1])
}
.tmp.f <- function(){
# from recover at fitting E0
temp <- temperatureKelvin - 273.15
plot( REco ~ temp )
lines( fLloydTaylor(RRefFit, E0, temperatureKelvin, T_ref.n=TRefFit) ~ temp)
}
partGLSmoothTempSens <- function(
### Smoothes time development of E0
E0Win ##<< data.frame with columns E0 and sdE0, RRefFit, and TRefFit with one row for each window
){
#return(E0Win)
#E0Win$E0[1] <- NA
E0Win$E0[c(FALSE,(diff(E0Win$E0) == 0))] <- NA # TODO return NA in the first place where the previous window was used
# mlegp does not work for long time series (Trevors bug report with Havard data)
# hence, smooth one year at a time
# testing: E0Win <- do.call( rbind, lapply(0:4,function(i){tmp<-E0Win; tmp$dayStart <- tmp$dayStart+i*365; tmp}))
E0Win$year <- ceiling(E0Win$dayStart/365)
yr <- E0Win$year[1]
resWin <- lapply( unique(E0Win$year), function(yr){
E0WinYr <- E0Win[ E0Win$year == yr, ,drop=FALSE]
isFiniteE0 <- is.finite(E0WinYr$E0)
E0WinFinite <- E0WinYr[ isFiniteE0, ]
output <- capture.output(
gpFit <- mlegp(X=E0WinFinite$iCentralRec, Z=E0WinFinite$E0, nugget = E0WinFinite$sdE0^2)
#gpFit <- mlegp(X=E0WinFinite$iCentralRec, Z=E0WinFinite$E0, nugget = (E0WinFinite$sdE0*2)^2, nugget.known=1L)
)
pred1 <- predict(gpFit, matrix(E0WinYr$iCentralRec, ncol=1), se.fit=TRUE)
nuggetNewObs <- quantile(gpFit$nugget, 0.9)
nugget <- rep(nuggetNewObs, nrow(E0WinYr))
nugget[isFiniteE0] <- gpFit$nugget
E0WinYr$E0 <- pred1$fit
E0WinYr$sdE0 <- pred1$se.fit + sqrt(nugget)
E0WinYr
})
##value<< dataframe E0Win with updated columns E0 and sdE0
#recover()
ans <- do.call(rbind,resWin)
}
.tmp.f <- function(){
plot( E0WinFinite$E0 ~ E0WinFinite$iCentralRec )
#arrows( E0WinFinite$iCentralRec, E0WinFinite$E0-1.96*E0WinFinite$sdE0, y1=E0WinFinite$E0+1.96*E0WinFinite$sdE0, length=0, col="grey" )
arrows( E0WinFinite$iCentralRec, E0WinFinite$E0-E0WinFinite$sdE0, y1=E0WinFinite$E0+E0WinFinite$sdE0, length=0, col="grey" )
#
E0Win$day <- (E0Win$iCentralRec-1) / 48 +1
E0WinFinite$day <- (E0WinFinite$iCentralRec-1) / 48 +1
plot( E0WinFinite$E0 ~ E0WinFinite$day )
plot( E0Win$E0Fit ~ E0Win$day )
points( E0Win$E0 ~ E0Win$day, col="blue", type="b", lty="dotted" )
#points( E0Win$E0 ~ E0Win$day, col="blue" )
#arrows( E0Win$day, E0Win$E0Fit, y1=E0Win$E0, col="grey", length=0.1)
lines( I(E0Win$E0 + 1.06*E0Win$sdE0) ~ E0Win$day, col="lightblue" )
lines( I(E0Win$E0 - 1.06*E0Win$sdE0) ~ E0Win$day, col="lightblue" )
#
E0Win$E0 <- E0Win$E0Fit
E0Win$sdE0 <- E0Win$sdE0Fit
E0Win$E0[c(FALSE,(diff(E0Win$E0) == 0))] <- NA # TODO return NA in the first place
#
E0Win$isFiniteE0 <- isFiniteE0
E0Win[ E0Win$E0 < 80,]
E0Win[ 65:75,]
subset(E0WinFinite, iWindow %in% 65:75)
E0Win$sdE0 / E0Win$E0
}
partGLFitNightRespRefOneWindow=function(
### Estimate Reference temperature from nicht-time and given temperature sensitivity E0
dss ##<< data.frame with numeric columns NEE, isNight, Temp, Rg
,winInfo ##<< one-row data.frame with window information, including iWindow
,prevRes=list() ##<< component prevRes from previous result, here not used.
,E0Win ##<< data.frame with columns E0 and sdE0, RRefFit, and TRefFit with one row for each window
,isVerbose=TRUE ##<< set to FALSE to suppress messages
,nRecInDay=48L ##<< number of records within one day (for half-hourly data its 48)
,controlGLPart=partGLControl() ##<< list of further default parameters
){
##author<<
## TW
##seealso<< \code{\link{partitionNEEGL}}
##description<<
## Estimation of respiration at reference temperature (R_Ref) and temperature (E_0) for one window.
isValid <- isValidNightRecord(dss)
# check that there are enough night and enough day-values for fitting, else continue with next window
if( sum(isValid) < controlGLPart$minNRecInDayWindow ) return(NULL)
dssNight <- dss[isValid,]
# tmp <- dss[!is.na(dss$isNight) & dss$isNight & !is.na(dss$NEE), ]; plot(NEE ~ Temp, tmp)
# points(NEE ~ Temp, dssNight, col="blue" )
##seealso<< \code{\link{partGLEstimateTempSensInBoundsE0Only}}
REco <- dssNight$NEE
E0 <- E0Win$E0[winInfo$iWindow]
TRef15 <- 273.15+15 # 15degC in Kelvin
R_ref <- if( length(REco) >= 3L ){
temperatureKelvin <- 273.15+dssNight$Temp
T_0.n=227.13 ##<< Regression temperature as fitted by LloydTaylor (1994) in Kelvin (degK)
TFacLloydTaylor <- exp(E0 * ( 1/(TRef15-T_0.n) - 1/(temperatureKelvin-T_0.n) ) )
lm15 <- lm(REco ~ TFacLloydTaylor -1)
coef(lm15)
} else
fLloydTaylor( E0Win$RRefFit[winInfo$iWindow], E0, TRef15, T_ref.n=E0Win$TRefFit[winInfo$iWindow])
R_refBounded <- max(0, R_ref)
##value<< list with entries
return(list(
NULL
,data.frame(RRef=R_refBounded)
,prevRes
))
}
partGLFitLRCOneWindow=function(
### Estimate parameters of the Rectangular Hyperbolic Light Response Curve function (a,b,R_ref, k) for successive periods
ds ##<< data.frame with numeric columns NEE, sdNEE, Temp (degC), VPD, Rg, and logical columns isNight and isDay
,winInfo ##<< one-row data.frame with window information, including iWindow
,prevRes ##<< component prevRes from previous result, here with item prevE0
,E0Win ##<< data.frame with columns E0, sdE0, RRef from nighttime, one row for each window
,controlGLPart=partGLControl() ##<< list of further default parameters
){
##author<< TW
##seealso<< \code{\link{partitionNEEGL}}
##description<<
## Estimation as in Lasslop et al., 2010 for successive periods, i.e. windows.
#requiredCols <- c("NEE", "sdNEE", "Temp", "VPD", "Rg", "isNight", "isDay")
#iMissing <- which( is.na(match( requiredCols, names(ds) )))
#if( length(iMissing) ) stop("missing columns: ",paste0(requiredCols[iMissing],collapse=","))
isValidDayRec <- !is.na(ds$isDay) & ds$isDay & !is.na(ds$NEE) & !is.na(ds$Temp) & !is.na(ds$VPD)
# check that there are enough night and enough day-values for fitting, else continue with next window
if( sum(isValidDayRec) < controlGLPart$minNRecInDayWindow ) return( NULL );
dsDay <- ds[isValidDayRec,]
##details<<
## Each window estimate is associated with a time or equivalently with a record.
## The first record, i.e. row number, of the day-window is reported.
## Moreover, the mean of all valid records numbers in the daytime window is reported for interpolation.
iMeanRecInDayWindow <- as.integer(round(mean(which(isValidDayRec))))
#TODO firstRecInDayWindow.i <- which(SubsetDayPeriod.b)[1] # the rownumber of the first record inside the day window
##seealso<< \code{\link{partGLEstimateTempSensInBoundsE0Only}}
E0 <- E0Win$E0[winInfo$iWindow]
# if no temperature - respiration relationship could be found, indicate no-fit
if( is.na(E0) ) return(NULL)
# if( DayStart.i > 72 ) recover()
sdE0 <- E0Win$sdE0[winInfo$iWindow]
RRefNight <- E0Win$RRef[winInfo$iWindow]
# RRefNight might not be determined due to lack of night-time data, still set to reasonable starting value
# better fill forward RRefNight, so that last estimate is used
# if( is.na(RRefNight) ) RRefNight<-mean(E0Win$RRef,na.rm = TRUE)
#
##seealso<< \code{\link{partGLFitLRC}}
resOpt <- resOpt0 <- partGLFitLRC(dsDay, E0=E0, sdE0=sdE0, RRefNight=RRefNight
, controlGLPart=controlGLPart, lastGoodParameters=prevRes$lastGoodParameters)
#
if( !is.finite(resOpt$opt.parms.V[1]) ) return(NULL)
# check the Beta bounds that depend on uncertainty, set to NA fit
sdParms <- resOpt$opt.parms.V; sdParms[] <- NA
sdParms[resOpt$iOpt] <- sqrt(diag(resOpt$covParms)[resOpt$iOpt])
if(isTRUE(as.vector( (resOpt$opt.parms.V[2] > 100) && (sdParms[2] >= resOpt$opt.parms.V[2]) ))){
return(NULL)
}
# check that R_ref estimated from daytime is not both:
# larger than twice the estimate from nighttime and more than 0.7 in absolute terms
# else this indicates a bad fit
# this is additional to Table A1 in Lasslop 2010
if( (resOpt$opt.parms.V[4L] > 2*RRefNight) &&
((resOpt$opt.parms.V[4L]-RRefNight) > 0.7)
){
return(NULL)
}
#
#recover()
prevRes$lastGoodParameters <- resOpt$opt.parms.V
# record valid fits results
ans <- list(
resOpt=resOpt
,summary = data.frame(
nValidRec=nrow(dsDay)
,iMeanRec=iMeanRecInDayWindow
,E_0=E0, E_0_sd=sdE0
,R_ref=resOpt$opt.parms.V[4], R_ref_SD=sdParms[4]
,a=resOpt$opt.parms.V[3], a_SD=sdParms[3]
,b=resOpt$opt.parms.V[2], b_SD=sdParms[2]
,k=resOpt$opt.parms.V[1], k_SD=sdParms[1]
,parms_out_range=as.integer(!identical(resOpt$iOpt,1:5))
)
,prevRes=prevRes
)
return(ans)
}
partGLFitLRC <- function(
### Optimize rectangular hyperbolic light response curve against data in one window and estimate uncertainty
dsDay ##<< data.frame with columns NEE, Rg, Temp_C, VPD, and no NAs in NEE
, E0 ##<< temperature sensitivity of respiration
, sdE0 ##<< standard deviation of E_0.n
, RRefNight ##<< basal respiration estimated from night time data, used as starring value and prior estimate
, controlGLPart=partGLControl() ##<< further default parameters (see \code{\link{partGLControl}})
, lastGoodParameters ##<< numeric vector of last good theta
){
##author<< TW, MM
##seealso<< \code{\link{partGLFitLRCWindows}}
#
#Definition of initial guess theta, theta2 and theta3. Three initial guess vectors are defined according to Lasslop et al., 2010
namesPars <- c("k","beta0", "alfa", "Rb","E0" )
nPar <- length(namesPars)
##details<<
## Optimization is performed for three initial parameter sets that differ by beta0 (*1.3, *0.8).
## From those three, the optimization result is selected that yielded the lowest misfit.
## Starting values are: k=0, beta=interpercentileRange(0.03,0.97) of respiration, alpha=0.1, R_ref
## from nightTime estimate.
## E0 is fixed to the night-time estimate, but varied for estimating parameter uncertainty.
theta.V.n<-matrix(NA, 3,nPar, dimnames=list(NULL,namesPars))
parameterPrior <- theta.V.n[1,] <- theta.V.n[2,] <- theta.V.n[3,] <- c(
k=0
#,beta=as.vector(abs(quantile(dsDay$NEE, 0.03)-quantile(dsDay$NEE, 0.97)))
,beta=as.vector(abs(quantile(dsDay$NEE, 0.03, na.rm=TRUE)-quantile(dsDay$NEE, 0.97, na.rm=TRUE)))
,alpha=0.1
#,R_ref=mean(NEENight.V.n, na.rm=T)
,R_ref=if( is.finite(RRefNight) ) as.vector(RRefNight) else stop("must provide finite RRefNight") #mean(NEENight.V.n, na.rm=T)
,E_0=as.vector(E0)
) #theta [numeric] -> parameter vector (theta[1]=kVPD, theta[2]-beta0, theta[3]=alfa, theta[4]=Rref)
#twutz: beta is quite well defined, so try not changing it too much
theta.V.n[2,2] <- parameterPrior[2]*1.3
theta.V.n[3,2] <- parameterPrior[2]*0.8
#
##seealso<< \code{\link{parmGLOptimLRCBounds}}
resOpt3 <- apply( theta.V.n, 1, function(theta0){
resOpt <- parmGLOptimLRCBounds(theta0, parameterPrior, dsDay=dsDay, ctrl=controlGLPart, lastGoodParameters.V.n=lastGoodParameters )
})
iValid <- which(sapply(resOpt3,function(resOpt){ is.finite(resOpt$theta[1]) }))
resOpt3Valid <- resOpt3[iValid]
optSSE <- sapply(resOpt3Valid, "[[", "value")
if( sum(!is.na(optSSE)) == 0L ){
# none of the intial fits yielded valid result, create NA-return values
opt.parms.V <- structure( rep(NA_real_, nPar), names=colnames(theta.V.n) )
covParms <- matrix(NA_real_, nPar, nPar, dimnames=list(colnames(theta.V.n),colnames(theta.V.n)))
resOpt <- list(iOpt = integer(0))
} else {
resOpt <- resOpt3Valid[[iBest <- which.min(optSSE)]] # select the one with the least cost
opt.parms.V<-resOpt$theta
if(controlGLPart$nBootUncertainty == 0L) {
##details<< If \code{controlGLPart.l$nBootUncertainty == 0L} then the covariance matrix of the
## parameters is estimated by the Hessian of the LRC curve at optimum.
## Then, the additional uncertainty and covariance with uncertaint E0 is neglected.
#seParmsHess <- seParmsHess0 <- sqrt(abs(diag(solve(resOpt$hessian))))
covParmsLRC <- if( (resOpt$hessian[1L,1L] < 1e-8) ){
# case where k=0 and not varying: cov(k,:) = cov(:,) = 0
covParmsLRC <- structure( diag(0,nrow=nrow(resOpt$hessian)), dimnames=dimnames(resOpt$hessian))
covParmsLRC[-1L,-1L] <- solve(resOpt$hessian[-1L,-1L])
covParmsLRC
} else {
solve(resOpt$hessian)
}
covParms <- structure( diag(0,nrow=length(resOpt$theta)), dimnames=list(namesPars,namesPars))
covParms[5L,5L] <- sdE0^2
covParms[resOpt$iOpt,resOpt$iOpt] <- covParmsLRC
if( any(diag(covParms) < 0)) opt.parms.V[] <- NA # no real if covariance negative
} else {
##details<<
## If \code{controlGLPart.l$nBootUncertainty > 0L} then the covariance matrix of the
## parameters is estimated by a bootstrap of the data.
## In each draw, E0 is drawn from N ~ (E_0, sdE_0).
# #seealso<< \code{\link{.bootStrapLRCFit}}
resBoot <- .bootStrapLRCFit(resOpt$theta, resOpt$iOpt, dsDay, sdE0, parameterPrior, controlGLPart)
#resBoot <- .bootStrapLRCFit(resOpt$theta, resOpt$iOpt, dsDay, sdE_0.n, parameterPrior, controlGLPart.l=within(controlGLPart.l,nBootUncertainty <- 30L))
covParms <- cov(resBoot)
#better not to average parameters
#opt.parms.V <- apply(resBoot, 2, median, na.rm=TRUE)
#se.parms.V <- apply(resBoot, 2, sd, na.rm=TRUE)
}
}
##value<< a list, If none of the optimizations from different starting conditions converged,
## the parameters are NA
ans <- list(
opt.parms.V=opt.parms.V ##<< numeric vector of optimized parameters including
## the fixed ones and including E0
,iOpt=c(resOpt$iOpt,5L) ##<< index of parameters that have been optimized
## , here including E0, which has been optimized prior to this function.
,initialGuess.parms.V.n=theta.V.n[1,] ##<< the initial guess from data
,covParms=covParms ##<< numeric matrix of the covariance matrix of parameters, including E0
)
}
.tmp.f <- function(){
plot(dsDay$Rg,-1*dsDay$NEE)
thetaB <- resOpt$theta
#thetaB <- opt.parms.V
#thetaB <- opt[[iBest <- which.min(optSSE)]]
#
#thetaB <- resOpt$par
#thetaB <- c(resOpt3[[1]]$par, E_0=as.vector(E_0.n))
#thetaB <- theta
#thetaB <- thetaOrig
points(
#lines(
dsDay$Rg, tmp <- partGL_RHLightResponse(thetaB
,Rg = dsDay$Rg
,Temp=dsDay$Temp
,VPD = dsDay$VPD
#,E0 = E_0.n
,fixVPD = FALSE
)$NEP, col="red")
points(dsDay$Rg, tmp <- partGL_RHLightResponse(thetaB
,Rg = dsDay$Rg
,Temp=dsDay$Temp
,VPD = dsDay$VPD
,E0 = E_0.n
,fixVPD = FALSE
)$Reco, col="green")
#)$GPP, col="blue")
}
parmGLOptimLRCBounds <- function(
### Optimize parameters of light response curve and refit with some fixed parameters if fitted parameters are outside bounds
theta0 ##<< initial parameter estimate
,parameterPrior ##<< prior estimate of model parameters
, ... ##<< further parameters to \code{.optimLRC}, such as \code{dsDay}, and \code{ctrl}
,lastGoodParameters.V.n ##<< parameters vector of last successful fit
){
##author<< TW, MM
##seealso<< \code{\link{partGLFitLRC}}
if( !is.finite(lastGoodParameters.V.n[3L]) ) lastGoodParameters.V.n[3L] <- 0.22 # twutz 161014: default alpha
isAlphaFix <- FALSE
isKFix <- FALSE
resOpt <- resOpt0 <- .optimLRC(theta0, isUsingFixedVPD=isKFix, isUsingFixedAlpha=isAlphaFix, parameterPrior = parameterPrior, ... )
# positions in theta0: "k" "beta0" "alfa" "Rb" "E0"
# IF kVPD parameter less or equal zero then estimate the parameters withouth VPD effect
##details<<
## If parameters alpha or k are outside bounds (Table A1 in Lasslop 2010), refit with some parameters fixed
## to values from fit of previous window.
theta0Adj <- theta0 # intial parameter estimate with some parameters adjusted to bounds
#dsDay <- list(...)$dsDay
if (resOpt$theta[1L] < 0){
isKFix <- TRUE
theta0Adj[1L] <- 0
resOpt <- .optimLRC(theta0Adj, isUsingFixedVPD=isKFix, isUsingFixedAlpha=isAlphaFix, parameterPrior = parameterPrior, ... )
# check alpha, if less than zero estimate parameters with fixed alpha of last window
if ( (resOpt$theta[3L] > 0.22) && is.finite(lastGoodParameters.V.n[3L]) ){
isAlphaFix <- TRUE
theta0Adj[3L] <- lastGoodParameters.V.n[3L]
resOpt <- .optimLRC(theta0Adj, isUsingFixedVPD=isKFix, isUsingFixedAlpha=isAlphaFix, parameterPrior = parameterPrior, ... )
}
} else {
# check alpha, if gt 0.22 estimate parameters with fixed alpha of last window
# if not last window exists, let alpha > 0.22
if ( (resOpt$theta[3L] > 0.22) && is.finite(lastGoodParameters.V.n[3L]) ){
isAlphaFix <- TRUE
theta0Adj[3L] <- lastGoodParameters.V.n[3L]
resOpt <- .optimLRC(theta0Adj, isUsingFixedVPD=isKFix, isUsingFixedAlpha=isAlphaFix, parameterPrior = parameterPrior, ... )
# check k, if less than zero estimate parameters without VPD effect and with fixed alpha of last window
if (resOpt$theta[1L] < 0){
isKFix <- TRUE
theta0Adj[1L] <- 0
resOpt <- .optimLRC(theta0Adj, isUsingFixedVPD=isKFix, isUsingFixedAlpha=isAlphaFix, parameterPrior = parameterPrior, ... )
}
}
}
##details<<
## No parameters are reported if alpha<0 or Rb < 0 or beta0 < 0 or beta0 > 250
# positions in theta0: "k" "beta0" "alfa" "Rb" "E0"
if( !is.na(resOpt$theta[1L]) && ((resOpt$theta[3L] < 0) || (resOpt$theta[4L] < 0) || (resOpt$theta[2L] < 0) || (resOpt$theta[2L] >= 250)) ){
# TODO estimate Rb from daytime data?
#LloydT_E0fix
#stop("case with alpha or beta < 0")
resOpt$theta[] <- NA
}
##details<<
## No parameters are reported if beta0 > 4*initialEstimate, to avoid cases where data is far away from saturation.
if( isTRUE(as.vector(resOpt$theta[2L] > 4*parameterPrior[2L])) ){
resOpt$theta[] <- NA
}
##value<< list result of optimization as of \code{.optimLRC} with entries
## \item{theta}{ numeric parameter vector that includes the fixed components}
## \item{iOpt}{ integer vector of indices of the vector that have been optimized}
resOpt
}
.bootStrapLRCFit <- function(
### Compute parameters uncertainty by bootstrap
theta0, iOpt, dsDay, sdE_0.n, parameterPrior, controlGLPart.l
){
##value<<
## matrix with each row a parameter estimate on a different bootstrap sample
ans <-matrix(NA, nrow=controlGLPart.l$nBootUncertainty, ncol=length(theta0), dimnames=list(NULL,names(theta0)))
isUsingFixedVPD <- !(1L %in% iOpt)
isUsingFixedAlpha <- !(3L %in% iOpt)
##details<<
## In addition to resampling the original data, also the temperature sensitivity is resampled
## from its uncertainty distribution.
E0r <- rnorm( controlGLPart.l$nBootUncertainty, theta0[5L], sdE_0.n )
E0 <- pmax(50,pmin(400,E0r))
theta <- theta0
#iBoot <- 1L
for (iBoot in c(1:controlGLPart.l$nBootUncertainty)){
idx <- sample(nrow(dsDay), replace=TRUE)
dsDayB <- dsDay[idx,]
theta[5L] <- E0[iBoot]
resOptBoot <- .optimLRC(theta, isUsingFixedVPD=isUsingFixedVPD, isUsingFixedAlpha=isUsingFixedAlpha
, dsDayB, parameterPrior, controlGLPart.l )
if( resOptBoot$convergence == 0L )
#TODO: also remove the bery bad cases?
ans[iBoot,]<-resOptBoot$theta
}
ans
}
.optimLRC <- function(
###<< one fit of the light response curve with a subset of parameters depending on which are fixed
theta ##<< numeric vector of starting values
, isUsingFixedVPD=FALSE, isUsingFixedAlpha=FALSE ##<< selection of which parameters to optimize
, dsDay ##<< dataframe of NEE, sdNEE and predictors Rg, VPD and Temp
#, E_0.n ##<< temperature sensitivity
, parameterPrior ##<< prior parameter estimates
, ctrl ##<< list of further controls
){
if( !all(is.finite(theta))) stop("need to provide finite starting values.")
##details<<
## Only those records are used for optimization where both NEE and sdNEE are finite.
dsDayFinite <- dsDay[ is.finite(dsDay$NEE) & is.finite(dsDay$sdNEE), ]
##details<<
## Optimization of Light-response curve parameters takes into account the uncertainty of the flux values.
## In order to avoid very strong leverage, values with a very low uncertainty (< median) are assigned
## the median of the uncertainty.
## This procedure downweighs records with a high uncertainty, but does not apply a large leverage for
## records with a very low uncertainty. Avoid this correction by suppyling setting \code{ctrl$isBoundLowerNEEUncertainty = FALSE}
Fc_unc <- if( isTRUE(ctrl$isBoundLowerNEEUncertainty) ){
pmax( dsDayFinite$sdNEE, quantile(dsDayFinite$sdNEE, 0.3) ) #twutz: avoid excessive weights by small uncertainties (of 1/unc^2)
} else {
dsDayFinite$sdNEE
}
#plot( Fc_unc ~ dsDayFinite$sdNEE)
##details<<
## The beta parameter is quite well defined. Hence use a prior with a standard deviation.
## The specific results are sometimes a bit sensitive to the uncertainty of the beta prior.
## This uncertainty is set corresponding to 10 times the median relative flux uncertainty.
## The prior is weighted n times the observations in the cost.
## Hence, overall it is using a weight of 1/10 of the weight of all observations.
sdParameterPrior <- if(ctrl$isLasslopPriorsApplied){
c(k=50, beta=600, alpha=10, Rb=80, E0=NA)
} else {
medianRelFluxUncertainty <- abs(median(Fc_unc/dsDayFinite$NEE))
sdBetaPrior <- 10*medianRelFluxUncertainty*parameterPrior[2]/sqrt(nrow(dsDayFinite))
c(k=NA, beta=as.vector(sdBetaPrior), alpha=NA, Rb=NA, E0=NA)
}
isUsingHessian <- (ctrl$nBootUncertainty==0L)
iOpt <-
if( !isUsingFixedVPD & !isUsingFixedAlpha ) 1:4 else
if( isUsingFixedVPD & !isUsingFixedAlpha ) 2:4 else
if( !isUsingFixedVPD & isUsingFixedAlpha ) c(1L,2L,4L) else
if( isUsingFixedVPD & isUsingFixedAlpha ) c(2L,4L)
sdParameterPrior[-iOpt] <- NA
# encountered the case where fit ran into a local optimum at very low (even negative) beta
# hence first do a fit with log(beta) or a fit with very strong prior on beta and alpha
# and relax priors afterwards.
# thetaLog <- theta; thetaLog[2] <- log(theta[2])
# parameterPriorLog <- parameterPrior; parameterPriorLog[2] <- log(parameterPrior[2])
# # for the log-Fit assume a strong prior on beta to avoid local minimum at beta=0
# sdParameterPriorLog <- sdParameterPrior #; sdParameterPriorLog[2] <- log(parameterPrior[2]) - log(parameterPrior[2] - sdParameterPrior[2])
# resOptimLog <- optim(thetaLog[iOpt], .partGLRHLightResponseCostLogBeta
# #tmp <- .partGLRHLightResponseCost( theta[iOpt],
# ,theta=thetaLog
# ,iOpt=iOpt
# ,flux = -dsDayFinite$NEE
# ,sdFlux = Fc_unc
# ,parameterPrior = parameterPriorLog
# ,sdParameterPrior = sdParameterPriorLog
# ,Rg = dsDayFinite$Rg
# ,VPD = dsDayFinite$VPD
# ,Temp=dsDayFinite$Temp
# #,E0 = E_0.n
# ,control=list(reltol=ctrl$LRCFitConvergenceTolerance)
# ,method="BFGS", hessian=isUsingHessian)
# thetaOrig <- theta; thetaOrig[iOpt] <- resOptimLog$par; thetaOrig[2] <- exp(thetaOrig[2])
#
thetaOrig <- theta
sdStrongPrior <- sdParameterPrior; sdStrongPrior[2] <- sdParameterPrior[2]/10; sdStrongPrior[3] <- 0.5
resOptimStrongPrior <- optim(thetaOrig[iOpt], .partGLRHLightResponseCost
#tmp <- .partGLRHLightResponseCost( theta[iOpt],
,theta=thetaOrig
,iOpt=iOpt
,flux = -dsDayFinite$NEE
,sdFlux = Fc_unc
,parameterPrior = parameterPrior
,sdParameterPrior = sdStrongPrior
,Rg = dsDayFinite$Rg
,VPD = dsDayFinite$VPD
,Temp=dsDayFinite$Temp
#,E0 = E_0.n
,control=list(reltol=ctrl$LRCFitConvergenceTolerance)
,method="BFGS", hessian=isUsingHessian)
thetaOrig[iOpt] <- resOptimStrongPrior$par
#
resOptim <- optim(thetaOrig[iOpt], .partGLRHLightResponseCost
#resOptim <- optim(theta, .partGLRHLightResponseCost
#tmp <- .partGLRHLightResponseCost( theta[iOpt]
,theta=thetaOrig
,iOpt=iOpt
,flux = -dsDayFinite$NEE
,sdFlux = Fc_unc
,parameterPrior = parameterPrior
,sdParameterPrior = sdParameterPrior
,Rg = dsDayFinite$Rg
,VPD = dsDayFinite$VPD
,Temp=dsDayFinite$Temp
#,E0 = E_0.n
,control=list(reltol=ctrl$LRCFitConvergenceTolerance)
,method="BFGS", hessian=isUsingHessian)
##value<<
## list of restult of \code{\link{optim}} amended with list
thetaOpt <- theta; thetaOpt[iOpt] <- resOptim$par
ans <- list(
theta = thetaOpt ##<< numeric vector: optimized parameter vector including the fixed components
,iOpt = iOpt ##<< integer vector: position of parameters that have been optimized
)
c(resOptim, ans)
}
.tmp.f <- function(){
tmp <- .partGLRHLightResponseCost(theta
,flux = -dsDay$NEE
,sdFlux = Fc_unc
,betaPrior = betaPrior
,sdBetaPrior = sdBetaPrior
,Rg = dsDay$Rg
,VPD = dsDay$VPD
,Temp=dsDay$Temp
,E_0.n = E_0.n
)
}
partGL_RHLightResponse <- function(
###Rectungular Hyperbolic Light Response function: (Xu & Baldocchi, 2004; Falge et al., 2001; Lasslop et al., 2010)
theta ##<< theta [numeric] -> parameter vector (theta[1]=kVPD (k), theta[2]=beta0 (beta), theta[3]=alfa, theta[4]=Rref (rb), theta[5]=E0)
##<< E0: Temperature sensitivity ("activation energy") in Kelvin (degK) #get("testparams", envir=environment(foo)
,Rg ##<< ppfd [numeric] -> photosynthetic flux density [umol/m2/s] or Global Radiation
,VPD ##<< VPD [numeric] -> Vapor Pressure Deficit [hPa]
,Temp ##<< Temp [degC] -> Temperature [degC]
#,E0 ##<< Temperature sensitivity ("activation energy") in Kelvin (degK) #get("testparams", envir=environment(foo)
,VPD0 = 10 ##<< VPD0 [hPa] -> Parameters VPD0 fixed to 10 hPa according to Lasslop et al 2010
,fixVPD = FALSE ##<< fixVPD TRUE or FALSE -> if TRUE the VPD effect is not considered
) {
##details<<
## with fixes E_0.n for the estimation of Rd
## VPD effect included according to Lasslop et al., 2010
##details<<
## If theta is a matrix, a different row of parameters is used for different entries of other inputs
if( is.matrix(theta) ){
kVPD<-theta[,1]
beta0<-theta[,2]
alfa<-theta[,3]
Rref<-theta[,4]
E0<-theta[,5]
} else {
kVPD<-theta[1]
beta0<-theta[2]
alfa<-theta[3]
Rref<-theta[4]
E0<-theta[5]
}
Amax <- if( isTRUE(fixVPD) ) beta0 else {
ifelse(VPD > VPD0, beta0*exp(-kVPD*(VPD-VPD0)), beta0)
}
#Reco<-Rref*exp(E0*(1/((273.15+10)-227.13)-1/(Temp+273.15-227.13))) # formerly at 10degC
Reco<-Rref*exp(E0*(1/((273.15+15)-227.13)-1/(Temp+273.15-227.13)))
GPP <- (Amax*alfa*Rg)/(alfa*Rg+Amax)
NEP <- GPP - Reco
## a data.frame of length of Rg of computed
ans <- list(
NEP=NEP ##<< Net ecosystem production (-NEE), vector of length(Rg)
,Reco=Reco ##<< Ecosystem respiration
,GPP=GPP ##<< Gross primary production
)
}
partGL_RHLightResponseGrad <- function(
### Gradient of \code{\link{partGL_RHLightResponse}}
theta ##<< theta [numeric] -> parameter vector (theta[1]=kVPD (k), theta[2]=beta0 (beta), theta[3]=alfa, theta[4]=Rref (rb), theta[4]=E0)
##<< E0: Temperature sensitivity ("activation energy") in Kelvin (degK)
,Rg ##<< ppfd [numeric] -> photosynthetic flux density [umol/m2/s] or Global Radiation
,VPD ##<< VPD [numeric] -> Vapor Pressure Deficit [hPa]
,Temp ##<< Temp [degC] -> Temperature [degC]
#,E0 ##<< Temperature sensitivity ("activation energy") in Kelvin (degK) #get("testparams", envir=environment(foo)
,VPD0 = 10 ##<< VPD0 [hPa] -> Parameters VPD0 fixed to 10 hPa according to Lasslop et al 2010
,fixVPD = FALSE ##<< fixVPD TRUE or FALSE -> if TRUE the VPD effect is not considered
) {
if( is.matrix(theta) ){
kVPD<-theta[,1]
beta0<-theta[,2]
alfa<-theta[,3]
Rref<-theta[,4]
E0<-theta[,5]
} else {
kVPD<-theta[1]
beta0<-theta[2]
alfa<-theta[3]
Rref<-theta[4]
E0<-theta[5]
}
Amax <- if( isTRUE(fixVPD) ) beta0 else {
ifelse(VPD > VPD0, beta0*exp(-kVPD*(VPD-VPD0)), beta0)
}
#ex <- expression( beta0*exp(-kVPD*(VPD-VPD0)) ); deriv(ex,c("beta0","kVPD"))
dAmax_dkVPD <- if( isTRUE(fixVPD) ) 0 else {
ifelse(VPD > VPD0, beta0*-(VPD-VPD0)*exp(-kVPD*(VPD-VPD0)), 0)
}
dAmax_dbeta0 <- if( isTRUE(fixVPD) ) 0 else {
ifelse(VPD > VPD0, exp(-kVPD*(VPD-VPD0)), 1)
}
#Reco<-Rref*exp(E0*(1/((273.15+10)-227.13)-1/(Temp+273.15-227.13)))
#ex <- expression( Rref*exp(E0*(1/((273.15+10)-227.13)-1/(Temp+273.15-227.13))) ); deriv(ex,c("Rref","E0"))
#.expr7 <- 1/(273.15 + 10 - 227.13) - 1/(Temp + 273.15 - 227.13)
.expr7 <- 1/(273.15 + 15 - 227.13) - 1/(Temp + 273.15 - 227.13)
.expr9 <- exp(E0 * .expr7)
gradReco <- matrix(0, ncol=2L, nrow=length(.expr9), dimnames=list(NULL,c("Rref","E0")))
gradReco[,"Rref"] <- dReco_dRRef <- .expr9
gradReco[,"E0"] <- dReco_dE0 <- Rref * (.expr9 * .expr7)
#GPP <- (Amax*alfa*Rg)/(alfa*Rg+Amax)
#ex <- expression( (Amax*alfa*Rg)/(alfa*Rg+Amax) ); deriv(ex,c("Amax","alfa"))
.expr2 <- Amax * alfa * Rg
.expr3 <- alfa * Rg
.expr4 <- .expr3 + Amax
.expr7 <- .expr4^2
.value <- .expr2/.expr4
gradGPP <- array(0, c(length(.value), 3L), list(NULL, c("kVPD","beta0","alfa")))
dGPP_dAMax <- .expr3/.expr4 - .expr2/.expr7
gradGPP[, "beta0"] <- dGPP_dAMax * dAmax_dbeta0
gradGPP[, "kVPD"] <- dGPP_dAMax * dAmax_dkVPD
gradGPP[, "alfa"] <- Amax * Rg/.expr4 - .expr2 * Rg/.expr7
#NEP <- GPP - Reco
gradNEP <- cbind(gradGPP, -gradReco)
## list with gradient matrices. For each record (length(Rg)), c("kVPD","beta0","alfa","Rref")
ans <- list(
NEP=gradNEP
,Reco=gradReco
,GPP=gradGPP
)
}
.partGLRHLightResponseCost <- function(
### Computing residual sum of sqares for predictions vs. data of NEE
thetaOpt ##<< parameter vecotr with components of theta0 that are optimized
,theta ##<< parameter vector with positions as in argument of \code{\link{partGL_RHLightResponse}}
,iOpt ##<< position in theta0 that are optimized
,flux=NA ##<< numeric: NEP (-NEE) or GPP time series [umolCO2/m2/s], should not contain NA
,sdFlux=NA ##<< numeric: standard deviation of Flux [umolCO2/m2/s], should not contain NA
,parameterPrior ##<< numeric vector along theta: prior estimate of parameter (range of values)
,sdParameterPrior ##<< standard deviation of parameterPrior
,... ##<< other arguments to \code{\link{partGL_RHLightResponse}}
,VPD0 = 10 ##<< VPD0 [hPa] -> Parameters VPD0 fixed to 10 hPa according to Lasslop et al 2010
,fixVPD = FALSE ##<< fixVPD TRUE or FALSE -> if TRUE the VPD effect is not considered
,useCVersion=TRUE ##<< set to FALSE to use R instead of fast C version, e.g. for debugging
# also tried to pass dsDay data.frame instead of all the variables separately, but accessing data.frame
# columns in the cost function was a severe penalty (here double execution time)
) {
theta[iOpt] <- thetaOpt
#if( FALSE ){
if( useCVersion){
# generated in RcppExports: RHLightResponseCostC <- function(theta, flux, sdFlux, parameterPrior, sdParameterPrior, Rg, VPD, Temp, VPD0, fixVPD) {
# when generating add description
### Compute the cost of RLightResponse prediction versus observations by using fast C-code
RHLightResponseCostC(theta, flux, sdFlux, parameterPrior, sdParameterPrior, ..., VPD0=VPD0, fixVPD=fixVPD)
} else {
resPred <- partGL_RHLightResponse(theta, ..., VPD0=VPD0, fixVPD=fixVPD)
NEP_mod <- resPred$NEP
misFitPrior <- (((theta - parameterPrior))/(sdParameterPrior))^2
misFitObs <- sum(((NEP_mod-flux)/sdFlux)^2)
RSS <- misFitObs + sum(misFitPrior, na.rm=TRUE)
#if( !is.finite(RSS) ) recover() # debugging the fit
RSS
}
}
.partGLRHLightResponseCostLogBeta <- function(
### Computing residual sum of sqares for predictions vs. data of NEE, with parameter beta given at log-Scale
thetaOpt ##<< parameter vecotr with components of theta0 that are optimized
,theta ##<< parameter vector with positions as in argument of \code{\link{partGL_RHLightResponse}}
,iOpt ##<< position in theta0 that are optimized
,flux=NA ##<< numeric: NEP (-NEE) or GPP time series [umolCO2/m2/s], should not contain NA
,sdFlux=NA ##<< numeric: standard deviation of Flux [umolCO2/m2/s], should not contain NA
,parameterPrior ##<< numeric vector along theta: prior estimate of parameter (range of values)
,sdParameterPrior ##<< standard deviation of parameterPrior
,... ##<< other arguments to \code{\link{partGL_RHLightResponse}}
,VPD0 = 10 ##<< VPD0 [hPa] -> Parameters VPD0 fixed to 10 hPa according to Lasslop et al 2010
,fixVPD = FALSE ##<< fixVPD TRUE or FALSE -> if TRUE the VPD effect is not considered
,useCVersion=TRUE ##<< set to FALSE to use R instead of fast C version, e.g. for debugging
# also tried to pass dsDay data.frame instead of all the variables separately, but accessing data.frame
# columns in the cost function was a severe penalty (here double execution time)
) {
thetaOrig <- theta
thetaOrig[iOpt] <- thetaOpt
thetaOrig[2] <- exp(thetaOrig[2])
#if( FALSE ){
if( useCVersion){
RHLightResponseCostC(thetaOrig, flux, sdFlux, parameterPrior, sdParameterPrior, ..., VPD0=VPD0, fixVPD=fixVPD)
} else {
resPred <- partGL_RHLightResponse(thetaOrig, ..., VPD0=VPD0, fixVPD=fixVPD)
NEP_mod <- resPred$NEP
misFitPrior <- (((thetaOrig - parameterPrior))/(sdParameterPrior))^2
misFitObs <- sum(((NEP_mod-flux)/sdFlux)^2)
RSS <- misFitObs + sum(misFitPrior, na.rm=TRUE)
#if( !is.finite(RSS) ) recover() # debugging the fit
RSS
}
}
partGLBoundParameters <- function(
### Check if parameters are in the range and correct to bounds
resOpt ##<< list with first entry vector of parameters, second entry its uncertainty
, last_good ##<< vector of windows most recent good fit, make sure to set it to reasonable values, e.g. initial guess before optimiztion
){
##details<<
## adjustments according to TAble A1 in Lasslop 2010
## vector positions: 1: k, 2: beta, 3: alpha, 4: rb
#
# Table A1 Initial guess; Valid Range; If outside range
# E_0.n 100; 50-400; Set to value of previous window, if no previous window exists estimates <50 were set to 50, estimates >400 were set to 400
# rb Mean of nighttime NEE >0; Whole parameter set is not used
# alpha 0.01; ?0,<0.22; Set to value of previous window, if no previous window exists and <0, set to zero
# beta0 Abs (0.03quantile - 0.97quantile) of NEE ?0; <250; If >100 then ? (?)<? If negative set to zero, else the whole parameter set is not used
# kVPD ;0 ?0; Set to zero
#parms_out_range<-0 #IF set to 1 means that the parameters are outside range and no computation uncertainties
opt.parms.V <- resOpt[[1]]
se.parms.V <- resOpt[[2]]
isGoodParameterSet <- TRUE # FALSE means parameters are outside range and no uncertainties are reported
if( any(is.na(opt.parms.V)) ){
isGoodParameterSet <- FALSE
}
# isTRUE(as.vector()) returns FALSE also for NA
if ( isTRUE((opt.parms.V[1] < 0)) ){ # k
opt.parms.V[1]<-0
isGoodParameterSet <- FALSE
}
if ( isTRUE(as.vector(opt.parms.V[2] < 0)) ){ # beta
opt.parms.V[2]<-0
isGoodParameterSet <- FALSE
}
if ( isTRUE(as.vector(opt.parms.V[3] <= 0 || opt.parms.V[3] > 0.22))){ #set to alpha values of the latest good parameters set
opt.parms.V[3]<-last_good[3]
isGoodParameterSet <- FALSE
}
#whole par set not used, if beta > 250 or (beta > 100 and sdBeta >= beta) or rb < 0
if ( isTRUE(as.vector(opt.parms.V[2] > 250 || (opt.parms.V[2] > 100 && se.parms.V[2] >= opt.parms.V[2] ) ))){
opt.parms.V[]<-NA
isGoodParameterSet <- FALSE
}
if ( isTRUE(as.vector(opt.parms.V[4] < 0 ))){
opt.parms.V[]<-NA
isGoodParameterSet <- FALSE
}
if( !isGoodParameterSet){
se.parms.V[]<-NA
}
##value<< list with entries
list(
opt.parms.V = opt.parms.V ##<< numeric vector of bounded optimized parameters
,se.parms.V = se.parms.V ##<< numeric vector of uncertainties of parameters
,isGoodParameterSet = isGoodParameterSet ##<< boolean scalar that is FALSE if original parameters were outside bounds
)
}
partGLInterpolateFluxes <- function(
### Interpolate ecoystem respiration (Reco) and Gross primary production (GPP) and associated uncertainty from two neighboring parameter sets of Light response curves
Rg ##<< photosynthetic flux density [umol/m2/s] or Global Radiation
,VPD ##<< Vapor Pressure Deficit [hPa]
,Temp ##<< Temperature [degC]
,resParms ##<< data frame with results of \code{\link{partGLFitLRCWindows}} of fitting the light-response-curve for several windows
,controlGLPart=partGLControl() ##<< further default parameters, see \code{\link{partGLControl}}
){
##author<< TW
##seealso<< \code{link{partitionNEEGL}}
##details<<
## \code{resLRC$iFirstRecInCentralDay} must denote the row for which the LRC parameters are representative,
## here, the first record of the center day
# create a dataframe with index of rows of estimates before and after and correponding weights
summaryLRC <- resParms$summary
nLRC <- nrow(summaryLRC)
nRec <- length(Rg)
Temp_Kelvin <- Temp+273.15
if( isTRUE(controlGLPart$isAssociateParmsToMeanOfValids) ){
# there might be several rows with the same iMeanRec, omit those rows unless the first of each reoccuring iMeanRec
tabMeanRec <- table(summaryLRC$iMeanRec)
iRecsDouble <- as.integer(names(tabMeanRec[ tabMeanRec > 1L ]))
iRecsOmit <- do.call(c, lapply( iRecsDouble, function(iRecDouble){
which(summaryLRC$iMeanRec==iRecDouble)[-1L]
}))
if( length(iRecsOmit)) summaryLRC <- summaryLRC[-iRecsOmit,]
}
##details<<
## Parameter estimates are reported for the first record of the window, or
## if \code{isTRUE(controlGLPart.l$isAssociateParmsToMeanOfValids)} for the mean time of all valid records within the window
# for each original record merge parameters assicated with previous fit or next fit respectively
colNameAssoc <- if( isTRUE(controlGLPart$isAssociateParmsToMeanOfValids) ) "iMeanRec" else "iCentralRec"
dsAssoc <- .partGPAssociateSpecialRows(summaryLRC[[colNameAssoc]],nRec)
# now we have iBefore and iAfter
dsBefore <- merge(
structure(data.frame(dsAssoc$iSpecialBefore, dsAssoc$iBefore),names=c("iParRec",colNameAssoc))
, summaryLRC[,c(colNameAssoc,"R_ref","E_0","a","b","k")]
)
dsAfter <- merge( structure(data.frame(dsAssoc$iSpecialAfter, dsAssoc$iAfter),names=c("iParRec",colNameAssoc)), summaryLRC[,c(colNameAssoc,"R_ref","E_0","a","b","k")])
if( (nrow(dsBefore) != nRec) || (nrow(dsAfter) != nRec)) stop("error in merging parameters to original records.")
Reco2 <- lapply( list(dsBefore,dsAfter), function(dsi){
tmp <- fLloydTaylor(dsi$R_ref, dsi$E_0, Temp_Kelvin, T_ref.n=273.15+15)
})
#dsi <- dsBefore
GPP2 <- lapply( list(dsBefore,dsAfter), function(dsi){
theta <- as.matrix(dsi[,c("k","b","a","R_ref","E_0")])
tmp <- partGL_RHLightResponse(theta, Rg, VPD, Temp=Temp)$GPP
})
# interpolate between previous and next fit, weights already sum to 1
Reco <- (dsAssoc$wBefore*Reco2[[1]] + dsAssoc$wAfter*Reco2[[2]]) #/ (dsAssoc$wBefore+dsAssoc$wAfter)
GPP <- (dsAssoc$wBefore*GPP2[[1]] + dsAssoc$wAfter*GPP2[[2]]) #/ (dsAssoc$wBefore+dsAssoc$wAfter)
#recover() # to inspect the deviations between successive estimates
ans <- ansPred <- data.frame(
Reco = Reco
,GPP = GPP
)
if( isTRUE(controlGLPart$isSdPredComputed)){
varPred2 <- lapply( list(dsBefore,dsAfter), function(dsi){
theta <- as.matrix(dsi[,c("k","b","a","R_ref","E_0")])
grad <- partGL_RHLightResponseGrad(theta, Rg, VPD, Temp)
varPred <- matrix(NA_real_, nrow=nrow(dsi), ncol=2L, dimnames=list(NULL,c("varGPP","varReco")))
#iRec <- 1L
for( iRec in 1:nrow(dsi)){
iParRec <- dsi$iParRec[iRec]
# get the fitting object, TODO better document
resOpt <- resParms$resOptList[[ iParRec ]]
gradGPP <- grad$GPP[iRec,]
gradReco <- grad$Reco[iRec,]
# make sure parameter names match postions in covParms
varPred[iRec,1L] <- varGPP <- gradGPP %*% resOpt$covParms[1:3,1:3] %*% gradGPP
varPred[iRec,2L] <- varReco <- gradReco %*% resOpt$covParms[4:5,4:5] %*% gradReco
}
varPred
})
sdGPP <- sqrt(dsAssoc$wBefore^2*varPred2[[1]][,"varGPP"] + dsAssoc$wAfter^2*varPred2[[1]][,"varGPP"])
sdReco <- sqrt(dsAssoc$wBefore^2*varPred2[[1]][,"varReco"] + dsAssoc$wAfter^2*varPred2[[1]][,"varReco"])
ans <- cbind(ansPred, data.frame(sdGPP=sdGPP, sdReco=sdReco))
}
# compute differences in rows to next parameter estimate
dRecBefore <- dsBefore[[colNameAssoc]] - 1:nrow(dsBefore)
dRecAfter <- dsAfter[[colNameAssoc]] - 1:nrow(dsBefore)
isBeforeCloser <- abs(dRecBefore) <= abs(dRecAfter)
ans$dRecNextEstimate <- ifelse(isBeforeCloser, dRecBefore, dRecAfter)
##value<< data.frame with nrow() rows and columns GPP, Reco, varGPP, varReco, and dRecNextEstimate
ans
}
.tmp.f <- function(){
# execute on recover jsut before ans <- in .partGPInterploateFluxes, from test_that("interpolate Fluxes"
# there is only one period with changing parametrization
x <- seq_along(Rg)
plot( Reco ~ x )
lines( Reco2[[1]] ~ x, col="blue" )
lines( Reco2[[2]] ~ x, col="orange" )
# note how the black dots first follow the blue line, then converge to the orange line
plot( GPP ~ x )
lines( GPP2[[1]] ~ x, col="blue" )
lines( GPP2[[2]] ~ x, col="orange" )
# not much difference between blue and orange (difference only seen if plot window is very large)
}
.partGPAssociateSpecialRows <- function(
### associate each row with the previous and next row from a subset of special rows
iRowsSpecial ##<< ordered unique integer vector specifying the rows for which some special data is available
,nRec ##<< integer scalar of number of rows in the full data.frame
){
##details<<
## When only for a subset of rows some more data available, i.e parameter estimates for
## only a subset of rows, this function creates
## columns that refer to the previous and next row that are in the subset.
## E.g. if some more data is available for rows 3 and 7, then rows 4:6 will indicate
## \code{iBefore=3, iAfter=7}.
##value<< a dataframe with index of previous and next rows inside the subset
ans <- data.frame(
iRec=1:nRec ##<< the original row number
, iSpecialBefore=NA_integer_ ##<< index within \code{iRowsSpecial}
, iSpecialAfter=NA_integer_ ##<< index within \code{iRowsSpecial}
, iBefore=NA_integer_ ##<< index of the previous special row
, iAfter=NA_integer_ ##<< index of the next special row
, wBefore=NA_real_ ##<< weight of the previous, inverse of the distance in records
, wAfter=NA_real_) ##<< weight of the next, inverse of the distance in records
##details<<
## The subset rows inside the subset refer both (before and after) to the same subset rows, with weights 1/2
nRecS <- length(iRowsSpecial)
if( 0 == nRecS ) stop("cannot associate special rows, if length of argument iRowsSpecial is zero.")
ans[iRowsSpecial,"iSpecialBefore"] <- ans[iRowsSpecial,"iSpecialAfter"] <- seq_along(iRowsSpecial)
ans[iRowsSpecial,"iBefore"] <- ans[iRowsSpecial,"iAfter"] <- iRowsSpecial
ans[iRowsSpecial,c("wBefore","wAfter")] <- 0.5
#iS <- 2L
for( iS in 1:nRecS ){
currRec <- iRowsSpecial[iS]
# before and after last special row will be treated afterwards
prevRec <- if( iS==1L) currRec else iRowsSpecial[iS-1L]
nextRec <- if( iS==nRecS) currRec else iRowsSpecial[iS+1L]
#c(prevRec,currRec,nextRec)
##details<<
## The weight is inversely proportional to the distance in rows
## The two weights wBefore and wAfter always sum to 1
distPrev <- currRec-prevRec
if( distPrev > 1L){
ans[(prevRec+1L):(currRec-1L),"iSpecialAfter"] <- iS
ans[(prevRec+1L):(currRec-1L),"iAfter"] <- currRec
ans[(prevRec+1L):(currRec-1L),"wAfter"] <- (1:(distPrev-1))/distPrev
}
distNext <- nextRec-currRec
if( distNext > 1L){
ans[(currRec+1L):(nextRec-1L),"iSpecialBefore"] <- iS
ans[(currRec+1L):(nextRec-1L),"iBefore"] <- currRec
ans[(currRec+1L):(nextRec-1L),"wBefore"] <- ((distNext-1):1)/distNext
}
}
##details<<
## the rows before the first subset row refer bot (after and before) to the first subset row with weights 1/2
## similar the rows after the last subset row refer to the last subset row
ans[1:iRowsSpecial[1],c("iSpecialBefore","iSpecialAfter")] <- 1L
ans[1:iRowsSpecial[1],c("iBefore","iAfter")] <- iRowsSpecial[1L]
ans[1:iRowsSpecial[1],c("wBefore","wAfter")] <- 0.5
ans[iRowsSpecial[nRecS]:nrow(ans),c("iSpecialBefore","iSpecialAfter")] <- nRecS
ans[iRowsSpecial[nRecS]:nrow(ans),c("iBefore","iAfter")] <- iRowsSpecial[nRecS]
ans[iRowsSpecial[nRecS]:nrow(ans),c("wBefore","wAfter")] <- 0.5
ans
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.