R/PartitioningLasslop10.R
In bgctw/REddyProc: Post Processing of (Half-)Hourly Eddy-Covariance Measurements
Defines functions
.computeRecoNight
replaceMissingSdByPercentage
.partGPAssociateSpecialRows
.tmp.f
.tmp.f
partGLInterpolateFluxes
.bootStrapLRCFit
partGLFitLRCOneWindow
partGLFitLRCWindows
partGLExtractStandardData
partGLControlLasslopCompatible
partGLControl
partitionNEEGL
Documented in
partGLControl
partGLControlLasslopCompatible
partGLExtractStandardData
partitionNEEGL
#' @export
partitionNEEGL <- function(
### Partition NEE fluxes into GP and Reco using the daytime method.
ds ##<< dataset with all the specified input columns
## and full days in equidistant times
, NEEVar = if (!missing(NEEVar.s)) NEEVar.s else paste0('NEE', suffixDash, '_f') ##<< Variable of NEE
, TempVar = if (!missing(TempVar.s)) TempVar.s else 'Tair_f' ##<< Filled air or soil temperature variable (degC)
, VPDVar = if (!missing(VPDVar.s)) VPDVar.s else 'VPD_f' ##<< Filled Vapor Pressure Deficit - VPD - (hPa)
, RadVar = if (!missing(RadVar.s)) RadVar.s else 'Rg_f' ##<< Filled radiation variable
, suffix = if (!missing(Suffix.s)) Suffix.s else "" ##<< string inserted into column names before
## identifier for NEE column defaults
## (see \code{\link{sEddyProc_sMDSGapFillAfterUstar}}).
, NEEVar.s ##<< deprecated
, TempVar.s ##<< deprecated
, VPDVar.s ##<< deprecated
, RadVar.s ##<< deprecated
, Suffix.s ##<< deprecated
## identifier for NEE column defaults
## (see \code{\link{sEddyProc_sMDSGapFillAfterUstar}}).
, ... ##<< further arguments to
## \code{\link{partGLExtractStandardData}}, such as \code{PotRadVar}
, 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)
, lrcFitter = RectangularLRCFitter() ##<< R5 class instance
## responsible for fitting the light response curve.
## Current possibilities are \code{RectangularLRCFitter()},
## \code{NonrectangularLRCFitter()},
## and \code{LogisticSigmoidLRCFitter()}.
) {
varNamesDepr <- c(
"NEEVar.s","TempVar.s","VPDVar.s","RadVar.s","Suffix.s")
varNamesNew <- c(
"NEEVar","TempVar","VPDVar","RadVar","suffix")
iDepr = which(!c(
missing(NEEVar.s),missing(TempVar.s),missing(VPDVar.s)
,missing(RadVar.s),missing(Suffix.s)))
if (length(iDepr)) warning(
"Argument names ",varNamesDepr[iDepr]," have been deprecated."
," Please, use instead ", varNamesNew[iDepr])
##details<<
## Daytime-based partitioning of measured net ecosystem fluxes into
## gross primary production (GPP) and ecosystem respiration (Reco)
##
## The fit to the light-response-curve is done by default using the Rectangular
## hyperbolic function, as in Lasslop et al. (2010)
## Alternative fittings can be used by providing the corresponding subclass of
## \code{\link{LightResponseCurveFitter-class}} to \code{lrcFitter} argument.
## (see \code{\link{LightResponseCurveFitter_predictGPP}})
#
##author<< 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
suffixDash <- paste( (if (fCheckValString(suffix)) "_" else "")
, suffix, sep = "") # used to compute default NEEVar
if (isVerbose) message('Start daytime flux partitioning for variable '
, NEEVar, ' with temperature ', TempVar, '.')
dsR <- partGLExtractStandardData(ds, NEEVar = NEEVar, TempVar = TempVar
, VPDVar = VPDVar, RadVar = RadVar , suffix = suffix
, ..., controlGLPart = controlGLPart)
##value<<
## \item{Reco_DT_<suffix>}{predicted ecosystem respiration: mumol CO2/m2/s}
## \item{GPP_DT_<suffix>}{predicted gross primary production mumol CO2/m2/s}
## \item{<LRC>}{Further light response curve (LRC) parameters and
## their standard deviation depend on the used LRC
## (e.g. for the non-rectangular LRC
## see \code{\link{NonrectangularLRCFitter_getParameterNames}}).
## They are estimated for windows and are reported with the first record
## of the window}
dsAns0 <- tibble(
FP_VARnight = rep(NA_real_, nrow(ds)) ##<< NEE filtered
## for nighttime records (others NA)
, FP_VARday = NA_real_ ##<< NEE filtered for daytime
## records (others NA)
, NEW_FP_Temp = dsR$Temp ##<< temperature after filtering for quality
## flag degree Celsius
, NEW_FP_VPD = dsR$VPD ##<< vapour pressure deficit after filtering
## for quality flag, hPa
, FP_RRef_Night = NA_real_ ##<< basal respiration estimated
## from nighttime (W / m2)
, 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
, FP_errorcode = NA_integer_ ##<< information why LRC-fit was not
## successful or was rejected, see result of
## \code{\link{LightResponseCurveFitter_fitLRC}}
, FP_GPP2000 = NA_real_ ##<< predicted GPP at VPD = 0 and PAR = 2000:
## a surrogate for maximum photosynthetic capacity
, FP_OPT_VPD = vector(mode = "list", length = nrow(ds)) ##<< list object
## of fitting results including iOpt and covParms
, FP_OPT_NoVPD = vector(mode = "list", length = nrow(ds)) ##<< same as
## FP_OPT_VPD holding optimization results with fit
## neglecting the VPD effect
)
#
# append LRC parameter result columns
lrcParNames <- lrcFitter$getParameterNames()
lrcParNames <- c(lrcParNames, paste0(lrcParNames, "_sd"))
FP_lrcParNames <- paste0("FP_", lrcParNames)
tmp <- matrix(NA_real_, nrow = nrow(ds), ncol = length(FP_lrcParNames)
, dimnames = list(NULL, FP_lrcParNames) )
dsAns <- cbind(dsAns0, tmp)
dsAns$FP_VARnight <- ifelse(dsR$isNight, dsR$NEE, NA)
attr(dsAns$FP_VARnight, 'varnames') <- paste(attr(dsR$NEE, 'varnames')
, '_night', sep = '')
attr(dsAns$FP_VARnight, 'units') <- attr(dsR$NEE, 'units')
dsAns$FP_VARday <- ifelse(dsR$isDay, dsR$NEE, NA)
attr(dsAns$FP_VARday, 'varnames') <- paste(attr(dsR$NEE, 'varnames')
, '_day', sep = '')
attr(dsAns$FP_VARday, 'units') <- attr(dsR$NEE, 'units')
#
#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")
#
##seealso<< \code{partGLFitNightTimeTRespSens}
isUsingFixedTempSens <- length(controlGLPart$fixedTempSens) &&
all(is.finite(controlGLPart$fixedTempSens$E0)) &&
all(is.finite(controlGLPart$fixedTempSens$sdE0))
dsTempSens <- if (!isUsingFixedTempSens ) {
dsTempSens <- partGLFitNightTimeTRespSens(dsR
, nRecInDay = nRecInDay
, controlGLPart = controlGLPart
, isVerbose = isVerbose
)
} else {
dsTempSens <- controlGLPart$fixedTempSens
nWindow <- length(getStartRecsOfWindows(nrow(dsR), nRecInDay = nRecInDay))
# if one row is given, use it for all windows
if (nrow(dsTempSens) == 1L) dsTempSens <- do.call(rbind,
lapply(1:nWindow, function(i) dsTempSens))
if (nrow(dsTempSens) != nWindow) stop(
"when providing controlGLPart$fixedTempSens then it must be 1 row"
,"or rows must macht the number of windows (", nWindow, "), but was "
, nrow(dsTempSens))
# if no RRef is given, estimate from nighttime data
if (!length(dsTempSens$RRef) || all(is.na(dsTempSens$RRef)) ) {
resRef15 <- simplifyApplyWindows(applyWindows(dsR
, partGLFitNightRespRefOneWindow
, nRecInDay = nRecInDay
, E0Win = dsTempSens
, controlGLPart = controlGLPart
, isVerbose = isVerbose
))
dsTempSens$RRef <- resRef15$RRef
}
# fill missing RRef by previous window
dsTempSens$RRef <- fillNAForward(dsTempSens$RRef, firstValue =
# on NA at the start of the series, take the first occuring
# finite value
dsTempSens$RRef[which(is.finite(dsTempSens$RRef))[1] ] )
dsTempSens
}
##seealso<< \code{partGLFitLRCWindows}
resParms <- resParmsWithVPD <- partGLFitLRCWindows(dsR
, nRecInDay = nRecInDay
, dsTempSens = dsTempSens
, controlGLPart = controlGLPart
, lrcFitter = lrcFitter
, isVerbose = isVerbose
)
# if no windows was fitted, return error. Else error on missing columns
# (no parameter columns returned)
if (sum(resParms$convergence == 0) == 0) stop(
"could not fit a single window in dayTime partitioning. Check the data.")
# 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$iCentralRec, c("FP_RRef_Night", "FP_qc", "FP_errorcode"
, FP_lrcParNames, "FP_GPP2000", "FP_OPT_VPD")] <-
resParms[, c("RRef_night", "parms_out_range", "convergence"
, lrcParNames, "GPP2000", "resOpt")]
##seealso<< \code{partGLInterpolateFluxes}
dsAnsFluxes <- partGLInterpolateFluxes(
#dsR$Rg
ds[[RadVar]]
#, dsAns$NEW_FP_VPD, dsAns$NEW_FP_Temp
, ds[[VPDVar]]
, ds[[TempVar]] # prediction using non-filtered, i.e. gap-filled
, resParms
, controlGLPart = controlGLPart
, lrcFitter = lrcFitter
, isVerbose = isVerbose
)
if (!controlGLPart$isNeglectVPDEffect &&
controlGLPart$isRefitMissingVPDWithNeglectVPDEffect) {
##details<<
## While the extrapolation uses filled data, the parameter optimization
## may use only measured data, i.e. with specified quality flag.
## Even with using filled VPD, there may be large gaps that have not been
## filled.
## With the common case where VPD is missing for fitting the LRC, by default
## (with \code{controlGLPart$isRefitMissingVPDWithNeglectVPDEffect = TRUE})
## is to redo the estimation of LRC parameters with neglecting the VPD-effect.
## Next, in the predictions (rows) with missing VPD are then replaced
## with predictions
## based on LRC-fits that neglected the VPD effect.
iNAVPD <- which(is.na(ds[[VPDVar]] & is.na(dsAnsFluxes$GPP)))
if (length(iNAVPD)) {
message(" could not predict GPP in ", length(iNAVPD)
, " cases due to missing VPD.")
message(" Therefore refitting LightResponseCurve with option "
,"isNeglectVPDEffect = TRUE")
ctrlNeglectVPD <- within(controlGLPart, isNeglectVPDEffect <- TRUE)
resParms <- resParmsWithoutVPD <- partGLFitLRCWindows(dsR
, nRecInDay = nRecInDay
, dsTempSens = dsTempSens
, controlGLPart = ctrlNeglectVPD
, lrcFitter = lrcFitter
)
dsAns$FP_OPT_NoVPD[resParmsWithoutVPD$iCentralRec] <-
resParmsWithoutVPD$resOpt
dsAnsFluxes2 <- partGLInterpolateFluxes(dsR$Rg
#, dsAns$NEW_FP_VPD, dsAns$NEW_FP_Temp
, ds[[VPDVar]]
, ds[[TempVar]] # do prediction also using gap-Filled values
, resParmsWithoutVPD
, controlGLPart = ctrlNeglectVPD
, lrcFitter = lrcFitter
)
dsAnsFluxes[iNAVPD, ] <- dsAnsFluxes2[iNAVPD, ]
}
}
if (isTRUE(controlGLPart$useNightimeBasalRespiration))
dsAnsFluxes <- .computeRecoNight(dsAnsFluxes, dsR, dsTempSens)
nNAGPP <- sum(is.na(dsAnsFluxes$GPP))
if (nNAGPP) warning("could not predict GPP in ", nNAGPP, " cases.")
#
dsAns$FP_dRecPar <- dsAnsFluxes$dRecNextEstimate
# copy quality flag from parameter row
# first copy from iCentralRec to colNameAssoc
# next transfer by dRecNextEstimate to each row
iFiniteMeanRec <- which(is.finite(resParms[[colNameAssoc]]))
dsAns$FP_qc[resParms[[colNameAssoc]][iFiniteMeanRec] ] <-
dsAns$FP_qc[resParms$iCentralRec[iFiniteMeanRec] ]
dsAns$FP_qc <- dsAns$FP_qc[1:nrow(dsAns) + dsAns$FP_dRecPar]
#set quality flag to 2 where next parameter estimate is more than 14 days away
dsAns$FP_qc[abs(dsAns$FP_dRecPar) > (14 * nRecInDay)] <- 2L
#dsAns[is.finite(dsAns$FP_beta), ]
#
RecoDTVar.s <- paste0('Reco_DT', suffixDash)
GPPDTVar.s <- paste0('GPP_DT', suffixDash)
RecoDTSdVar.s <- paste0(RecoDTVar.s, "_SD")
GPPDTSdVar.s <- paste0(GPPDTVar.s, "_SD")
dsAns[[RecoDTVar.s]] <- dsAnsFluxes$Reco
attr(dsAns[[RecoDTVar.s]], 'varnames') <- RecoDTVar.s
attr(dsAns[[RecoDTVar.s]], 'units') <- attr(dsR$NEE, 'units')
dsAns[[GPPDTVar.s]] <- dsAnsFluxes$GPP
attr(dsAns[[GPPDTVar.s]], 'varnames') <- GPPDTVar.s
attr(dsAns[[GPPDTVar.s]], 'units') <- attr(dsR$NEE, 'units')
if (controlGLPart$isSdPredComputed) {
dsAns[[RecoDTSdVar.s]] <- dsAnsFluxes$sdReco
attr(dsAns[[RecoDTSdVar.s]], 'varnames') <- RecoDTSdVar.s
attr(dsAns[[RecoDTSdVar.s]], 'units') <- attr(dsR$NEE, 'units')
dsAns[[GPPDTSdVar.s]] <- dsAnsFluxes$sdGPP
attr(dsAns[[GPPDTSdVar.s]], 'varnames') <- GPPDTSdVar.s
attr(dsAns[[GPPDTSdVar.s]], 'units') <- attr(dsR$NEE, '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)
}
#' @export
partGLControl <- function(
### Default list of parameters for Lasslop 2010 daytime flux partitioning
LRCFitConvergenceTolerance = 1e-3 ##<< convergence criterion for rectangular
## light response curve 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.
, nLRCFitConvergenceTolerance = 1e-3 ##<< convergence criterion for
## nonrectangular light response curve fit.
## Here its a factor of machine tolerance.
, 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 = TRUE ##<< 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
, isUsingLasslopQualityConstraints = FALSE ##<< set to TRUE to avoid
## quality constraints additional to Lasslop 2010
, 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 (radiation, temperature, 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.
, fixedTRefAtNightTime = NA ##<< if a finite value (degree Centigrade)
## is given, it is used instead of median data temperature as reference
## temperature in estimation of temperature sensitivity from night data
, isExtendTRefWindow = TRUE ##<< set to FALSE to avoid successively
## extending the night-time window
## in order to estimate a temperature sensitivity where previous estimates
## failed
, 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
, isNeglectPotRadForNight = FALSE ##<< set to TRUE to not use potential
## radiation in determining night-time data.
, NRHRfunction = FALSE ##<< deprecated: Flag if TRUE use the NRHRF
## for partitioning; Now use \code{lrcFitter = NonrectangularLRCFitter()}
, isNeglectVPDEffect = FALSE ##<< set to TRUE to avoid using VPD in the
## computations. This may help when VPD is rarely measured.
, isRefitMissingVPDWithNeglectVPDEffect = TRUE ##<< set to FALSE to avoid
## repeating estimation
## with \code{isNeglectVPDEffect = TRUE} trying to predict when VPD
## is missing
, fixedTempSens = data.frame( ##<< data.frame
## of one row or nRow = nWindow
## corresponding to return value of \code{partGLFitNightTimeTRespSens}
## While column \code{RRef} is used only as a prior and initial value for
## the daytime-fitting and can be NA,
## \code{E0} is used as given temperature sensitivity and varied according
## to \code{sdE0} in the bootstrap.
E0 = NA_real_, sdE0 = NA_real_, RRef = NA_real_)
, replaceMissingSdNEEParms = c(perc = 0.2, minSd = 0.7) ##<< parameters for
## replacing missing standard deviation of NEE.
## see \code{replaceMissingSdByPercentage}.
## Default sets missing uncertainty to 20% of NEE but at least 0.7
## flux-units (usually mumol CO2 / m2 / s).
## Specify c(NA, NA) to avoid replacing missings in standard deviation of
## NEE and to omit those records from LRC fit.
, neglectNEEUncertaintyOnMissing = FALSE ##<< If set to TRUE: if there are
## records with missing uncertainty of NEE inside one window,
## set all uncertainties to 1.
## This overrules option replaceMissingSdNEEParms.
, minPropSaturation = NA ##<< quality criterion for sufficient data
## in window. If GPP prediction of highest PAR of window is less than
## minPropSaturation * (GPP at light-saturation, i.e. beta)
## this indicates that PAR is not sufficiently high to constrain the
## shape of the LRC
, useNightimeBasalRespiration = FALSE ##<< set to TRUE to estimate
## nighttime respiration based on basal respiration estimated on
## nighttime data instead of basal respiration estimated from daytime
## data. This implements the modified daytime method from
## Keenan 2019 (doi:10.1038/s41559-019-0809-2)
) {
##author<< TW
##seealso<< \code{\link{partitionNEEGL}}
##description<<
## For highest compatibility to the pvWave code of G.Lasslop
## (used by first BGC-online tool)
## see function \code{\link{partGLControlLasslopCompatible}}.
if (NRHRfunction) stop("option 'NRHRfunction' is deprecated."
," Use instead in partitionNEEGL argument:"
,"lrcFitter = NonrectangularLRCFitter()")
if (isTRUE(neglectNEEUncertaintyOnMissing) )
replaceMissingSdNEEParms <- c(NA, NA)
ctrl <- list(
LRCFitConvergenceTolerance = LRCFitConvergenceTolerance
, nLRCFitConvergenceTolerance = nLRCFitConvergenceTolerance
, nBootUncertainty = nBootUncertainty
, minNRecInDayWindow = minNRecInDayWindow
, isAssociateParmsToMeanOfValids = isAssociateParmsToMeanOfValids
, isLasslopPriorsApplied = isLasslopPriorsApplied
, isUsingLasslopQualityConstraints = isUsingLasslopQualityConstraints
, isSdPredComputed = isSdPredComputed
, isFilterMeteoQualityFlag = isFilterMeteoQualityFlag
, isBoundLowerNEEUncertainty = isBoundLowerNEEUncertainty
, fixedTRefAtNightTime = fixedTRefAtNightTime
, isExtendTRefWindow = isExtendTRefWindow
, smoothTempSensEstimateAcrossTime = smoothTempSensEstimateAcrossTime
, isNeglectPotRadForNight = isNeglectPotRadForNight
, isNeglectVPDEffect = isNeglectVPDEffect
, isRefitMissingVPDWithNeglectVPDEffect = isRefitMissingVPDWithNeglectVPDEffect
, fixedTempSens = fixedTempSens
, replaceMissingSdNEEParms = replaceMissingSdNEEParms
, neglectNEEUncertaintyOnMissing = neglectNEEUncertaintyOnMissing
, minPropSaturation = minPropSaturation
, useNightimeBasalRespiration = useNightimeBasalRespiration
)
##value<< list with entries of given arguments.
ctrl
}
attr(partGLControl, "ex") <- function() {
partGLControl(nBootUncertainty = 40L)
}
#' @export
partGLControlLasslopCompatible <- function(
### Daytime flux partitioning parms compatible with with the pvWave
nBootUncertainty = 0L ##<< 0: Derive uncertainty from
## curvature of a single fit, neglecting the uncertainty of previously
## estimated temperature sensitivity, E0
, minNRecInDayWindow = 10L ##<< Minimum number of 10 valid records
## for regression in a single window
, isAssociateParmsToMeanOfValids = FALSE ##<< associate parameters to
## the first record of the window for interpolation instead of mean across
## valid records inside a window
, isLasslopPriorsApplied = TRUE ##<< Apply fixed Lasslop priors
## in LRC fitting.
, isUsingLasslopQualityConstraints = TRUE ##<< avoid quality constraints
## additional to the ones in Lasslop 2010
, isBoundLowerNEEUncertainty = FALSE ##<< FALSE: avoid adjustment of very
## low uncertainties before
## day-Time fitting that avoids the high leverage those records with
## unreasonable low uncertainty.
, fixedTRefAtNightTime = 15 ##<< use fixed (degree Centigrade)
## temperature sensitivity
## instead of median data temperature as reference temperature in
## estimation of temperature sensitivity from night data
, isExtendTRefWindow = FALSE ##<< avoid successively extending the
##night-time window
## in order to estimate a temperature sensitivity where previous
## estimates failed
, smoothTempSensEstimateAcrossTime = FALSE ##<< FALSE: use independent
## estimates of temperature
## sensitivity on each windows instead of a vector of E0 that is
## smoothed over time
, isRefitMissingVPDWithNeglectVPDEffect = FALSE ##<< FALSE: avoid
## repeating estimation with \code{isNeglectVPDEffect = TRUE}
, minPropSaturation = NA ##<< NA: avoid quality constraint of
## sufficient saturation in data
## This option is overruled, i.e. not considered, if option
## isUsingLasslopQualityConstraints = TRUE.
, isNeglectVPDEffect = FALSE ##<< FALSE: do not neglect VPD effect
, replaceMissingSdNEEParms = c(NA, NA) ##<< do not replace missing NEE,
## but see option
, neglectNEEUncertaintyOnMissing = TRUE ##<< if there are records with
## missing uncertainty of NEE inside one window, set all sdNEE to 1.
## This overrules option replaceMissingSdNEEParms.
, ... ##<< further arguments to \code{\link{partGLControl}}
) {
##author<< TW
##seealso<< \code{\link{partGLControl}}
partGLControl(
nBootUncertainty = nBootUncertainty
, minNRecInDayWindow = minNRecInDayWindow
, isAssociateParmsToMeanOfValids = isAssociateParmsToMeanOfValids
, isLasslopPriorsApplied = isLasslopPriorsApplied
, isUsingLasslopQualityConstraints = isUsingLasslopQualityConstraints
, isBoundLowerNEEUncertainty = isBoundLowerNEEUncertainty
, fixedTRefAtNightTime = fixedTRefAtNightTime
, isExtendTRefWindow = isExtendTRefWindow
, smoothTempSensEstimateAcrossTime = smoothTempSensEstimateAcrossTime
, isNeglectVPDEffect = isNeglectVPDEffect
, isRefitMissingVPDWithNeglectVPDEffect = isRefitMissingVPDWithNeglectVPDEffect
, replaceMissingSdNEEParms = replaceMissingSdNEEParms
, neglectNEEUncertaintyOnMissing = neglectNEEUncertaintyOnMissing
, minPropSaturation = minPropSaturation
, ...
)
}
attr(partGLControlLasslopCompatible, "ex") <- function() {
partGLControlLasslopCompatible()
}
#' @export
partGLExtractStandardData <- function(
### Relevant columns from original input with defined names
ds ##<< dataset with all the specified input columns and
## full days in equidistant times
, NEEVar = paste0('NEE', suffixDash, '_f') ##<< Variable of NEE
, QFNEEVar = if (!missing(QFNEEVar.s)) QFNEEVar.s else paste0('NEE', suffixDash, '_fqc') ##<< Quality
## flag of variable
, QFNEEValue = if (!missing(QFNEEValue.n)) QFNEEValue.n else 0 ##<< Value of quality flag for
## _good_ (original) data
, NEESdVar = if (!missing(NEESdVar.s)) NEESdVar.s else paste0('NEE', suffixDash, '_fsd') ##<< Variable of
## standard deviation of net ecosystem fluxes
, TempVar = paste0('Tair_f') ##<< Filled air or soil
## temperature variable (degC)
, QFTempVar = if (!missing(QFTempVar.s)) QFTempVar.s else paste0('Tair_fqc') ##<< Quality flag of
## filled temperature variable
, QFTempValue = if (!missing(QFTempValue.n)) QFTempValue.n else 0 ##<< Value of temperature quality flag
##for _good_ (original) data
, VPDVar = if (!missing(VPDVar.s)) VPDVar.s else paste0('VPD_f') ##<< Filled Vapor Pressure Deficit, VPD (hPa)
, QFVPDVar = if (!missing(QFVPDVar.s)) QFVPDVar.s else paste0('VPD_fqc') ##<< Quality flag of filled VPD variable
, QFVPDValue = if (!missing(QFVPDValue.n)) QFVPDValue.n else 0 ##<< Value of VPD quality flag for
## _good_ (original) data
, RadVar = if (!missing(RadVar.s)) RadVar.s else 'Rg_f' ##<< Filled radiation variable
, QFRadVar = if (!missing(QFRadVar.s)) QFRadVar.s else paste0('Rg_fqc') ##<< Quality flag of filled radiation variable
, QFRadValue = if (!missing(QFRadValue.n)) QFRadValue.n else 0 ##<< Value of radiation quality flag for
## _good_ (original) data
, PotRadVar = if (!missing(PotRadVar.s)) PotRadVar.s else "PotRad_NEW" ##<< Variable name of potential rad. (W / m2)
, suffix = if (!missing(Suffix.s)) Suffix.s else "" ##<< string inserted into column names before
## identifier for NEE column defaults
## (see \code{\link{sEddyProc_sMDSGapFillAfterUstar}}).
, NEEVar.s ##<< deprecated
, QFNEEVar.s ##<< deprecated
, QFNEEValue.n ##<< deprecated
, NEESdVar.s ##<< deprecated
, TempVar.s ##<< deprecated
, QFTempVar.s ##<< deprecated
, QFTempValue.n ##<< deprecated
, VPDVar.s ##<< deprecated
, QFVPDVar.s ##<< deprecated
, QFVPDValue.n ##<< deprecated
, RadVar.s ##<< deprecated
, QFRadVar.s ##<< deprecated
, QFRadValue.n ##<< deprecated
, PotRadVar.s ##<< deprecated
, Suffix.s ##<< deprecated
, controlGLPart = partGLControl() ##<< further default parameters,
## see \code{\link{partGLControl}}
) {
if (!missing(NEEVar.s)) NEEVar <- NEEVar.s # in default, lines too wide
if (!missing(TempVar.s)) TempVar <- TempVar.s
varNamesDepr <- c(
"NEEVar.s","TempVar.s","VPDVar.s","RadVar.s","Suffix.s"
,"QFNEEVar.s" ,"QFNEEValue.n" ,"QFTempVar.s" ,"QFTempValue.n"
,"QFVPDVar.s" ,"QFVPDValue.n" ,"QFRadVar.s" ,"QFRadValue.n" ,"PotRadVar.s"
)
varNamesNew <- c(
"NEEVar","TempVar","VPDVar","RadVar","suffix"
,"QFNEEVar" ,"QFNEEValue" ,"QFTempVar" ,"QFTempValue"
,"QFVPDVar" ,"QFVPDValue" ,"QFRadVar" ,"QFRadValue" ,"PotRadVar"
)
iDepr = which(!c(
missing(NEEVar.s),missing(TempVar.s),missing(VPDVar.s),missing(RadVar.s)
,missing(Suffix.s),missing(QFNEEVar.s),missing(QFNEEValue.n)
,missing(QFTempVar.s),missing(QFTempValue.n),missing(QFVPDVar.s)
,missing(QFVPDValue.n),missing(QFRadVar.s),missing(QFRadValue.n)
,missing(PotRadVar.s)))
if (length(iDepr)) warning(
"Argument names ",varNamesDepr[iDepr]," have been deprecated."
," Please, use instead ", varNamesNew[iDepr])
# Check if specified columns exist in sDATA or sTEMP and if numeric
# and plausible. Then apply quality flag
suffixDash <- paste( (if (fCheckValString(suffix)) "_" else "")
, suffix, sep = "")
# Check if specified columns exist in sDATA or sTEMP and if numeric
# and plausible. Then apply quality flag
fCheckColNames(ds, c(NEEVar, QFNEEVar, TempVar, QFTempVar, RadVar, VPDVar
, QFRadVar, PotRadVar, NEESdVar), 'sGLFluxPartition')
fCheckColNum(ds, c(NEEVar, QFNEEVar, TempVar, QFTempVar, RadVar, VPDVar
, QFRadVar, PotRadVar, NEESdVar), 'sGLFluxPartition')
fCheckColPlausibility(ds, c(NEEVar, QFNEEVar, TempVar, QFTempVar, VPDVar
, RadVar, QFRadVar, PotRadVar), 'sGLFluxPartition')
NEEFiltered <- fSetQF(ds, NEEVar, QFNEEVar, QFNEEValue, 'sGLFluxPartition')
#
# Apply quality flag for temperature and VPD
# TODO: docu meteo filter, standard FALSE
NEW_FP_Temp <- if (isTRUE(controlGLPart$isFilterMeteoQualityFlag) ) fSetQF(
ds, TempVar, QFTempVar, QFTempValue, 'partGLExtractStandardData') else
ds[[TempVar]]
NEW_FP_VPD <- if (isTRUE(controlGLPart$isFilterMeteoQualityFlag) )
fSetQF(ds, VPDVar, QFVPDVar, QFVPDValue, 'partGLExtractStandardData') else
ds[[VPDVar]]
NEW_FP_Rg <- if (isTRUE(controlGLPart$isFilterMeteoQualityFlag) )
fSetQF(ds, RadVar, QFRadVar, QFRadValue, 'partGLExtractStandardData') else
ds[[RadVar]]
#
# 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
#! New code: PotRad in sGLFluxPartition: Slightly different subset than
#PV-Wave due to time zone correction (avoids timezone offset between Rg and PotRad)
isPotRadZero <- if (isTRUE(controlGLPart$isNeglectPotRadForNight) )
TRUE else (ds[[PotRadVar]] <= 0)
isNotPotRadZero <- if (isTRUE(controlGLPart$isNeglectPotRadForNight) )
TRUE else (ds[[PotRadVar]] > 0)
isNight <- (ds[[RadVar]] <= 4 & isPotRadZero)
# 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, twutz does not understand this comment
isDay = (ds[[RadVar]] > 4 & isNotPotRadZero)
#
##value<< a data.frame with columns
dsR <- data.frame(
sDateTime = ds[[1]] ##<< first column of \code{ds},
## usually the time stamp
## not used, but usually first column is a DateTime is kept
## for aiding debug
, NEE = NEEFiltered ##<< NEE filtered for quality flay
, sdNEE = ds[[NEESdVar]] ##<< standard deviation of NEE
## with missing values replaced
, Temp = NEW_FP_Temp ##<< Temperature, quality filtered
## if isTRUE(controlGLPart$isFilterMeteoQualityFlag)
, VPD = NEW_FP_VPD ##<< Water pressure deficit, quality
## filtered if isTRUE(controlGLPart$isFilterMeteoQualityFlag)
, Rg = NEW_FP_Rg ##<< Incoming radiation
, isDay = isDay ##<< Flag that is true for daytime records
, isNight = isNight ##<< Flag that is true for nighttime records
)
##details<<
## The LRC fit usually weights NEE records by its uncertainty. In order
## to also use
## records with missing \code{NEESdVar}, uncertainty of the missing values
## is by default set
## to a conservatively high value, parameterized by
## \code{controlGLPart$replaceMissingSdNEEParms)}.
## Controlled by argument \code{replaceMissingSdNEEParms} in
## \code{\link{partGLControl}}, but overruled
## by argument \code{neglectNEEUncertaintyOnMissing}.
if (!controlGLPart$neglectNEEUncertaintyOnMissing)
dsR$sdNEE <- replaceMissingSdByPercentage(dsR$sdNEE, dsR$NEE
, controlGLPart$replaceMissingSdNEEParms[1]
, controlGLPart$replaceMissingSdNEEParms[2] )
dsR
}
partGLFitLRCWindows <- function(
### Estimate successive Rectangular Hyperbolic Light Response Curve parameters
ds ##<< data.frame with numeric columns NEE, sdNEE,
## Temp (degC), VPD, Rg, and logical columns isNight and isDay
, winSizeDays = 4L ##<< Window size in days for daytime fits
, strideInDays = 2L ##<< step in days for shifting the windows
, nRecInDay = 48L ##<< number of records within one day
## (for half-hourly data its 48)
, dsTempSens ##<< data.frame that reports for each window
## temperature sensitivity parameters E0 and RRef
, isVerbose = TRUE ##<< set to FALSE to suppress messages
, controlGLPart = partGLControl() ##<< list of further default parameters
, lrcFitter ##<< R5 class instance responsible for fitting the
## light response curve
) {
##seealso<< \code{\link{partGLFitLRCOneWindow}}
if (isVerbose) message(
" Estimating light response curve parameters from day time NEE "
, appendLF = FALSE)
resLRC <- applyWindows(ds, partGLFitLRCOneWindow
, prevRes = list(resOpt = list(thetaOpt = rep(NA_real_, 6L)))
, winSizeInDays = 4L
, isVerbose = isVerbose
, nRecInDay = nRecInDay
#
, E0Win = dsTempSens
, controlGLPart = controlGLPart
, lrcFitter = lrcFitter
)
lrcSummary <- lapply(resLRC$resFUN, "[[", "summary")
iNoSummary <- which(sapply(lrcSummary, length) == 0)
if (length(iNoSummary) ) {
stop("expected summary returned by all fits, but found missing summaries.")
# # put dummy NA data.frame where no fit was obtained
# dummySummary <- lrcSummary[-iNoSummary][[1]]
# dummySummary[] <- NA
# lrcSummary[iNoSummary] <- list(dummySummary)
}
resOptList <- lapply(resLRC$resFUN, "[[", "resOpt")
resOptDf <- tibble(resOpt = resOptList)
resParms <- as_tibble(cbind(resLRC$winInfo, bind_rows(lrcSummary), resOptDf))
#table(resParms$convergence)
#E0_night equals E0, but uncertainty might differ
resParms$E0_bootstrap_sd <- resParms$E0_sd # due to bootstrap,
# this may differ, save before overriding by night-time estimate
resParms$E0 <- dsTempSens$E0 # report E0 and E0_sd from nighttime
# even when LRC fit did not converge
resParms$E0_sd <- dsTempSens$sdE0
resParms$RRef_night <- dsTempSens$RRef
# summary$iMeanRec yet based on window instead of entire time, need to
# add beginning of window
resParms$iMeanRec <- resParms$iRecStart - 1L + resParms$iMeanRec
# # omit records where NULL was returned
# iWinNoFit <- which(is.na(resParms$parms_out_range) )
# if (length(iWinNoFit) ) {
# resParms <- resParms[-iWinNoFit, ]
# resParms$resOptList <- resParms$resOptList[-iWinNoFit]
# }
##value<< a tibble with a row for each window
## providing information on window size, parameters estimated
## from light response curve (LRC),
## and their standard deviation.
## Moreover, and estimated from night-time data after smoothing
## and forward-filling:
## the uncertainty of temperature sensitivity \code{E0_night_sd},
## And the respiration at reference temperature \code{RRef_night}.
## Finally it contains a list column\code{resOpt} with the optimization results.
resParms
}
partGLFitLRCOneWindow <- function(
### Estimate Rectangular Hyperbolic Light Response Curve parameters
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
, lrcFitter = RectangularLRCFitter() ##<< R5 class instance responsible
## for fitting the light response curve
) {
##author<< TW
##seealso<< \code{\link{partitionNEEGL}}
##details<<
## Estimation as in Lasslop et al., 2010 for successive periods, i.e. windows.
## Fitting is done on a subset of the data where isDay, NEE, sdNEE, Temp,
## Rg, and VPD are all non-NA
## and isDay is TRUE
isValidDayRecNoSdNEEConstraint <- !is.na(ds$isDay) & ds$isDay &
!is.na(ds$NEE) & !is.na(ds$Temp) & !is.na(ds$Rg)
##details<<
## With option \code{controlGLPart$neglectNEEUncertaintyOnMissing = TRUE}
## all sdNEE are set to 1 if there is any missing sdNEE at
## otherwise valid records
if (isTRUE(controlGLPart$neglectNEEUncertaintyOnMissing) &&
any(!is.finite(ds$sdNEE[isValidDayRecNoSdNEEConstraint])) ) {
ds$sdNEE <- 1L
}
isValidDayRecNoVPDConstraint <- isValidDayRecNoSdNEEConstraint &
!is.na(ds$sdNEE)
##details<<
## If \code{controlGLPart$isNeglectVPDEffect = TRUE}, also records
## with VPD = NA maybe valid
isValidDayRec <- if (isTRUE(controlGLPart$isNeglectVPDEffect) )
isValidDayRecNoVPDConstraint else
isValidDayRecNoVPDConstraint & !is.na(ds$VPD)
iMeanRecInDayWindow <- as.integer(round(mean(which(isValidDayRec))))
#
# if too few records or
# if no temperature-respiration relationship could be found,
# indicate no-fit, but report Window properties
E0 <- E0Win$E0[winInfo$iWindow]
# when provided fixed E0 the column sdEO does not exist
sdE0 <- if (length(E0Win$sdE0)) E0Win$sdE0[winInfo$iWindow] else NA
RRefNight <- E0Win$RRef[winInfo$iWindow]
getNAResult <- function(convergenceCode) {
list(
resOpt = NULL
, summary = data.frame(
nValidRec = sum(isValidDayRec)
, iMeanRec = iMeanRecInDayWindow
, convergence = convergenceCode
, E0 = E0
, E0_sd = sdE0
, RRefNight = RRefNight
)
, isValid = FALSE
)}
if (is.na(E0) ) return(getNAResult(1010L))
#
##details<< If there are too few records (<
## \code{controlGLPart$minNRecInDayWindow}), then
## the constraint on non-NA VPD is neglected and
## \code{controlGLPart$isNeglectVPDEffect} is set to TRUE
## If there are still too few records, an NA-result with convergence
## code 1011L is returned.
if ( (sum(isValidDayRec) < controlGLPart$minNRecInDayWindow) ) {
controlGLPart$isNeglectVPDEffect <- TRUE
isValidDayRec <- isValidDayRecNoVPDConstraint
iMeanRecInDayWindow <- as.integer(round(mean(which(isValidDayRec))))
if ( (sum(isValidDayRec) < controlGLPart$minNRecInDayWindow) )
return(getNAResult(1011L))
}
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.
#better report NA and care for it properly: if (is.na(iMeanRecInDayWindow))
# iMeanRecInDayWindow <- as.integer(nrow(ds)%/%2)
#TODO firstRecInDayWindow.i <- which(SubsetDayPeriod.b)[1]
# the rownumber of the first record inside the day window
##seealso<< \code{\link{partGLEstimateTempSensInBoundsE0Only}}
#
##seealso<< \code{\link{LightResponseCurveFitter_fitLRC}}
#lrcFitter <- RectangularLRCFitter()
resOpt <- resOpt0 <- lrcFitter$fitLRC(dsDay, E0 = E0, sdE0 = sdE0
, RRefNight = RRefNight
, controlGLPart = controlGLPart
, lastGoodParameters = prevRes$resOpt$thetaOpt
)
if (!is.finite(resOpt$thetaOpt[1]) ) {
return(getNAResult(resOpt$convergence))
}
sdTheta <- resOpt$thetaOpt; sdTheta[] <- NA
sdTheta[resOpt$iOpt] <- sqrt(diag(resOpt$covParms)[resOpt$iOpt])
#
# record valid fits results
#as.data.frame(t(resOpt$thetaOpt))
#if (as.POSIXlt(dsDay$sDateTime[1])$mday + 2L >= 11) recover()
#if (as.POSIXlt(dsDay$sDateTime[1])$mday + 2L >= 27) recover()
# save(dsDay, file = "tmp / dsDayDebug.RData")
#
# compute GPP2000 without VPD effect
GPP2000 <- lrcFitter$predictLRC(
resOpt$thetaOpt, 2000, VPD = NA, Temp = NA, fixVPD = TRUE)$GPP
ans <- list(
resOpt = resOpt
, summary = cbind(data.frame(
nValidRec = nrow(dsDay)
, iMeanRec = iMeanRecInDayWindow
, convergence = resOpt$convergence
, parms_out_range = as.integer(!all(1:5 %in% resOpt$iOpt) )
)
, as.data.frame(t(resOpt$thetaOpt))
, as.data.frame(t(structure(sdTheta, names = paste0(names(sdTheta), "_sd"))))
, GPP2000 = GPP2000
, isValid = TRUE
)
)
return(ans)
}
.bootStrapLRCFit <- function(
### Compute parameters uncertainty by bootstrap
theta0, iOpt, dsDay, sdE_0.n, parameterPrior, controlGLPart
, lrcFitter ##<< Light Response Curve R5 instance
, iPosE0 = 5L ##<< position (integer scalar) of temperature
##sensitivity in parameter vector
#should be dealt with iOpt: isNeglectVPD = FALSE
) {
##value<<
## Matrix with each row a parameter estimate on a different bootstrap sample.
## Some of the rows might be NA, if the subsampled dataset could not be fitted.
ans <- matrix(NA, nrow = controlGLPart$nBootUncertainty, ncol = length(theta0)
, dimnames = list(NULL, names(theta0)))
##details<<
## In addition to resampling the original data, also the temperature sensitivity
## is resampled
## from its uncertainty distribution but bounded in between 50 to 400.
E0r <- rnorm(controlGLPart$nBootUncertainty, theta0[5L], sdE_0.n )
E0 <- pmax(50, pmin(400, E0r))
theta <- theta0
#iBoot <- 1L
for (iBoot in c(1:controlGLPart$nBootUncertainty)) {
idx <- sample(nrow(dsDay), replace = TRUE)
dsDayB <- dsDay[idx, ]
theta[iPosE0] <- E0[iBoot]
resOptBoot <- try(lrcFitter$optimLRCOnAdjustedPrior(theta, iOpt = iOpt
, dsDay = dsDayB
, parameterPrior = parameterPrior, ctrl = controlGLPart), silent = TRUE)
if (!inherits(resOptBoot, "try-error") && resOptBoot$convergence == 0L) {
#TODO: also remove the very bad cases?
ans[iBoot, ]<- resOptBoot$theta
}else{
#recover()
}
}
ans
}
partGLInterpolateFluxes <- function(
### Interpolate Reco and GPP from two neighboring parameter sets
Rg ##<< numeric vector of photosynthetic flux density [umol / m2 / s]
## or Global Radiation
, VPD ##<< numeric vector of Vapor Pressure Deficit [hPa]
, Temp ##<< numeric vector of Temperature [degC]
, resParms ##<< data frame with results of \code{partGLFitLRCWindows}
## of fitting the light-response-curve for several windows
, controlGLPart = partGLControl() ##<< further default parameters,
## see \code{\link{partGLControl}}
, lrcFitter ##<< R5 class instance responsible for fitting the LRC curve
, isVerbose = TRUE ##<< set to FALSE to suppress messages
) {
##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
# corresponding weights
iValidWin <- which(is.finite(resParms$parms_out_range))
summaryLRC <- resParms %>% select(!"resOpt") %>% slice(iValidWin)
resOptList <- resParms$resOpt[iValidWin]
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
# they all will be based on the same valid data, if iMeanRec is equal
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 central 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 columns iBefore and iAfter, which can be used to merge the
# parameter estimate columns to each row
parNames <- lrcFitter$getParameterNames()
dsBefore <- merge(
structure(data.frame(dsAssoc$iSpecialBefore, dsAssoc$iBefore)
, names = c("iParRec", colNameAssoc))
, summaryLRC[, c(colNameAssoc, parNames)]
)
dsAfter <- merge(structure(data.frame(dsAssoc$iSpecialAfter, dsAssoc$iAfter)
, names = c("iParRec", colNameAssoc))
, summaryLRC[, c(colNameAssoc, lrcFitter$getParameterNames())])
if ( (nrow(dsBefore) != nRec) || (nrow(dsAfter) != nRec)) stop(#
"error in merging parameters to original records.")
Reco2 <- lapply(list(dsBefore, dsAfter), function(dsi) {
#twutz170316: fLloydTaylor gives unreasonable values with very
# low temperatures, hence constrain lower temperature
tmp <- fLloydTaylor(dsi$RRef, dsi$E0, pmax(-40, Temp_Kelvin)
, TRef = 273.15 + 15)
})
#dsi <- dsBefore
GPP2 <- lapply(list(dsBefore, dsAfter), function(dsi) {
theta <- as.matrix(dsi[, parNames])
tmp <- lrcFitter$predictLRC(theta, Rg, VPD
, Temp = pmax(-40, 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)
ans <- ansPred <- data.frame(
Reco = Reco
, GPP = GPP
)
if (isTRUE(controlGLPart$isSdPredComputed)) {
#dsi <- dsBefore
varPred2 <- lapply(list(dsBefore, dsAfter), function(dsi) {
theta <- as.matrix(dsi[, parNames])
grad <- lrcFitter$computeLRCGradient(theta, Rg, VPD, pmax(-40, 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] # row within sequence of valid
# parameters (at centralREc or meanRec)
# get the fitting object, TODO better document
resOpt <- resOptList[[iParRec]]
gradGPP <- grad$GPP[iRec, ]
gradReco <- grad$Reco[iRec, ]
# make sure parameter names match postions in covParms
varPred[iRec, 2L] <- varReco <- gradReco %*%
resOpt$covParms[names(gradReco), names(gradReco)] %*% gradReco
varPred[iRec, 1L] <- varGPP <- gradGPP %*%
resOpt$covParms[names(gradGPP), names(gradGPP)] %*% gradGPP
}
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() { # omit declaration for now
computeAggregatedCovariance <- function(
### Covariances between Reco and GPP predictions due same uncertain model coefficients
dsPred ##<< data.frame with predictors (Rg, VPD, Temp)
, resParms ##<< data.frame with results of \code{partGLFitLRCWindows} of
## fitting the light-response-curve for several windows
, resParmsNoVPD ##<< data.frame like resParms, but was fitted with option
## isNeglectVPDEffect = TRUE, for predicting if VPD is missing
, controlGLPart = partGLControl() ##<< further default parameters,
## see \code{\link{partGLControl}} with entry "isAssociateParmsToMeanOfValids"
, lrcFitter ##<< R5 class instance responsible for fitting the light response
## curve, with method getParameterNames()
, iAggregate = 1:nrow(dsPred) ##<< row indices about which to
## sum over, must be contiguous
) {
sumCovGPP <- 0
sumCovReco <- 0
# merge parameters to predictors
isValidWin <- is.finite(resParms$parms_out_range)
summaryLRC <- resParms[isValidWin, , drop = FALSE]
resOptList <- resParms$resOptList[isValidWin]
colNameAssoc <- if (isTRUE(controlGLPart$isAssociateParmsToMeanOfValids) )
"iMeanRec" else "iCentralRec"
dsAssoc <- .partGPAssociateSpecialRows(summaryLRC[[colNameAssoc]], nRec)
# now we have columns iBefore and iAfter, which can be used to
# merge the parameter estimate columns to each row
parNames <- lrcFitter$getParameterNames()
dsBefore <- merge(
structure(data.frame(dsAssoc$iSpecialBefore, dsAssoc$iBefore)
, names = c("iParRec", colNameAssoc))
, summaryLRC[, c(colNameAssoc, parNames)]
)
dsAfter <- merge(structure(data.frame(dsAssoc$iSpecialAfter, dsAssoc$iAfter)
, names = c("iParRec", colNameAssoc))
, summaryLRC[, c(colNameAssoc, lrcFitter$getParameterNames())])
#
# compute gradients (each entry is a list with components vectors Reco and GPP
grad2 <- lapply(list(dsBefore, dsAfter), function(dsi) {
theta <- as.matrix(dsi[, parNames])
grad <- lrcFitter$computeLRCGradient(theta, dsPred$Rg, dsPred$VPD
, pmax(-40, dsPred$Temp))
})
for (iRec in iAggregate[-1]) { # in first line only entry with i == j
iParRecBefore <- dsAssoc$iSpecialBefore[iRec]
resOptBefore <- resOptList[[iParRecBefore]]
iParRecAfter <- dsAssoc$iSpecialAfter[iRec]
resOptAfter <- resOptList[[iParRecAfter]]
gradGPPIBefore <- grad[[1]]$GPP[iRec, ]
gradRecoIBefore <- grad[[1]]$Reco[iRec, ]
gradGPPIAfter <- grad[[2]]$GPP[iRec, ]
gradRecoIAfter <- grad[[2]]$Reco[iRec, ]
for (jRec in iAggregate[1]:(iRec - 1) ) { # iter from start to one before iRec
# assume covariance larger than zero only in same interval between estimates
if ( (dsAssoc$iSpecialBefore[jRec] == iParRecBefore) &&
(dsAssoc$iSpecialAfter[jRec] == iParRecAfter)
) {
gradGPPJBefore <- grad[[1]]$GPP[jRec, ]
gradRecoJBefore <- grad[[1]]$Reco[jRec, ]
gradGPPJAfter <- grad[[2]]$GPP[jRec, ]
gradRecoJAfter <- grad[[2]]$Reco[jRec, ]
covRecoBefore <- gradRecoIBefore %*%
resOptBefore$covParms[names(gradRecoIBefore), names(gradRecoIBefore)] %*%
gradRecoJBefore
covGPPBefore <- gradGPPIBefore %*%
resOptBefore$covParms[names(gradGPPIBefore), names(gradGPPIBefore)] %*%
gradGPPJBefore
covRecoAfter <- gradRecoIAfter %*%
resOptAfter$covParms[names(gradRecoIAfter), names(gradRecoIAfter)] %*%
gradRecoJAfter
covGPPBAfter <- gradGPPIAfter %*%
resOptAfter$covParms[names(gradGPPIAfter), names(gradGPPIAfter)] %*%
gradGPPJAfter
covRecoMix <- dsAssoc$wBefore[iRec] * dsAssoc$wBefore[jRec] * covRecoBefore +
dsAssoc$wAfter[iRec] * dsAssoc$wAfter[jRec] * covRecoAfter
covGPPMix <- dsAssoc$wBefore[iRec] * dsAssoc$wBefore[jRec] * covGPPBefore +
dsAssoc$wAfter[iRec] * dsAssoc$wAfter[jRec] * covGPPAfter
sumCovReco <- sumCovReco + 2 * covRecoMix # symm, add also from transposed
sumCovGPP <- sumCovGPP + 2 * covGPPMix
} # match of iParRecBefore and iParRecAfter
} # for jRec
} # for iRec
##value<< named numeri vector with two components:
ans <- c(sumCovGPP = sumCovGPP ##<< sum_{i<>j} cov(GPP_i, GPP_j)
, sumCovReco = sumCovReco ##<< sum_{i<>j} cov(Reco_i, Reco_j)
)
}
} # .tmp.f
.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 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
}
replaceMissingSdByPercentage <- function(
### replace missing standard deviation of a measure x by a percentage of x
sdX ##<< numeric vector: with missing to be repalce
, x ##<< numeric vector of length(sdX): value form which percentage is computed
, perc = 0.2 ##<< numeric scalar: sdX = perc * x
, minSdX = 0.7 ##<< numeric scalar: minimum of sdX to be applied for low x
) {
##details<<
## If either perc or inSdX is NA then only the other criterion is applied.
## If both are NA then all missings are set to NA.
## \code{sdX[iToFill] <- pmax(minSdX, abs(x[iToFill] * perc), na.rm = TRUE)}
iToFill <- !is.finite(sdX)
sdX[iToFill] <- pmax(minSdX, abs(x[iToFill] * perc), na.rm = TRUE)
##value<< sdX with non-finite values replaced.
sdX
}
.computeRecoNight <- function(
### recompute ecosystem respiration at nighttime with nighttime reference
dsAnsFluxes ##<< data.frame with columsn Reco and sdReco
, dsR ##<< data.frame with columsn isNight, and Temp
, dsTempSens ##<< data.frame with column iCentralRec and E0
) {
##details<< if \code{partGLControl$useNightimeBasalRespiration == TRUE}
## then nighttime respiration is recomputed based on basal respiration
## inferred from nighttime data.
# Calculate the ecosystem respiration Reco by LlyodAndTaler
TRefK = 273.15 + 15; T0 = 227.13
RRef <- approx(
dsTempSens$iCentralRec, dsTempSens$RRef, 1:nrow(dsR), rule = 2)$y
E0 <- approx(dsTempSens$iCentralRec, dsTempSens$E0, 1:nrow(dsR), rule = 2)$y
dsAnsFluxes$Reco[dsR$isNight] <- respNight <- fLloydTaylor(
RRef[dsR$isNight], E0[dsR$isNight] , fConvertCtoK(dsR$Temp[dsR$isNight])
, TRef = TRefK)
##details<< With given uncertainty of temperature sensitivity E0, the
## logarithm of respraiton R is normally distributed, i.e. R is lognormally
## distributed. With low uncertainty (sigma*) < 1.2 this is well
## approximated by a normal distribution. Hence the sqrt(second momemt)
## of the lognormal is reported as sdReco
# log(R) = log(RRef) + tempFac*E0
sdE0 <- sqrt(approx(
dsTempSens$iCentralRec, dsTempSens$sdE0^2, 1:nrow(dsR), rule = 2)$y)
tempFacLloydTaylor <- 1/(TRefK - T0) - 1/(fConvertCtoK(dsR$Temp) - T0)
sigmaLogR = (abs(tempFacLloydTaylor)*sdE0)[dsR$isNight]
# formula of variance of lognormal with mu replaced by mean
# see lognrom/inst/docu/varianceBySigmaAndExpected.Rmd
dsAnsFluxes$sdReco[dsR$isNight] <- sqrt(exp(sigmaLogR^2) - 1)*respNight
##value<< \code{dsAnsfluxes} with updated night-time Reco and sdReco
dsAnsFluxes
}
bgctw/REddyProc documentation built on March 26, 2024, 11:35 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.
Defines functions .computeRecoNight replaceMissingSdByPercentage .partGPAssociateSpecialRows .tmp.f .tmp.f partGLInterpolateFluxes .bootStrapLRCFit partGLFitLRCOneWindow partGLFitLRCWindows partGLExtractStandardData partGLControlLasslopCompatible partGLControl partitionNEEGL
Documented in partGLControl partGLControlLasslopCompatible partGLExtractStandardData partitionNEEGL
#' @export
partitionNEEGL <- function(
### Partition NEE fluxes into GP and Reco using the daytime method.
ds ##<< dataset with all the specified input columns
## and full days in equidistant times
, NEEVar = if (!missing(NEEVar.s)) NEEVar.s else paste0('NEE', suffixDash, '_f') ##<< Variable of NEE
, TempVar = if (!missing(TempVar.s)) TempVar.s else 'Tair_f' ##<< Filled air or soil temperature variable (degC)
, VPDVar = if (!missing(VPDVar.s)) VPDVar.s else 'VPD_f' ##<< Filled Vapor Pressure Deficit - VPD - (hPa)
, RadVar = if (!missing(RadVar.s)) RadVar.s else 'Rg_f' ##<< Filled radiation variable
, suffix = if (!missing(Suffix.s)) Suffix.s else "" ##<< string inserted into column names before
## identifier for NEE column defaults
## (see \code{\link{sEddyProc_sMDSGapFillAfterUstar}}).
, NEEVar.s ##<< deprecated
, TempVar.s ##<< deprecated
, VPDVar.s ##<< deprecated
, RadVar.s ##<< deprecated
, Suffix.s ##<< deprecated
## identifier for NEE column defaults
## (see \code{\link{sEddyProc_sMDSGapFillAfterUstar}}).
, ... ##<< further arguments to
## \code{\link{partGLExtractStandardData}}, such as \code{PotRadVar}
, 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)
, lrcFitter = RectangularLRCFitter() ##<< R5 class instance
## responsible for fitting the light response curve.
## Current possibilities are \code{RectangularLRCFitter()},
## \code{NonrectangularLRCFitter()},
## and \code{LogisticSigmoidLRCFitter()}.
) {
varNamesDepr <- c(
"NEEVar.s","TempVar.s","VPDVar.s","RadVar.s","Suffix.s")
varNamesNew <- c(
"NEEVar","TempVar","VPDVar","RadVar","suffix")
iDepr = which(!c(
missing(NEEVar.s),missing(TempVar.s),missing(VPDVar.s)
,missing(RadVar.s),missing(Suffix.s)))
if (length(iDepr)) warning(
"Argument names ",varNamesDepr[iDepr]," have been deprecated."
," Please, use instead ", varNamesNew[iDepr])
##details<<
## Daytime-based partitioning of measured net ecosystem fluxes into
## gross primary production (GPP) and ecosystem respiration (Reco)
##
## The fit to the light-response-curve is done by default using the Rectangular
## hyperbolic function, as in Lasslop et al. (2010)
## Alternative fittings can be used by providing the corresponding subclass of
## \code{\link{LightResponseCurveFitter-class}} to \code{lrcFitter} argument.
## (see \code{\link{LightResponseCurveFitter_predictGPP}})
#
##author<< 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
suffixDash <- paste( (if (fCheckValString(suffix)) "_" else "")
, suffix, sep = "") # used to compute default NEEVar
if (isVerbose) message('Start daytime flux partitioning for variable '
, NEEVar, ' with temperature ', TempVar, '.')
dsR <- partGLExtractStandardData(ds, NEEVar = NEEVar, TempVar = TempVar
, VPDVar = VPDVar, RadVar = RadVar , suffix = suffix
, ..., controlGLPart = controlGLPart)
##value<<
## \item{Reco_DT_<suffix>}{predicted ecosystem respiration: mumol CO2/m2/s}
## \item{GPP_DT_<suffix>}{predicted gross primary production mumol CO2/m2/s}
## \item{<LRC>}{Further light response curve (LRC) parameters and
## their standard deviation depend on the used LRC
## (e.g. for the non-rectangular LRC
## see \code{\link{NonrectangularLRCFitter_getParameterNames}}).
## They are estimated for windows and are reported with the first record
## of the window}
dsAns0 <- tibble(
FP_VARnight = rep(NA_real_, nrow(ds)) ##<< NEE filtered
## for nighttime records (others NA)
, FP_VARday = NA_real_ ##<< NEE filtered for daytime
## records (others NA)
, NEW_FP_Temp = dsR$Temp ##<< temperature after filtering for quality
## flag degree Celsius
, NEW_FP_VPD = dsR$VPD ##<< vapour pressure deficit after filtering
## for quality flag, hPa
, FP_RRef_Night = NA_real_ ##<< basal respiration estimated
## from nighttime (W / m2)
, 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
, FP_errorcode = NA_integer_ ##<< information why LRC-fit was not
## successful or was rejected, see result of
## \code{\link{LightResponseCurveFitter_fitLRC}}
, FP_GPP2000 = NA_real_ ##<< predicted GPP at VPD = 0 and PAR = 2000:
## a surrogate for maximum photosynthetic capacity
, FP_OPT_VPD = vector(mode = "list", length = nrow(ds)) ##<< list object
## of fitting results including iOpt and covParms
, FP_OPT_NoVPD = vector(mode = "list", length = nrow(ds)) ##<< same as
## FP_OPT_VPD holding optimization results with fit
## neglecting the VPD effect
)
#
# append LRC parameter result columns
lrcParNames <- lrcFitter$getParameterNames()
lrcParNames <- c(lrcParNames, paste0(lrcParNames, "_sd"))
FP_lrcParNames <- paste0("FP_", lrcParNames)
tmp <- matrix(NA_real_, nrow = nrow(ds), ncol = length(FP_lrcParNames)
, dimnames = list(NULL, FP_lrcParNames) )
dsAns <- cbind(dsAns0, tmp)
dsAns$FP_VARnight <- ifelse(dsR$isNight, dsR$NEE, NA)
attr(dsAns$FP_VARnight, 'varnames') <- paste(attr(dsR$NEE, 'varnames')
, '_night', sep = '')
attr(dsAns$FP_VARnight, 'units') <- attr(dsR$NEE, 'units')
dsAns$FP_VARday <- ifelse(dsR$isDay, dsR$NEE, NA)
attr(dsAns$FP_VARday, 'varnames') <- paste(attr(dsR$NEE, 'varnames')
, '_day', sep = '')
attr(dsAns$FP_VARday, 'units') <- attr(dsR$NEE, 'units')
#
#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")
#
##seealso<< \code{partGLFitNightTimeTRespSens}
isUsingFixedTempSens <- length(controlGLPart$fixedTempSens) &&
all(is.finite(controlGLPart$fixedTempSens$E0)) &&
all(is.finite(controlGLPart$fixedTempSens$sdE0))
dsTempSens <- if (!isUsingFixedTempSens ) {
dsTempSens <- partGLFitNightTimeTRespSens(dsR
, nRecInDay = nRecInDay
, controlGLPart = controlGLPart
, isVerbose = isVerbose
)
} else {
dsTempSens <- controlGLPart$fixedTempSens
nWindow <- length(getStartRecsOfWindows(nrow(dsR), nRecInDay = nRecInDay))
# if one row is given, use it for all windows
if (nrow(dsTempSens) == 1L) dsTempSens <- do.call(rbind,
lapply(1:nWindow, function(i) dsTempSens))
if (nrow(dsTempSens) != nWindow) stop(
"when providing controlGLPart$fixedTempSens then it must be 1 row"
,"or rows must macht the number of windows (", nWindow, "), but was "
, nrow(dsTempSens))
# if no RRef is given, estimate from nighttime data
if (!length(dsTempSens$RRef) || all(is.na(dsTempSens$RRef)) ) {
resRef15 <- simplifyApplyWindows(applyWindows(dsR
, partGLFitNightRespRefOneWindow
, nRecInDay = nRecInDay
, E0Win = dsTempSens
, controlGLPart = controlGLPart
, isVerbose = isVerbose
))
dsTempSens$RRef <- resRef15$RRef
}
# fill missing RRef by previous window
dsTempSens$RRef <- fillNAForward(dsTempSens$RRef, firstValue =
# on NA at the start of the series, take the first occuring
# finite value
dsTempSens$RRef[which(is.finite(dsTempSens$RRef))[1] ] )
dsTempSens
}
##seealso<< \code{partGLFitLRCWindows}
resParms <- resParmsWithVPD <- partGLFitLRCWindows(dsR
, nRecInDay = nRecInDay
, dsTempSens = dsTempSens
, controlGLPart = controlGLPart
, lrcFitter = lrcFitter
, isVerbose = isVerbose
)
# if no windows was fitted, return error. Else error on missing columns
# (no parameter columns returned)
if (sum(resParms$convergence == 0) == 0) stop(
"could not fit a single window in dayTime partitioning. Check the data.")
# 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$iCentralRec, c("FP_RRef_Night", "FP_qc", "FP_errorcode"
, FP_lrcParNames, "FP_GPP2000", "FP_OPT_VPD")] <-
resParms[, c("RRef_night", "parms_out_range", "convergence"
, lrcParNames, "GPP2000", "resOpt")]
##seealso<< \code{partGLInterpolateFluxes}
dsAnsFluxes <- partGLInterpolateFluxes(
#dsR$Rg
ds[[RadVar]]
#, dsAns$NEW_FP_VPD, dsAns$NEW_FP_Temp
, ds[[VPDVar]]
, ds[[TempVar]] # prediction using non-filtered, i.e. gap-filled
, resParms
, controlGLPart = controlGLPart
, lrcFitter = lrcFitter
, isVerbose = isVerbose
)
if (!controlGLPart$isNeglectVPDEffect &&
controlGLPart$isRefitMissingVPDWithNeglectVPDEffect) {
##details<<
## While the extrapolation uses filled data, the parameter optimization
## may use only measured data, i.e. with specified quality flag.
## Even with using filled VPD, there may be large gaps that have not been
## filled.
## With the common case where VPD is missing for fitting the LRC, by default
## (with \code{controlGLPart$isRefitMissingVPDWithNeglectVPDEffect = TRUE})
## is to redo the estimation of LRC parameters with neglecting the VPD-effect.
## Next, in the predictions (rows) with missing VPD are then replaced
## with predictions
## based on LRC-fits that neglected the VPD effect.
iNAVPD <- which(is.na(ds[[VPDVar]] & is.na(dsAnsFluxes$GPP)))
if (length(iNAVPD)) {
message(" could not predict GPP in ", length(iNAVPD)
, " cases due to missing VPD.")
message(" Therefore refitting LightResponseCurve with option "
,"isNeglectVPDEffect = TRUE")
ctrlNeglectVPD <- within(controlGLPart, isNeglectVPDEffect <- TRUE)
resParms <- resParmsWithoutVPD <- partGLFitLRCWindows(dsR
, nRecInDay = nRecInDay
, dsTempSens = dsTempSens
, controlGLPart = ctrlNeglectVPD
, lrcFitter = lrcFitter
)
dsAns$FP_OPT_NoVPD[resParmsWithoutVPD$iCentralRec] <-
resParmsWithoutVPD$resOpt
dsAnsFluxes2 <- partGLInterpolateFluxes(dsR$Rg
#, dsAns$NEW_FP_VPD, dsAns$NEW_FP_Temp
, ds[[VPDVar]]
, ds[[TempVar]] # do prediction also using gap-Filled values
, resParmsWithoutVPD
, controlGLPart = ctrlNeglectVPD
, lrcFitter = lrcFitter
)
dsAnsFluxes[iNAVPD, ] <- dsAnsFluxes2[iNAVPD, ]
}
}
if (isTRUE(controlGLPart$useNightimeBasalRespiration))
dsAnsFluxes <- .computeRecoNight(dsAnsFluxes, dsR, dsTempSens)
nNAGPP <- sum(is.na(dsAnsFluxes$GPP))
if (nNAGPP) warning("could not predict GPP in ", nNAGPP, " cases.")
#
dsAns$FP_dRecPar <- dsAnsFluxes$dRecNextEstimate
# copy quality flag from parameter row
# first copy from iCentralRec to colNameAssoc
# next transfer by dRecNextEstimate to each row
iFiniteMeanRec <- which(is.finite(resParms[[colNameAssoc]]))
dsAns$FP_qc[resParms[[colNameAssoc]][iFiniteMeanRec] ] <-
dsAns$FP_qc[resParms$iCentralRec[iFiniteMeanRec] ]
dsAns$FP_qc <- dsAns$FP_qc[1:nrow(dsAns) + dsAns$FP_dRecPar]
#set quality flag to 2 where next parameter estimate is more than 14 days away
dsAns$FP_qc[abs(dsAns$FP_dRecPar) > (14 * nRecInDay)] <- 2L
#dsAns[is.finite(dsAns$FP_beta), ]
#
RecoDTVar.s <- paste0('Reco_DT', suffixDash)
GPPDTVar.s <- paste0('GPP_DT', suffixDash)
RecoDTSdVar.s <- paste0(RecoDTVar.s, "_SD")
GPPDTSdVar.s <- paste0(GPPDTVar.s, "_SD")
dsAns[[RecoDTVar.s]] <- dsAnsFluxes$Reco
attr(dsAns[[RecoDTVar.s]], 'varnames') <- RecoDTVar.s
attr(dsAns[[RecoDTVar.s]], 'units') <- attr(dsR$NEE, 'units')
dsAns[[GPPDTVar.s]] <- dsAnsFluxes$GPP
attr(dsAns[[GPPDTVar.s]], 'varnames') <- GPPDTVar.s
attr(dsAns[[GPPDTVar.s]], 'units') <- attr(dsR$NEE, 'units')
if (controlGLPart$isSdPredComputed) {
dsAns[[RecoDTSdVar.s]] <- dsAnsFluxes$sdReco
attr(dsAns[[RecoDTSdVar.s]], 'varnames') <- RecoDTSdVar.s
attr(dsAns[[RecoDTSdVar.s]], 'units') <- attr(dsR$NEE, 'units')
dsAns[[GPPDTSdVar.s]] <- dsAnsFluxes$sdGPP
attr(dsAns[[GPPDTSdVar.s]], 'varnames') <- GPPDTSdVar.s
attr(dsAns[[GPPDTSdVar.s]], 'units') <- attr(dsR$NEE, '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)
}
#' @export
partGLControl <- function(
### Default list of parameters for Lasslop 2010 daytime flux partitioning
LRCFitConvergenceTolerance = 1e-3 ##<< convergence criterion for rectangular
## light response curve 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.
, nLRCFitConvergenceTolerance = 1e-3 ##<< convergence criterion for
## nonrectangular light response curve fit.
## Here its a factor of machine tolerance.
, 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 = TRUE ##<< 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
, isUsingLasslopQualityConstraints = FALSE ##<< set to TRUE to avoid
## quality constraints additional to Lasslop 2010
, 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 (radiation, temperature, 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.
, fixedTRefAtNightTime = NA ##<< if a finite value (degree Centigrade)
## is given, it is used instead of median data temperature as reference
## temperature in estimation of temperature sensitivity from night data
, isExtendTRefWindow = TRUE ##<< set to FALSE to avoid successively
## extending the night-time window
## in order to estimate a temperature sensitivity where previous estimates
## failed
, 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
, isNeglectPotRadForNight = FALSE ##<< set to TRUE to not use potential
## radiation in determining night-time data.
, NRHRfunction = FALSE ##<< deprecated: Flag if TRUE use the NRHRF
## for partitioning; Now use \code{lrcFitter = NonrectangularLRCFitter()}
, isNeglectVPDEffect = FALSE ##<< set to TRUE to avoid using VPD in the
## computations. This may help when VPD is rarely measured.
, isRefitMissingVPDWithNeglectVPDEffect = TRUE ##<< set to FALSE to avoid
## repeating estimation
## with \code{isNeglectVPDEffect = TRUE} trying to predict when VPD
## is missing
, fixedTempSens = data.frame( ##<< data.frame
## of one row or nRow = nWindow
## corresponding to return value of \code{partGLFitNightTimeTRespSens}
## While column \code{RRef} is used only as a prior and initial value for
## the daytime-fitting and can be NA,
## \code{E0} is used as given temperature sensitivity and varied according
## to \code{sdE0} in the bootstrap.
E0 = NA_real_, sdE0 = NA_real_, RRef = NA_real_)
, replaceMissingSdNEEParms = c(perc = 0.2, minSd = 0.7) ##<< parameters for
## replacing missing standard deviation of NEE.
## see \code{replaceMissingSdByPercentage}.
## Default sets missing uncertainty to 20% of NEE but at least 0.7
## flux-units (usually mumol CO2 / m2 / s).
## Specify c(NA, NA) to avoid replacing missings in standard deviation of
## NEE and to omit those records from LRC fit.
, neglectNEEUncertaintyOnMissing = FALSE ##<< If set to TRUE: if there are
## records with missing uncertainty of NEE inside one window,
## set all uncertainties to 1.
## This overrules option replaceMissingSdNEEParms.
, minPropSaturation = NA ##<< quality criterion for sufficient data
## in window. If GPP prediction of highest PAR of window is less than
## minPropSaturation * (GPP at light-saturation, i.e. beta)
## this indicates that PAR is not sufficiently high to constrain the
## shape of the LRC
, useNightimeBasalRespiration = FALSE ##<< set to TRUE to estimate
## nighttime respiration based on basal respiration estimated on
## nighttime data instead of basal respiration estimated from daytime
## data. This implements the modified daytime method from
## Keenan 2019 (doi:10.1038/s41559-019-0809-2)
) {
##author<< TW
##seealso<< \code{\link{partitionNEEGL}}
##description<<
## For highest compatibility to the pvWave code of G.Lasslop
## (used by first BGC-online tool)
## see function \code{\link{partGLControlLasslopCompatible}}.
if (NRHRfunction) stop("option 'NRHRfunction' is deprecated."
," Use instead in partitionNEEGL argument:"
,"lrcFitter = NonrectangularLRCFitter()")
if (isTRUE(neglectNEEUncertaintyOnMissing) )
replaceMissingSdNEEParms <- c(NA, NA)
ctrl <- list(
LRCFitConvergenceTolerance = LRCFitConvergenceTolerance
, nLRCFitConvergenceTolerance = nLRCFitConvergenceTolerance
, nBootUncertainty = nBootUncertainty
, minNRecInDayWindow = minNRecInDayWindow
, isAssociateParmsToMeanOfValids = isAssociateParmsToMeanOfValids
, isLasslopPriorsApplied = isLasslopPriorsApplied
, isUsingLasslopQualityConstraints = isUsingLasslopQualityConstraints
, isSdPredComputed = isSdPredComputed
, isFilterMeteoQualityFlag = isFilterMeteoQualityFlag
, isBoundLowerNEEUncertainty = isBoundLowerNEEUncertainty
, fixedTRefAtNightTime = fixedTRefAtNightTime
, isExtendTRefWindow = isExtendTRefWindow
, smoothTempSensEstimateAcrossTime = smoothTempSensEstimateAcrossTime
, isNeglectPotRadForNight = isNeglectPotRadForNight
, isNeglectVPDEffect = isNeglectVPDEffect
, isRefitMissingVPDWithNeglectVPDEffect = isRefitMissingVPDWithNeglectVPDEffect
, fixedTempSens = fixedTempSens
, replaceMissingSdNEEParms = replaceMissingSdNEEParms
, neglectNEEUncertaintyOnMissing = neglectNEEUncertaintyOnMissing
, minPropSaturation = minPropSaturation
, useNightimeBasalRespiration = useNightimeBasalRespiration
)
##value<< list with entries of given arguments.
ctrl
}
attr(partGLControl, "ex") <- function() {
partGLControl(nBootUncertainty = 40L)
}
#' @export
partGLControlLasslopCompatible <- function(
### Daytime flux partitioning parms compatible with with the pvWave
nBootUncertainty = 0L ##<< 0: Derive uncertainty from
## curvature of a single fit, neglecting the uncertainty of previously
## estimated temperature sensitivity, E0
, minNRecInDayWindow = 10L ##<< Minimum number of 10 valid records
## for regression in a single window
, isAssociateParmsToMeanOfValids = FALSE ##<< associate parameters to
## the first record of the window for interpolation instead of mean across
## valid records inside a window
, isLasslopPriorsApplied = TRUE ##<< Apply fixed Lasslop priors
## in LRC fitting.
, isUsingLasslopQualityConstraints = TRUE ##<< avoid quality constraints
## additional to the ones in Lasslop 2010
, isBoundLowerNEEUncertainty = FALSE ##<< FALSE: avoid adjustment of very
## low uncertainties before
## day-Time fitting that avoids the high leverage those records with
## unreasonable low uncertainty.
, fixedTRefAtNightTime = 15 ##<< use fixed (degree Centigrade)
## temperature sensitivity
## instead of median data temperature as reference temperature in
## estimation of temperature sensitivity from night data
, isExtendTRefWindow = FALSE ##<< avoid successively extending the
##night-time window
## in order to estimate a temperature sensitivity where previous
## estimates failed
, smoothTempSensEstimateAcrossTime = FALSE ##<< FALSE: use independent
## estimates of temperature
## sensitivity on each windows instead of a vector of E0 that is
## smoothed over time
, isRefitMissingVPDWithNeglectVPDEffect = FALSE ##<< FALSE: avoid
## repeating estimation with \code{isNeglectVPDEffect = TRUE}
, minPropSaturation = NA ##<< NA: avoid quality constraint of
## sufficient saturation in data
## This option is overruled, i.e. not considered, if option
## isUsingLasslopQualityConstraints = TRUE.
, isNeglectVPDEffect = FALSE ##<< FALSE: do not neglect VPD effect
, replaceMissingSdNEEParms = c(NA, NA) ##<< do not replace missing NEE,
## but see option
, neglectNEEUncertaintyOnMissing = TRUE ##<< if there are records with
## missing uncertainty of NEE inside one window, set all sdNEE to 1.
## This overrules option replaceMissingSdNEEParms.
, ... ##<< further arguments to \code{\link{partGLControl}}
) {
##author<< TW
##seealso<< \code{\link{partGLControl}}
partGLControl(
nBootUncertainty = nBootUncertainty
, minNRecInDayWindow = minNRecInDayWindow
, isAssociateParmsToMeanOfValids = isAssociateParmsToMeanOfValids
, isLasslopPriorsApplied = isLasslopPriorsApplied
, isUsingLasslopQualityConstraints = isUsingLasslopQualityConstraints
, isBoundLowerNEEUncertainty = isBoundLowerNEEUncertainty
, fixedTRefAtNightTime = fixedTRefAtNightTime
, isExtendTRefWindow = isExtendTRefWindow
, smoothTempSensEstimateAcrossTime = smoothTempSensEstimateAcrossTime
, isNeglectVPDEffect = isNeglectVPDEffect
, isRefitMissingVPDWithNeglectVPDEffect = isRefitMissingVPDWithNeglectVPDEffect
, replaceMissingSdNEEParms = replaceMissingSdNEEParms
, neglectNEEUncertaintyOnMissing = neglectNEEUncertaintyOnMissing
, minPropSaturation = minPropSaturation
, ...
)
}
attr(partGLControlLasslopCompatible, "ex") <- function() {
partGLControlLasslopCompatible()
}
#' @export
partGLExtractStandardData <- function(
### Relevant columns from original input with defined names
ds ##<< dataset with all the specified input columns and
## full days in equidistant times
, NEEVar = paste0('NEE', suffixDash, '_f') ##<< Variable of NEE
, QFNEEVar = if (!missing(QFNEEVar.s)) QFNEEVar.s else paste0('NEE', suffixDash, '_fqc') ##<< Quality
## flag of variable
, QFNEEValue = if (!missing(QFNEEValue.n)) QFNEEValue.n else 0 ##<< Value of quality flag for
## _good_ (original) data
, NEESdVar = if (!missing(NEESdVar.s)) NEESdVar.s else paste0('NEE', suffixDash, '_fsd') ##<< Variable of
## standard deviation of net ecosystem fluxes
, TempVar = paste0('Tair_f') ##<< Filled air or soil
## temperature variable (degC)
, QFTempVar = if (!missing(QFTempVar.s)) QFTempVar.s else paste0('Tair_fqc') ##<< Quality flag of
## filled temperature variable
, QFTempValue = if (!missing(QFTempValue.n)) QFTempValue.n else 0 ##<< Value of temperature quality flag
##for _good_ (original) data
, VPDVar = if (!missing(VPDVar.s)) VPDVar.s else paste0('VPD_f') ##<< Filled Vapor Pressure Deficit, VPD (hPa)
, QFVPDVar = if (!missing(QFVPDVar.s)) QFVPDVar.s else paste0('VPD_fqc') ##<< Quality flag of filled VPD variable
, QFVPDValue = if (!missing(QFVPDValue.n)) QFVPDValue.n else 0 ##<< Value of VPD quality flag for
## _good_ (original) data
, RadVar = if (!missing(RadVar.s)) RadVar.s else 'Rg_f' ##<< Filled radiation variable
, QFRadVar = if (!missing(QFRadVar.s)) QFRadVar.s else paste0('Rg_fqc') ##<< Quality flag of filled radiation variable
, QFRadValue = if (!missing(QFRadValue.n)) QFRadValue.n else 0 ##<< Value of radiation quality flag for
## _good_ (original) data
, PotRadVar = if (!missing(PotRadVar.s)) PotRadVar.s else "PotRad_NEW" ##<< Variable name of potential rad. (W / m2)
, suffix = if (!missing(Suffix.s)) Suffix.s else "" ##<< string inserted into column names before
## identifier for NEE column defaults
## (see \code{\link{sEddyProc_sMDSGapFillAfterUstar}}).
, NEEVar.s ##<< deprecated
, QFNEEVar.s ##<< deprecated
, QFNEEValue.n ##<< deprecated
, NEESdVar.s ##<< deprecated
, TempVar.s ##<< deprecated
, QFTempVar.s ##<< deprecated
, QFTempValue.n ##<< deprecated
, VPDVar.s ##<< deprecated
, QFVPDVar.s ##<< deprecated
, QFVPDValue.n ##<< deprecated
, RadVar.s ##<< deprecated
, QFRadVar.s ##<< deprecated
, QFRadValue.n ##<< deprecated
, PotRadVar.s ##<< deprecated
, Suffix.s ##<< deprecated
, controlGLPart = partGLControl() ##<< further default parameters,
## see \code{\link{partGLControl}}
) {
if (!missing(NEEVar.s)) NEEVar <- NEEVar.s # in default, lines too wide
if (!missing(TempVar.s)) TempVar <- TempVar.s
varNamesDepr <- c(
"NEEVar.s","TempVar.s","VPDVar.s","RadVar.s","Suffix.s"
,"QFNEEVar.s" ,"QFNEEValue.n" ,"QFTempVar.s" ,"QFTempValue.n"
,"QFVPDVar.s" ,"QFVPDValue.n" ,"QFRadVar.s" ,"QFRadValue.n" ,"PotRadVar.s"
)
varNamesNew <- c(
"NEEVar","TempVar","VPDVar","RadVar","suffix"
,"QFNEEVar" ,"QFNEEValue" ,"QFTempVar" ,"QFTempValue"
,"QFVPDVar" ,"QFVPDValue" ,"QFRadVar" ,"QFRadValue" ,"PotRadVar"
)
iDepr = which(!c(
missing(NEEVar.s),missing(TempVar.s),missing(VPDVar.s),missing(RadVar.s)
,missing(Suffix.s),missing(QFNEEVar.s),missing(QFNEEValue.n)
,missing(QFTempVar.s),missing(QFTempValue.n),missing(QFVPDVar.s)
,missing(QFVPDValue.n),missing(QFRadVar.s),missing(QFRadValue.n)
,missing(PotRadVar.s)))
if (length(iDepr)) warning(
"Argument names ",varNamesDepr[iDepr]," have been deprecated."
," Please, use instead ", varNamesNew[iDepr])
# Check if specified columns exist in sDATA or sTEMP and if numeric
# and plausible. Then apply quality flag
suffixDash <- paste( (if (fCheckValString(suffix)) "_" else "")
, suffix, sep = "")
# Check if specified columns exist in sDATA or sTEMP and if numeric
# and plausible. Then apply quality flag
fCheckColNames(ds, c(NEEVar, QFNEEVar, TempVar, QFTempVar, RadVar, VPDVar
, QFRadVar, PotRadVar, NEESdVar), 'sGLFluxPartition')
fCheckColNum(ds, c(NEEVar, QFNEEVar, TempVar, QFTempVar, RadVar, VPDVar
, QFRadVar, PotRadVar, NEESdVar), 'sGLFluxPartition')
fCheckColPlausibility(ds, c(NEEVar, QFNEEVar, TempVar, QFTempVar, VPDVar
, RadVar, QFRadVar, PotRadVar), 'sGLFluxPartition')
NEEFiltered <- fSetQF(ds, NEEVar, QFNEEVar, QFNEEValue, 'sGLFluxPartition')
#
# Apply quality flag for temperature and VPD
# TODO: docu meteo filter, standard FALSE
NEW_FP_Temp <- if (isTRUE(controlGLPart$isFilterMeteoQualityFlag) ) fSetQF(
ds, TempVar, QFTempVar, QFTempValue, 'partGLExtractStandardData') else
ds[[TempVar]]
NEW_FP_VPD <- if (isTRUE(controlGLPart$isFilterMeteoQualityFlag) )
fSetQF(ds, VPDVar, QFVPDVar, QFVPDValue, 'partGLExtractStandardData') else
ds[[VPDVar]]
NEW_FP_Rg <- if (isTRUE(controlGLPart$isFilterMeteoQualityFlag) )
fSetQF(ds, RadVar, QFRadVar, QFRadValue, 'partGLExtractStandardData') else
ds[[RadVar]]
#
# 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
#! New code: PotRad in sGLFluxPartition: Slightly different subset than
#PV-Wave due to time zone correction (avoids timezone offset between Rg and PotRad)
isPotRadZero <- if (isTRUE(controlGLPart$isNeglectPotRadForNight) )
TRUE else (ds[[PotRadVar]] <= 0)
isNotPotRadZero <- if (isTRUE(controlGLPart$isNeglectPotRadForNight) )
TRUE else (ds[[PotRadVar]] > 0)
isNight <- (ds[[RadVar]] <= 4 & isPotRadZero)
# 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, twutz does not understand this comment
isDay = (ds[[RadVar]] > 4 & isNotPotRadZero)
#
##value<< a data.frame with columns
dsR <- data.frame(
sDateTime = ds[[1]] ##<< first column of \code{ds},
## usually the time stamp
## not used, but usually first column is a DateTime is kept
## for aiding debug
, NEE = NEEFiltered ##<< NEE filtered for quality flay
, sdNEE = ds[[NEESdVar]] ##<< standard deviation of NEE
## with missing values replaced
, Temp = NEW_FP_Temp ##<< Temperature, quality filtered
## if isTRUE(controlGLPart$isFilterMeteoQualityFlag)
, VPD = NEW_FP_VPD ##<< Water pressure deficit, quality
## filtered if isTRUE(controlGLPart$isFilterMeteoQualityFlag)
, Rg = NEW_FP_Rg ##<< Incoming radiation
, isDay = isDay ##<< Flag that is true for daytime records
, isNight = isNight ##<< Flag that is true for nighttime records
)
##details<<
## The LRC fit usually weights NEE records by its uncertainty. In order
## to also use
## records with missing \code{NEESdVar}, uncertainty of the missing values
## is by default set
## to a conservatively high value, parameterized by
## \code{controlGLPart$replaceMissingSdNEEParms)}.
## Controlled by argument \code{replaceMissingSdNEEParms} in
## \code{\link{partGLControl}}, but overruled
## by argument \code{neglectNEEUncertaintyOnMissing}.
if (!controlGLPart$neglectNEEUncertaintyOnMissing)
dsR$sdNEE <- replaceMissingSdByPercentage(dsR$sdNEE, dsR$NEE
, controlGLPart$replaceMissingSdNEEParms[1]
, controlGLPart$replaceMissingSdNEEParms[2] )
dsR
}
partGLFitLRCWindows <- function(
### Estimate successive Rectangular Hyperbolic Light Response Curve parameters
ds ##<< data.frame with numeric columns NEE, sdNEE,
## Temp (degC), VPD, Rg, and logical columns isNight and isDay
, winSizeDays = 4L ##<< Window size in days for daytime fits
, strideInDays = 2L ##<< step in days for shifting the windows
, nRecInDay = 48L ##<< number of records within one day
## (for half-hourly data its 48)
, dsTempSens ##<< data.frame that reports for each window
## temperature sensitivity parameters E0 and RRef
, isVerbose = TRUE ##<< set to FALSE to suppress messages
, controlGLPart = partGLControl() ##<< list of further default parameters
, lrcFitter ##<< R5 class instance responsible for fitting the
## light response curve
) {
##seealso<< \code{\link{partGLFitLRCOneWindow}}
if (isVerbose) message(
" Estimating light response curve parameters from day time NEE "
, appendLF = FALSE)
resLRC <- applyWindows(ds, partGLFitLRCOneWindow
, prevRes = list(resOpt = list(thetaOpt = rep(NA_real_, 6L)))
, winSizeInDays = 4L
, isVerbose = isVerbose
, nRecInDay = nRecInDay
#
, E0Win = dsTempSens
, controlGLPart = controlGLPart
, lrcFitter = lrcFitter
)
lrcSummary <- lapply(resLRC$resFUN, "[[", "summary")
iNoSummary <- which(sapply(lrcSummary, length) == 0)
if (length(iNoSummary) ) {
stop("expected summary returned by all fits, but found missing summaries.")
# # put dummy NA data.frame where no fit was obtained
# dummySummary <- lrcSummary[-iNoSummary][[1]]
# dummySummary[] <- NA
# lrcSummary[iNoSummary] <- list(dummySummary)
}
resOptList <- lapply(resLRC$resFUN, "[[", "resOpt")
resOptDf <- tibble(resOpt = resOptList)
resParms <- as_tibble(cbind(resLRC$winInfo, bind_rows(lrcSummary), resOptDf))
#table(resParms$convergence)
#E0_night equals E0, but uncertainty might differ
resParms$E0_bootstrap_sd <- resParms$E0_sd # due to bootstrap,
# this may differ, save before overriding by night-time estimate
resParms$E0 <- dsTempSens$E0 # report E0 and E0_sd from nighttime
# even when LRC fit did not converge
resParms$E0_sd <- dsTempSens$sdE0
resParms$RRef_night <- dsTempSens$RRef
# summary$iMeanRec yet based on window instead of entire time, need to
# add beginning of window
resParms$iMeanRec <- resParms$iRecStart - 1L + resParms$iMeanRec
# # omit records where NULL was returned
# iWinNoFit <- which(is.na(resParms$parms_out_range) )
# if (length(iWinNoFit) ) {
# resParms <- resParms[-iWinNoFit, ]
# resParms$resOptList <- resParms$resOptList[-iWinNoFit]
# }
##value<< a tibble with a row for each window
## providing information on window size, parameters estimated
## from light response curve (LRC),
## and their standard deviation.
## Moreover, and estimated from night-time data after smoothing
## and forward-filling:
## the uncertainty of temperature sensitivity \code{E0_night_sd},
## And the respiration at reference temperature \code{RRef_night}.
## Finally it contains a list column\code{resOpt} with the optimization results.
resParms
}
partGLFitLRCOneWindow <- function(
### Estimate Rectangular Hyperbolic Light Response Curve parameters
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
, lrcFitter = RectangularLRCFitter() ##<< R5 class instance responsible
## for fitting the light response curve
) {
##author<< TW
##seealso<< \code{\link{partitionNEEGL}}
##details<<
## Estimation as in Lasslop et al., 2010 for successive periods, i.e. windows.
## Fitting is done on a subset of the data where isDay, NEE, sdNEE, Temp,
## Rg, and VPD are all non-NA
## and isDay is TRUE
isValidDayRecNoSdNEEConstraint <- !is.na(ds$isDay) & ds$isDay &
!is.na(ds$NEE) & !is.na(ds$Temp) & !is.na(ds$Rg)
##details<<
## With option \code{controlGLPart$neglectNEEUncertaintyOnMissing = TRUE}
## all sdNEE are set to 1 if there is any missing sdNEE at
## otherwise valid records
if (isTRUE(controlGLPart$neglectNEEUncertaintyOnMissing) &&
any(!is.finite(ds$sdNEE[isValidDayRecNoSdNEEConstraint])) ) {
ds$sdNEE <- 1L
}
isValidDayRecNoVPDConstraint <- isValidDayRecNoSdNEEConstraint &
!is.na(ds$sdNEE)
##details<<
## If \code{controlGLPart$isNeglectVPDEffect = TRUE}, also records
## with VPD = NA maybe valid
isValidDayRec <- if (isTRUE(controlGLPart$isNeglectVPDEffect) )
isValidDayRecNoVPDConstraint else
isValidDayRecNoVPDConstraint & !is.na(ds$VPD)
iMeanRecInDayWindow <- as.integer(round(mean(which(isValidDayRec))))
#
# if too few records or
# if no temperature-respiration relationship could be found,
# indicate no-fit, but report Window properties
E0 <- E0Win$E0[winInfo$iWindow]
# when provided fixed E0 the column sdEO does not exist
sdE0 <- if (length(E0Win$sdE0)) E0Win$sdE0[winInfo$iWindow] else NA
RRefNight <- E0Win$RRef[winInfo$iWindow]
getNAResult <- function(convergenceCode) {
list(
resOpt = NULL
, summary = data.frame(
nValidRec = sum(isValidDayRec)
, iMeanRec = iMeanRecInDayWindow
, convergence = convergenceCode
, E0 = E0
, E0_sd = sdE0
, RRefNight = RRefNight
)
, isValid = FALSE
)}
if (is.na(E0) ) return(getNAResult(1010L))
#
##details<< If there are too few records (<
## \code{controlGLPart$minNRecInDayWindow}), then
## the constraint on non-NA VPD is neglected and
## \code{controlGLPart$isNeglectVPDEffect} is set to TRUE
## If there are still too few records, an NA-result with convergence
## code 1011L is returned.
if ( (sum(isValidDayRec) < controlGLPart$minNRecInDayWindow) ) {
controlGLPart$isNeglectVPDEffect <- TRUE
isValidDayRec <- isValidDayRecNoVPDConstraint
iMeanRecInDayWindow <- as.integer(round(mean(which(isValidDayRec))))
if ( (sum(isValidDayRec) < controlGLPart$minNRecInDayWindow) )
return(getNAResult(1011L))
}
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.
#better report NA and care for it properly: if (is.na(iMeanRecInDayWindow))
# iMeanRecInDayWindow <- as.integer(nrow(ds)%/%2)
#TODO firstRecInDayWindow.i <- which(SubsetDayPeriod.b)[1]
# the rownumber of the first record inside the day window
##seealso<< \code{\link{partGLEstimateTempSensInBoundsE0Only}}
#
##seealso<< \code{\link{LightResponseCurveFitter_fitLRC}}
#lrcFitter <- RectangularLRCFitter()
resOpt <- resOpt0 <- lrcFitter$fitLRC(dsDay, E0 = E0, sdE0 = sdE0
, RRefNight = RRefNight
, controlGLPart = controlGLPart
, lastGoodParameters = prevRes$resOpt$thetaOpt
)
if (!is.finite(resOpt$thetaOpt[1]) ) {
return(getNAResult(resOpt$convergence))
}
sdTheta <- resOpt$thetaOpt; sdTheta[] <- NA
sdTheta[resOpt$iOpt] <- sqrt(diag(resOpt$covParms)[resOpt$iOpt])
#
# record valid fits results
#as.data.frame(t(resOpt$thetaOpt))
#if (as.POSIXlt(dsDay$sDateTime[1])$mday + 2L >= 11) recover()
#if (as.POSIXlt(dsDay$sDateTime[1])$mday + 2L >= 27) recover()
# save(dsDay, file = "tmp / dsDayDebug.RData")
#
# compute GPP2000 without VPD effect
GPP2000 <- lrcFitter$predictLRC(
resOpt$thetaOpt, 2000, VPD = NA, Temp = NA, fixVPD = TRUE)$GPP
ans <- list(
resOpt = resOpt
, summary = cbind(data.frame(
nValidRec = nrow(dsDay)
, iMeanRec = iMeanRecInDayWindow
, convergence = resOpt$convergence
, parms_out_range = as.integer(!all(1:5 %in% resOpt$iOpt) )
)
, as.data.frame(t(resOpt$thetaOpt))
, as.data.frame(t(structure(sdTheta, names = paste0(names(sdTheta), "_sd"))))
, GPP2000 = GPP2000
, isValid = TRUE
)
)
return(ans)
}
.bootStrapLRCFit <- function(
### Compute parameters uncertainty by bootstrap
theta0, iOpt, dsDay, sdE_0.n, parameterPrior, controlGLPart
, lrcFitter ##<< Light Response Curve R5 instance
, iPosE0 = 5L ##<< position (integer scalar) of temperature
##sensitivity in parameter vector
#should be dealt with iOpt: isNeglectVPD = FALSE
) {
##value<<
## Matrix with each row a parameter estimate on a different bootstrap sample.
## Some of the rows might be NA, if the subsampled dataset could not be fitted.
ans <- matrix(NA, nrow = controlGLPart$nBootUncertainty, ncol = length(theta0)
, dimnames = list(NULL, names(theta0)))
##details<<
## In addition to resampling the original data, also the temperature sensitivity
## is resampled
## from its uncertainty distribution but bounded in between 50 to 400.
E0r <- rnorm(controlGLPart$nBootUncertainty, theta0[5L], sdE_0.n )
E0 <- pmax(50, pmin(400, E0r))
theta <- theta0
#iBoot <- 1L
for (iBoot in c(1:controlGLPart$nBootUncertainty)) {
idx <- sample(nrow(dsDay), replace = TRUE)
dsDayB <- dsDay[idx, ]
theta[iPosE0] <- E0[iBoot]
resOptBoot <- try(lrcFitter$optimLRCOnAdjustedPrior(theta, iOpt = iOpt
, dsDay = dsDayB
, parameterPrior = parameterPrior, ctrl = controlGLPart), silent = TRUE)
if (!inherits(resOptBoot, "try-error") && resOptBoot$convergence == 0L) {
#TODO: also remove the very bad cases?
ans[iBoot, ]<- resOptBoot$theta
}else{
#recover()
}
}
ans
}
partGLInterpolateFluxes <- function(
### Interpolate Reco and GPP from two neighboring parameter sets
Rg ##<< numeric vector of photosynthetic flux density [umol / m2 / s]
## or Global Radiation
, VPD ##<< numeric vector of Vapor Pressure Deficit [hPa]
, Temp ##<< numeric vector of Temperature [degC]
, resParms ##<< data frame with results of \code{partGLFitLRCWindows}
## of fitting the light-response-curve for several windows
, controlGLPart = partGLControl() ##<< further default parameters,
## see \code{\link{partGLControl}}
, lrcFitter ##<< R5 class instance responsible for fitting the LRC curve
, isVerbose = TRUE ##<< set to FALSE to suppress messages
) {
##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
# corresponding weights
iValidWin <- which(is.finite(resParms$parms_out_range))
summaryLRC <- resParms %>% select(!"resOpt") %>% slice(iValidWin)
resOptList <- resParms$resOpt[iValidWin]
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
# they all will be based on the same valid data, if iMeanRec is equal
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 central 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 columns iBefore and iAfter, which can be used to merge the
# parameter estimate columns to each row
parNames <- lrcFitter$getParameterNames()
dsBefore <- merge(
structure(data.frame(dsAssoc$iSpecialBefore, dsAssoc$iBefore)
, names = c("iParRec", colNameAssoc))
, summaryLRC[, c(colNameAssoc, parNames)]
)
dsAfter <- merge(structure(data.frame(dsAssoc$iSpecialAfter, dsAssoc$iAfter)
, names = c("iParRec", colNameAssoc))
, summaryLRC[, c(colNameAssoc, lrcFitter$getParameterNames())])
if ( (nrow(dsBefore) != nRec) || (nrow(dsAfter) != nRec)) stop(#
"error in merging parameters to original records.")
Reco2 <- lapply(list(dsBefore, dsAfter), function(dsi) {
#twutz170316: fLloydTaylor gives unreasonable values with very
# low temperatures, hence constrain lower temperature
tmp <- fLloydTaylor(dsi$RRef, dsi$E0, pmax(-40, Temp_Kelvin)
, TRef = 273.15 + 15)
})
#dsi <- dsBefore
GPP2 <- lapply(list(dsBefore, dsAfter), function(dsi) {
theta <- as.matrix(dsi[, parNames])
tmp <- lrcFitter$predictLRC(theta, Rg, VPD
, Temp = pmax(-40, 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)
ans <- ansPred <- data.frame(
Reco = Reco
, GPP = GPP
)
if (isTRUE(controlGLPart$isSdPredComputed)) {
#dsi <- dsBefore
varPred2 <- lapply(list(dsBefore, dsAfter), function(dsi) {
theta <- as.matrix(dsi[, parNames])
grad <- lrcFitter$computeLRCGradient(theta, Rg, VPD, pmax(-40, 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] # row within sequence of valid
# parameters (at centralREc or meanRec)
# get the fitting object, TODO better document
resOpt <- resOptList[[iParRec]]
gradGPP <- grad$GPP[iRec, ]
gradReco <- grad$Reco[iRec, ]
# make sure parameter names match postions in covParms
varPred[iRec, 2L] <- varReco <- gradReco %*%
resOpt$covParms[names(gradReco), names(gradReco)] %*% gradReco
varPred[iRec, 1L] <- varGPP <- gradGPP %*%
resOpt$covParms[names(gradGPP), names(gradGPP)] %*% gradGPP
}
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() { # omit declaration for now
computeAggregatedCovariance <- function(
### Covariances between Reco and GPP predictions due same uncertain model coefficients
dsPred ##<< data.frame with predictors (Rg, VPD, Temp)
, resParms ##<< data.frame with results of \code{partGLFitLRCWindows} of
## fitting the light-response-curve for several windows
, resParmsNoVPD ##<< data.frame like resParms, but was fitted with option
## isNeglectVPDEffect = TRUE, for predicting if VPD is missing
, controlGLPart = partGLControl() ##<< further default parameters,
## see \code{\link{partGLControl}} with entry "isAssociateParmsToMeanOfValids"
, lrcFitter ##<< R5 class instance responsible for fitting the light response
## curve, with method getParameterNames()
, iAggregate = 1:nrow(dsPred) ##<< row indices about which to
## sum over, must be contiguous
) {
sumCovGPP <- 0
sumCovReco <- 0
# merge parameters to predictors
isValidWin <- is.finite(resParms$parms_out_range)
summaryLRC <- resParms[isValidWin, , drop = FALSE]
resOptList <- resParms$resOptList[isValidWin]
colNameAssoc <- if (isTRUE(controlGLPart$isAssociateParmsToMeanOfValids) )
"iMeanRec" else "iCentralRec"
dsAssoc <- .partGPAssociateSpecialRows(summaryLRC[[colNameAssoc]], nRec)
# now we have columns iBefore and iAfter, which can be used to
# merge the parameter estimate columns to each row
parNames <- lrcFitter$getParameterNames()
dsBefore <- merge(
structure(data.frame(dsAssoc$iSpecialBefore, dsAssoc$iBefore)
, names = c("iParRec", colNameAssoc))
, summaryLRC[, c(colNameAssoc, parNames)]
)
dsAfter <- merge(structure(data.frame(dsAssoc$iSpecialAfter, dsAssoc$iAfter)
, names = c("iParRec", colNameAssoc))
, summaryLRC[, c(colNameAssoc, lrcFitter$getParameterNames())])
#
# compute gradients (each entry is a list with components vectors Reco and GPP
grad2 <- lapply(list(dsBefore, dsAfter), function(dsi) {
theta <- as.matrix(dsi[, parNames])
grad <- lrcFitter$computeLRCGradient(theta, dsPred$Rg, dsPred$VPD
, pmax(-40, dsPred$Temp))
})
for (iRec in iAggregate[-1]) { # in first line only entry with i == j
iParRecBefore <- dsAssoc$iSpecialBefore[iRec]
resOptBefore <- resOptList[[iParRecBefore]]
iParRecAfter <- dsAssoc$iSpecialAfter[iRec]
resOptAfter <- resOptList[[iParRecAfter]]
gradGPPIBefore <- grad[[1]]$GPP[iRec, ]
gradRecoIBefore <- grad[[1]]$Reco[iRec, ]
gradGPPIAfter <- grad[[2]]$GPP[iRec, ]
gradRecoIAfter <- grad[[2]]$Reco[iRec, ]
for (jRec in iAggregate[1]:(iRec - 1) ) { # iter from start to one before iRec
# assume covariance larger than zero only in same interval between estimates
if ( (dsAssoc$iSpecialBefore[jRec] == iParRecBefore) &&
(dsAssoc$iSpecialAfter[jRec] == iParRecAfter)
) {
gradGPPJBefore <- grad[[1]]$GPP[jRec, ]
gradRecoJBefore <- grad[[1]]$Reco[jRec, ]
gradGPPJAfter <- grad[[2]]$GPP[jRec, ]
gradRecoJAfter <- grad[[2]]$Reco[jRec, ]
covRecoBefore <- gradRecoIBefore %*%
resOptBefore$covParms[names(gradRecoIBefore), names(gradRecoIBefore)] %*%
gradRecoJBefore
covGPPBefore <- gradGPPIBefore %*%
resOptBefore$covParms[names(gradGPPIBefore), names(gradGPPIBefore)] %*%
gradGPPJBefore
covRecoAfter <- gradRecoIAfter %*%
resOptAfter$covParms[names(gradRecoIAfter), names(gradRecoIAfter)] %*%
gradRecoJAfter
covGPPBAfter <- gradGPPIAfter %*%
resOptAfter$covParms[names(gradGPPIAfter), names(gradGPPIAfter)] %*%
gradGPPJAfter
covRecoMix <- dsAssoc$wBefore[iRec] * dsAssoc$wBefore[jRec] * covRecoBefore +
dsAssoc$wAfter[iRec] * dsAssoc$wAfter[jRec] * covRecoAfter
covGPPMix <- dsAssoc$wBefore[iRec] * dsAssoc$wBefore[jRec] * covGPPBefore +
dsAssoc$wAfter[iRec] * dsAssoc$wAfter[jRec] * covGPPAfter
sumCovReco <- sumCovReco + 2 * covRecoMix # symm, add also from transposed
sumCovGPP <- sumCovGPP + 2 * covGPPMix
} # match of iParRecBefore and iParRecAfter
} # for jRec
} # for iRec
##value<< named numeri vector with two components:
ans <- c(sumCovGPP = sumCovGPP ##<< sum_{i<>j} cov(GPP_i, GPP_j)
, sumCovReco = sumCovReco ##<< sum_{i<>j} cov(Reco_i, Reco_j)
)
}
} # .tmp.f
.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 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
}
replaceMissingSdByPercentage <- function(
### replace missing standard deviation of a measure x by a percentage of x
sdX ##<< numeric vector: with missing to be repalce
, x ##<< numeric vector of length(sdX): value form which percentage is computed
, perc = 0.2 ##<< numeric scalar: sdX = perc * x
, minSdX = 0.7 ##<< numeric scalar: minimum of sdX to be applied for low x
) {
##details<<
## If either perc or inSdX is NA then only the other criterion is applied.
## If both are NA then all missings are set to NA.
## \code{sdX[iToFill] <- pmax(minSdX, abs(x[iToFill] * perc), na.rm = TRUE)}
iToFill <- !is.finite(sdX)
sdX[iToFill] <- pmax(minSdX, abs(x[iToFill] * perc), na.rm = TRUE)
##value<< sdX with non-finite values replaced.
sdX
}
.computeRecoNight <- function(
### recompute ecosystem respiration at nighttime with nighttime reference
dsAnsFluxes ##<< data.frame with columsn Reco and sdReco
, dsR ##<< data.frame with columsn isNight, and Temp
, dsTempSens ##<< data.frame with column iCentralRec and E0
) {
##details<< if \code{partGLControl$useNightimeBasalRespiration == TRUE}
## then nighttime respiration is recomputed based on basal respiration
## inferred from nighttime data.
# Calculate the ecosystem respiration Reco by LlyodAndTaler
TRefK = 273.15 + 15; T0 = 227.13
RRef <- approx(
dsTempSens$iCentralRec, dsTempSens$RRef, 1:nrow(dsR), rule = 2)$y
E0 <- approx(dsTempSens$iCentralRec, dsTempSens$E0, 1:nrow(dsR), rule = 2)$y
dsAnsFluxes$Reco[dsR$isNight] <- respNight <- fLloydTaylor(
RRef[dsR$isNight], E0[dsR$isNight] , fConvertCtoK(dsR$Temp[dsR$isNight])
, TRef = TRefK)
##details<< With given uncertainty of temperature sensitivity E0, the
## logarithm of respraiton R is normally distributed, i.e. R is lognormally
## distributed. With low uncertainty (sigma*) < 1.2 this is well
## approximated by a normal distribution. Hence the sqrt(second momemt)
## of the lognormal is reported as sdReco
# log(R) = log(RRef) + tempFac*E0
sdE0 <- sqrt(approx(
dsTempSens$iCentralRec, dsTempSens$sdE0^2, 1:nrow(dsR), rule = 2)$y)
tempFacLloydTaylor <- 1/(TRefK - T0) - 1/(fConvertCtoK(dsR$Temp) - T0)
sigmaLogR = (abs(tempFacLloydTaylor)*sdE0)[dsR$isNight]
# formula of variance of lognormal with mu replaced by mean
# see lognrom/inst/docu/varianceBySigmaAndExpected.Rmd
dsAnsFluxes$sdReco[dsR$isNight] <- sqrt(exp(sigmaLogR^2) - 1)*respNight
##value<< \code{dsAnsfluxes} with updated night-time Reco and sdReco
dsAnsFluxes
}
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.