Nothing
inst/scripts/EddyPartitioning.R
In REddyProc: Data processing and plotting utilities of (half-)hourly eddy-covariance measurements
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#+++ R script with sEddyProc methods for flux partitioning +++
#+++ Flux partitionig algorithm, adapted after the PV-Wave code and paper by Markus Reichstein +++
#+++ Dependencies: Eddy.R, DataFunctions.R
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#TEST: FluxVar.s <- 'NEE_f'; QFFluxVar.s <- 'NEE_fqc'; QFFluxValue.n <- 0; TempVar.s <- 'Tair_f'; QFTempVar.s <- 'Tair_fqc'; QFTempValue.n <- 0
#TEST: RadVar.s <- 'Rg'; Lat_deg.n <- 51.0; Long_deg.n <- 13.6; TimeZone_h.n <- 1.0; CallFunction.s='test'
#TEST: NightFlux.s='FP_VARnight'; TempVar.s='FP_Temp_NEW'; WinDays.i=7; DayStep.i=5; TempRange.n=5; NumE_0.n=3; Trim.n=5
#TEST: NightFlux.s='FP_VARnight'; TempVar.s='FP_Temp_NEW'; E_0.s='E_0_NEW'; WinDays.i=4; DayStep.i=4;
#TEST: sDATA <- EddyProc.C$sDATA; sTEMP <- EddyProc.C$sTEMP
fOptimSingleE0 <- function(
##title<<
## Estimate temperature sensitivity E_0 using a Newton type optimization
##description<<
## Estimate temperature sensitivity E_0 of \code{\link{fLloydTaylor}} for a single time series
## using a Newton type optimization.
NEEnight.V.n ##<< (Original) nighttime ecosystem carbon flux, i.e. respiration vector
,Temp_degK.V.n ##<< (Original) air or soil temperature vector (degC)
,T_ref.n #=273.15+15 ##<< Reference temperature in Kelvin (degK) used in \code{fLloydTaylor} for regressing Flux and Temperature
,Trim.n=5 ##<< Percentile to trim residual (%)
,recoverOnError=FALSE ##<< Set to TRUE to debug errors instead of catching them
,algorithm="default" ##<< optimization algorithm used (see \code{\link{nls}})
)
##author<<
## AMM, TW
{
# Original implementation by AMM
res <- tryCatch({
# Non-linear regression
NLS.L <- nls(formula=R_eco ~ fLloydTaylor(R_ref, E_0, Temp, T_ref.n=T_ref.n), trace=FALSE,
data=as.data.frame(cbind(R_eco=NEEnight.V.n,Temp=Temp_degK.V.n)), start=list(R_ref=2,E_0=200)
,algorithm=algorithm
)
# Remove points with residuals outside Trim.n quantiles
Residuals.V.n <- resid(NLS.L)
#Residuals.V.n <- fLloydTaylor(R_ref=coef(summary(NLS.L))['R_ref',1], E_0=coef(summary(NLS.L))['E_0',1],
# Temp_degK.V.n, T_ref.n=T_ref.n) - NEEnight.V.n
t.b <- Residuals.V.n >= quantile(Residuals.V.n, probs=c(Trim.n/100)) & Residuals.V.n <= quantile(Residuals.V.n, probs=c(1-Trim.n/100))
# Trimmed non-linear regression
NLS_trim.L <- nls(formula=R_eco ~ fLloydTaylor(R_ref, E_0, Temp, T_ref.n=T_ref.n), algorithm='default', trace=FALSE,
data=as.data.frame(cbind(R_eco=NEEnight.V.n[t.b], Temp=Temp_degK.V.n[t.b])), start=list(R_ref=2,E_0=200))
##value<< Numeric vector with following components:
res <- c(
R_ref=coef(summary(NLS.L))['R_ref',1] ##<< Estimate of respiration rate at reference temperature
, R_ref_SD=coef(summary(NLS.L))['R_ref',2] ##<< Standard deviation of R_ref
, E_0=coef(summary(NLS.L))['E_0',1] ##<< Estimate of temperature sensitivity ("activation energy") in Kelvin (degK) for untrimmed dataset
, E_0_SD=coef(summary(NLS.L))['E_0',2] ##<< Standard deviation of E_0
, E_0_trim=coef(summary(NLS_trim.L))['E_0',1] ##<< Estimate of temperature sensitivity ("activation energy") in Kelvin (degK) for trimmed dataset
, E_0_trim_SD=coef(summary(NLS_trim.L))['E_0',2] ##<< Standard deviation of E_0_trim
##end<<
)
# Note on other tested algorithms:
# Standard require(stats) nls with PORT algorithm and lower and upper bounds
# require(FME) for modFit and modCost, has PORT algorithm included (and other algorithms like MCMC)
# require(robustbase) for ltsReg but only linear regression
# require(nlme) for heteroscedastic and mixed NLR but no port algo with upper and lower bounds
# require(nlstools) for bootstrapping with nlsBoot(nls...)
}, error = function(e) {
if( isTRUE(recoverOnError) ) recover()
res <- c(R_ref=NA, R_ref_SD=NA, E_0=NA, E_0_SD=NA, E_0_trim=NA, E_0_trim_SD=NA)
} ) #Spaces between brackets required to avoid replacement on documentation generation
res
}
fOptimSingleE0_Lev <- function(
##title<<
## Estimate temperature sensitivity E_0 using Levenberg-Marquard optimization
##description<<
## Estimate temperature sensitivity E_0 of \code{\link{fLloydTaylor}} for a single time series
## using Levenberg-Marquard optimization.
NEEnight.V.n ##<< (Original) nighttime ecosystem carbon flux, i.e. respiration vector
,Temp_degK.V.n ##<< (Original) air or soil temperature vector (degC)
,T_ref.n #=273.15+15 ##<< Reference temperature in Kelvin (degK) used in \code{fLloydTaylor} for regressing Flux and Temperature
,Trim.n=5 ##<< Percentile to trim residual (%)
,recoverOnError=FALSE ##<< Set to TRUE to debug errors instead of catching them
,algorithm='LM' ##<< optimization algorithm used (see nlsLM from package minpack.lm)
)
##author<<
## TW
{
if( !requireNamespace('minpack.lm') ) stop("Need to install package minpack.lm before using LM optimization.")
res <- tryCatch({
# Non-linear regression
NLS.L <- minpack.lm::nlsLM(formula=R_eco ~ fLloydTaylor(R_ref, E_0, Temp, T_ref.n=T_ref.n), trace=FALSE,
data=as.data.frame(cbind(R_eco=NEEnight.V.n,Temp=Temp_degK.V.n)), start=list(R_ref=2,E_0=200)
,control=minpack.lm::nls.lm.control(maxiter = 20)
,algorithm=algorithm
)
# Remove points with residuals outside Trim.n quantiles
Residuals.V.n <- resid(NLS.L)
#plot( Residuals.V.n ~ Temp_degK.V.n )
t.b <- Residuals.V.n >= quantile(Residuals.V.n, probs=c(Trim.n/100)) & Residuals.V.n <= quantile(Residuals.V.n, probs=c(1-Trim.n/100))
#points( Residuals.V.n[!t.b] ~ Temp_degK.V.n[!t.b], col="red")
# Trimmed non-linear regression
NLS_trim.L <- minpack.lm::nlsLM(formula=R_eco ~ fLloydTaylor(R_ref, E_0, Temp, T_ref.n=T_ref.n), algorithm='default', trace=FALSE,
data=as.data.frame(cbind(R_eco=NEEnight.V.n[t.b], Temp=Temp_degK.V.n[t.b]))
, start=coef(NLS.L)
#, start=list(R_ref=2,E_0=200)
)
##value<< Numeric vector with following components:
res <- c(
R_ref=coef(summary(NLS.L))['R_ref',1] ##<< Estimate of espiration rate at reference temperature
, R_ref_SD=coef(summary(NLS.L))['R_ref',2] ##<< Standard deviation of R_ref
,E_0=coef(summary(NLS.L))['E_0',1] ##<< Estimate of temperature sensitivity ("activation energy") in Kelvin (degK) for untrimmed dataset
, E_0_SD=coef(summary(NLS.L))['E_0',2] ##<< Standard deviation of E_0
,E_0_trim=coef(summary(NLS_trim.L))['E_0',1] ##<< Estimate of temperature sensitivity ("activation energy") in Kelvin (degK) for trimmed dataset
, E_0_trim_SD=coef(summary(NLS_trim.L))['E_0',2] ##<< Standard deviation of E_0_trim
##end<<
)
}, error = function(e) {
if( isTRUE(recoverOnError) ) recover()
res <- c(R_ref=NA, R_ref_SD=NA, E_0=NA, E_0_SD=NA, E_0_trim=NA, E_0_trim_SD=NA)
} ) #Spaces between brackets required to avoid replacement on documentation generation
res
}
fRegrE0fromShortTerm = function(
##title<<
## Estimation of the temperature sensitivity E_0
##description<<
## Estimation of the temperature sensitivity E_0 from regression of \code{\link{fLloydTaylor}} for short periods
NightFlux.V.n ##<< (Original) nighttime ecosystem carbon flux, i.e. respiration vector
,TempVar.V.n ##<< (Original) air or soil temperature vector (degC)
,DayCounter.V.i ##<< Integer specifying the day of each record
,WinDays.i=7 ##<< Window size for \code{\link{fLloydTaylor}} regression in days
,DayStep.i=5 ##<< Window step for \code{\link{fLloydTaylor}} regression in days
,TempRange.n=5 ##<< Threshold temperature range to start regression (#! Could be larger for Tair)
,Trim.n=5 ##<< Percentile to trim residual (%)
,NumE_0.n=3 ##<< Number of best E_0's to average over
,MinE_0.n=30 ##<< Minimum E0 for validity check
,MaxE_0.n=450 ##<< Maximum E0 for validity check
,T_ref.n #=273.15+15 ##<< Reference temperature in Kelvin (degK) used in \code{fLloydTaylor} for regressing Flux and Temperature
,CallFunction.s='' ##<< Name of function called from
,optimAlgorithm='default' ##<< optimization algorithm used (see \code{\link{nls}} ) or 'LM' for Levenberg-Marquard (see nlsLM from package minpack.lm)
)
##author<<
## AMM
{
##details<<
## The coefficient E0 is estimated for windows with a length of \code{WinDays.i} days,
## for successive periods in steps of \code{DayStep.i} days.
## Only those windows with a sufficient number or records and with a sufficient temperature range \code{TempRange.n}
## are used for the \code{\link{fLloydTaylor}} regression of E0 using the optimization \code{\link{fOptimSingleE0}}.
## Unreasonable estimates are discarded (95% confidence interval inside \code{MinE_0.n} and \code{MaxE_0.n}) and
## the others are ordered by their standard deviations.
## The mean across the best (=lowest standard deviation) E0 estimates is reported
## with \code{NumE_0.n} defining the number of best estimates to use.
# Regression settings
#NLSRes.F <- data.frame(NULL) #Results of non-linear regression
#NLSRes_trim.F <- data.frame(NULL) #Results of non-linear regression
MinData.n <- 6 # Minimum number of data points
fOptim <- fOptimSingleE0
if( optimAlgorithm=='LM') fOptim <- fOptimSingleE0_Lev
#tw: better use rbind with a list instead of costly repeated extending a data.frame
NLSRes.F <- as.data.frame(do.call( rbind, NLSRes.l <- lapply( seq(WinDays.i+1, max(DayCounter.V.i), DayStep.i) ,function(DayMiddle.i){
#TEST: DayMiddle.i <- 8
DayStart.i <- DayMiddle.i-WinDays.i
DayEnd.i <- DayMiddle.i+WinDays.i
#! Window size of 7 days corresponds to full window length of 15 days as in paper, non-congruent with PV-Wave code of 14 days
#! New code: Last window has minimum width of WinDays.i
Subset.b <- DayCounter.V.i >= DayStart.i & DayCounter.V.i <= DayEnd.i & !is.na(NightFlux.V.n) & !is.na(TempVar.V.n)
NEEnight.V.n <- subset(NightFlux.V.n, Subset.b)
Temp.V.n <- subset(TempVar.V.n, Subset.b)
Temp_degK.V.n <- fConvertCtoK(Temp.V.n)
if( length(NEEnight.V.n) > MinData.n && diff(range(Temp_degK.V.n)) >= TempRange.n ) {
#CountRegr.i <- CountRegr.i+1
resOptim <- fOptim( NEEnight.V.n, Temp_degK.V.n, algorithm=optimAlgorithm, T_ref.n=T_ref.n)
NLSRes.F <- c(Start=DayStart.i, End=DayEnd.i, Num=length(NEEnight.V.n), TRange=diff(range(Temp_degK.V.n)),
resOptim)
} else NULL
} ) ))
Limits.b <- ( NLSRes.F$E_0_trim - NLSRes.F$E_0_trim_SD > MinE_0.n & NLSRes.F$E_0_trim + NLSRes.F$E_0_trim_SD < MaxE_0.n )
#! New code: Check validity with SD (standard deviation) limits, in PV-Wave without SD, in paper if E_0_SD < (E_0 * 50%)
NLSRes.F$E_0_trim_ok <- ifelse( Limits.b, NLSRes.F$E_0_trim, NA)
NLSRes.F$E_0_trim_SD_ok <- ifelse( Limits.b, NLSRes.F$E_0_trim_SD, NA)
#
# Sort data frame for smallest standard deviation
NLSsort.F <- NLSRes.F[order(NLSRes.F$E_0_trim_SD_ok),] # ordered data.frame
##details<<
## Take average of the three E_0 with lowest standard deviation
E_0_trim.n <- round(mean(NLSsort.F$E_0_trim_ok[1:NumE_0.n]), digits=2)
#
# Abort flux partitioning if regression of E_0 failed
if( is.na(E_0_trim.n) ) {
# twutz 150226: just warning and returning negative value gives problems later on, maybe better stop and catch exception
warning(CallFunction.s, ':::fRegrE0fromShortTerm::: Less than ', NumE_0.n, ' valid values for E_0 after regressing ',
nrow(NLSRes.F), ' periods! Aborting partitioning.')
return(-111)
}
E_0_trim.n
##value<<
## Estimated scalar temperature sensitivity (E_0, degK)
}
sEddyProc$methods(
sRegrE0fromShortTerm = function(
##title<<
## sEddyProc$sRegrE0fromShortTerm - Estimation of the temperature sensitivity E_0
##description<<
## Estimation of the temperature sensitivity E_0 from regression of \code{\link{fLloydTaylor}}
## for short periods by calling \code{\link{fRegrE0fromShortTerm}}
NightFlux.s ##<< Variable with (original) nighttime ecosystem carbon flux, i.e. respiration
,TempVar.s ##<< Variable with (original) air or soil temperature (degC)
,... ##<< Parameters passed to \code{\link{fRegrE0fromShortTerm}}
,CallFunction.s='' ##<< Name of function called from
,debug.l = list(fixedE0=NA) ##<< List with controls for debugging, see details
)
##author<<
## AMM, TW
{
'Estimation of the temperature sensitivity E_0 from regression of fLloydTaylor() for short periods by calling fRegrE0fromShortTerm()'
##details<< For further details see \code{\link{fRegrE0fromShortTerm}}.
#
# Check if specified columns are numeric
SubCallFunc.s <- paste(CallFunction.s,'sRegrE0fromShortTerm', sep=':::')
ds <- cbind(sDATA,sTEMP)
fCheckColNames(ds, c(NightFlux.s, TempVar.s), SubCallFunc.s)
fCheckColNum(ds, c(NightFlux.s, TempVar.s), SubCallFunc.s)
#
##details<< \describe{ \item{Debugging control}{
## When supplying a finite scalar value \code{debug.l$fixedE0}, then this value
## is used instead of the temperature sensitivity E_0 from short term data.
## }}
if( (length(debug.l$fixedE0) != 0) && is.finite(debug.l$fixedE0) ){
E_0_trim.n <- debug.l$fixedE0[1]
message('Using prescribed temperature sensitivity E_0 of: ', format(E_0_trim.n, digits=5), '.')
}else{
# Write into vectors since cbind of data frames is so slow (factor of ~20 (!))
NightFlux.V.n <- ds[,NightFlux.s]
TempVar.V.n <- ds[,TempVar.s]
DayCounter.V.i <- c(1:sINFO$DIMS) %/% sINFO$DTS
#
# Check for validity of E_0 regression results
if( grepl('Tair', TempVar.s) ) {
#Limits in PV-Wave code for Tair
MinE_0.n <- 30; MaxE_0.n <- 350
} else if( grepl('Tsoil', TempVar.s) ) {
#Limits in PV-Wave code for Tsoil
MinE_0.n <- 30; MaxE_0.n <- 550 # Higher values due to potentially high Q10 values
} else {
#Default limits taken from paper
MinE_0.n <- 30; MaxE_0.n <- 450
}
#
E_0_trim.n <- fRegrE0fromShortTerm( NightFlux.V.n, TempVar.V.n, DayCounter.V.i, ...,
MinE_0.n=MinE_0.n, MaxE_0.n=MaxE_0.n, CallFunction.s=SubCallFunc.s)
message('Estimate of the temperature sensitivity E_0 from short term data: ', format(E_0_trim.n, digits=5), '.')
}
# Add constant value of E_0 as column vector to sTEMP
E_0.V.n <- rep(E_0_trim.n, nrow(sTEMP))
attr(E_0.V.n, 'varnames') <- 'E_0'
attr(E_0.V.n, 'units') <- 'degK'
E_0.V.n
##value<<
## Data vector with (constant) temperature sensitivity (E_0, degK)
## On failure, all entries are set to -111
})
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
sEddyProc$methods(
sRegrRref = function(
##title<<
## sEddyProc$sRegrRref - Estimation of the time series of reference respiration Rref
##description<<
## Estimation of the reference respiration Rref of \code{\link{fLloydTaylor}} for successive periods
NightFlux.s ##<< Variable with (original) nighttime ecosystem carbon flux, i.e. respiration
,TempVar.s ##<< Variable with (original) air or soil temperature (degC)
,E_0.s ##<< Temperature sensitivity E_0 estimated with \code{\link{sRegrE0fromShortTerm}}
,T_ref.n #=273.15+15 ##<< Reference temperature in Kelvin (degK) used in \code{fLloydTaylor} for regressing Flux and Temperature
,WinDays.i=3 ##<< Window size for \code{\link{fLloydTaylor}} regression in days
,DayStep.i=4 ##<< Window step for \code{\link{fLloydTaylor}} regression in days
,CallFunction.s='' ##<< Name of function called from
)
##author<<
## AMM
##details<<
## While parameter E0 in the Temperature-Respiration relationship (\code{\link{fLloydTaylor}}) is kept konstant,
## parameter Rref is allowed to change with time.
## This method estimates Rref for a series of time windows of length 2*\code{WinDays.i}+1 days
## shifted by \code{DayStep.i} days.
##
## For some of the windows, it maybe not be possible to estimate Rref. These missing values are filled by linear
## interpolation by function \code{\link{fInterpolateGaps}}.
{
'Estimation of the reference respiration Rref of fLloydTaylor() for successive periods'
# Check if specified columns are numeric
SubCallFunc.s <- paste(CallFunction.s,'sRegrRref', sep=':::')
ds <- cbind(sDATA,sTEMP)
fCheckColNames(ds, c(NightFlux.s, TempVar.s, E_0.s), SubCallFunc.s)
fCheckColNum(ds, c(NightFlux.s, TempVar.s, E_0.s), SubCallFunc.s)
# Regression settings
LMRes.F <- data.frame(NULL) #Results of linear regression
MinData.n <- 2 # Minimum number of data points for regression
# Write into vectors since cbind of data frames is so slow (factor of ~20 (!))
NightFlux.V.n <- ds[,NightFlux.s]
TempVar.V.n <- ds[,TempVar.s]
# Loop regression periods
DayCounter.V.i <- c(1:sINFO$DIMS) %/% sINFO$DTS
CountRegr.i <- 0
for (DayMiddle.i in seq(WinDays.i+1, max(DayCounter.V.i), DayStep.i)) {
#TEST: DayMiddle.i <- 8
DayStart.i <- DayMiddle.i-WinDays.i
DayEnd.i <- DayMiddle.i+WinDays.i
#! Window size of 4 days corresponds to a full window length of 9 days, non-congruent with PV-Wave code of 8 days, in paper not mentioned
#! New code: Last window has minimum of window size
Subset.b <- DayCounter.V.i >= DayStart.i & DayCounter.V.i <= DayEnd.i & !is.na(NightFlux.V.n) & !is.na(TempVar.V.n)
MeanHour.i <- round(mean(which(Subset.b))) # Weighted middle of the time period
NEEnight.V.n <- subset(NightFlux.V.n, Subset.b)
Temp.V.n <- subset(TempVar.V.n, Subset.b)
Temp_degK.V.n <- fConvertCtoK(Temp.V.n)
E_0.V.n <- subset(sTEMP[,E_0.s], Subset.b) # (Constant value)
if( length(NEEnight.V.n) > MinData.n ) {
CountRegr.i <- CountRegr.i+1
tryCatch({
LM.L <- lm(R_eco ~ 0 + fLloydTaylor(R_ref, E_0, Temp_degK, T_ref.n=T_ref.n), data=as.data.frame(cbind(R_eco=NEEnight.V.n, R_ref=1, E_0=E_0.V.n, Temp_degK=Temp_degK.V.n)))
LMRes.F <- rbind(LMRes.F, cbind(Start=DayStart.i, End=DayEnd.i, Num=length(NEEnight.V.n), MeanH=MeanHour.i,
R_ref=coef(summary(LM.L))[1], R_ref_SD=coef(summary(LM.L))[2]))
#! Note on PV-Wave code: trimmed linear regression not used in the end, i.e. in online webtool
if( F ) { # Plot for testing
plot(NEEnight.V.n ~ fLloydTaylor(1, E_0.V.n, Temp_degK.V.n, T_ref.n=T_ref.n))
curve(coef(LM.L)[1] * x, add=T, col='green')
}
}, error = function(e) {
LMRes.F <- rbind(LMRes.F, cbind(Start=DayStart.i, End=DayEnd.i, Num=length(NEEnight.V.n), MeanH=MeanHour.i,
R_ref=NA, R_ref_SD=NA))
} ) #Spaces between brackets required to avoid replacement on documentation generation
}
}
# Check for validity of regression results
LMRes.F$R_ref_ok <- ifelse( LMRes.F$R_ref < 0 , NA, LMRes.F$R_ref )
#! New code: Omit regressions with R_ref <0, in PV-Wave smaller values are set to 0.000001, not mentioned in paper
#TODO later: Flag for long distances between R_refs, especially if long distance in the beginning
#TODO later: Provide some kind of uncertainty estimate from R_ref_SD
# Interpolate R_ref periods linearly between MeanHour.i and keep constant at beginning/end
Rref.V.n <- rep(NA, length(NightFlux.V.n))
Rref.V.n[LMRes.F$MeanH] <- LMRes.F$R_ref_ok
Rref.V.n <- fInterpolateGaps(Rref.V.n)
attr(Rref.V.n, 'varnames') <- 'R_ref'
attr(Rref.V.n, 'units') <- attr(sTEMP[,NightFlux.s], 'units')
message('Regression of reference temperature R_ref for ', sum(!is.na(LMRes.F$R_ref_ok)), ' periods.')
Rref.V.n
##value<<
## Data vector (length number of windows) with reference respiration (R_ref, flux units)
})
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
sEddyProc$methods(
sMRFluxPartition = function(
##title<<
## sEddyProc$sMRFluxPartition - Flux partitioning after Reichstein et al. (2005)
##description<<
## Nighttime-based partitioning of measured net ecosystem fluxes into gross primary production (GPP) and ecosystem respiration (Reco)
FluxVar.s='NEE_f' ##<< Variable name of column with original and filled net ecosystem fluxes (NEE)
,QFFluxVar.s='NEE_fqc' ##<< Quality flag of NEE variable
,QFFluxValue.n=0 ##<< Value of quality flag for _good_ (original) data
,TempVar.s='Tair_f' ##<< Filled air- or soil temperature variable (degC)
,QFTempVar.s='Tair_fqc'##<< Quality flag of filled temperature variable
,QFTempValue.n=0 ##<< Value of temperature quality flag for _good_ (original) data
,RadVar.s='Rg' ##<< Unfilled (original) radiation variable
,Lat_deg.n ##<< Latitude in (decimal) degrees
,Long_deg.n ##<< Longitude in (decimal) degrees
,TimeZone_h.n ##<< Time zone (in hours)
,T_ref.n=273.15+15 ##<< Reference temperature in Kelvin (degK) used in \code{fLloydTaylor} for regressing Flux and Temperature
,Suffix.s = '' ##<< String suffix needed for different processing setups on the same dataset (for explanations see below)
,debug.l=list( ##<< List with debugging control (passed also to \code{\link{sRegrE0fromShortTerm}}).
##describe<<
useLocaltime.b=FALSE ##<< see details on solar vs local time
##end<<
)
,parsE0Regression=list() ##<< list with further parameters passed down to \code{\link{sRegrE0fromShortTerm}} and \code{\link{fRegrE0fromShortTerm}}, such as \code{TempRange.n}
)
##author<<
## AMM,TW
##references<<
## Reichstein M, Falge E, Baldocchi D et al. (2005) On the separation of net ecosystem exchange
## into assimilation and ecosystem respiration: review and improved algorithm. Global Change Biology, 11, 1424-1439.
{
'Partitioning of measured net ecosystem fluxes into gross primary production (GPP) and ecosystem respiration (Reco)'
##details<< \describe{\item{Description of newly generated variables with partitioning results:}{
## \itemize{
## \item PotRad - Potential radiation \cr
## \item FP_NEEnight - Good (original) NEE nighttime fluxes used for flux partitioning \cr
## \item FP_Temp - Good (original) temperature measurements used for flux partitioning \cr
## \item E_0 - Estimated temperature sensitivity \cr
## \item R_ref - Estimated reference respiration \cr
## \item Reco - Estimated ecosystem respiration \cr
## \item GPP_f - Estimated gross primary production \cr
## }
## }}
##details<< \describe{\item{Background}{
## This partitioning is based on the regression of nighttime respiration with temperature
## using the Lloyd-Taylor-Function \code{\link{fLloydTaylor}}.
## First the temperature sensitivity E_0 is estimated from short term data, see \code{\link{sRegrE0fromShortTerm}}.
## Next the reference temperature R_ref is estimated for successive periods throughout the whole dataset (see \code{\link{sRegrRref}}).
## These estimates are then used to calculate the respiration during daytime and nighttime and with this GPP.
## Attention: Gap filling of the net ecosystem fluxes (NEE) and temperature measurements (Tair or Tsoil) is required
## prior to the partitioning!
## }}
# If the default variable names are used and a Suffix.s is provided, then the suffix is added to the variable name (see also comment below)
if( FluxVar.s=='NEE_f' && QFFluxVar.s=='NEE_fqc' && fCheckValString(Suffix.s) ) {
FluxVar.s <- paste('NEE_', Suffix.s, '_f', sep='')
QFFluxVar.s <- paste('NEE_', Suffix.s, '_fqc', sep='')
}
ds <- cbind(sDATA,sTEMP)
# Check if specified columns exist in sDATA or sTEMP and if numeric and plausible. Then apply quality flag
fCheckColNames(ds, c(FluxVar.s, QFFluxVar.s, TempVar.s, QFTempVar.s, RadVar.s), 'sMRFluxPartition')
fCheckColNum(ds, c(FluxVar.s, QFFluxVar.s, TempVar.s, QFTempVar.s, RadVar.s), 'sMRFluxPartition')
fCheckColPlausibility(ds, c(FluxVar.s, QFFluxVar.s, TempVar.s, QFTempVar.s, RadVar.s), 'sMRFluxPartition')
Var.V.n <- fSetQF(ds, FluxVar.s, QFFluxVar.s, QFFluxValue.n, 'sMRFluxPartition')
message('Start flux partitioning for variable ', FluxVar.s, ' with temperature ', TempVar.s, '.')
##details<< \describe{\item{Selection of daytime data based on solar time}{
## The respiration-temperature regression is very
## sensitive to the selection of night- and daytime data.
## Nighttime is selected by a combined threshold of current solar radiation and potential radiation.
## The current implementation calculates potential radiation based on exact solar time, based on latitude and longitude.
## (see \code{\link{fCalcPotRadiation}} )
## Therefore it might differ from implementations that use local winter clock time instead.
## }}
# Calculate potential radiation
#! New code: Local time and equation of time accounted for in potential radiation calculation
DoY.V.n <- as.numeric(format(sDATA$sDateTime, '%j'))
Hour.V.n <- as.numeric(format(sDATA$sDateTime, '%H')) + as.numeric(format(sDATA$sDateTime, '%M'))/60
sTEMP$PotRad_NEW <<- fCalcPotRadiation(DoY.V.n, Hour.V.n, Lat_deg.n, Long_deg.n, TimeZone_h.n
, useSolartime.b=!isTRUE(debug.l$useLocaltime.b) )
# Filter night time values only
#! Note: Rg <= 10 congruent with MR PV-Wave, in paper Rg <= 20
# Should be unfilled (original) radiation variable, therefore dataframe set to sDATA only
sTEMP$FP_VARnight <<- ifelse(sDATA[,RadVar.s] > 10 | sTEMP$PotRad_NEW > 0, NA, Var.V.n)
attr(sTEMP$FP_VARnight, 'varnames') <<- paste(attr(Var.V.n, 'varnames'), '_night', sep='')
attr(sTEMP$FP_VARnight, 'units') <<- attr(Var.V.n, 'units')
#! New code: Slightly different subset than PV-Wave due to time zone correction (avoids timezone offset between Rg and PotRad)
# Apply quality flag for temperature
sTEMP$FP_Temp_NEW <<- fSetQF(ds, TempVar.s, QFTempVar.s, QFTempValue.n, 'sMRFluxPartition')
# Estimate E_0 and R_ref (results are saved in sTEMP)
# twutz1508: changed to do.call in order to allow passing further parameters when calling sMRFluxPartition
sTEMP$E_0_NEW <<- do.call( .self$sRegrE0fromShortTerm, c(
list('FP_VARnight', 'FP_Temp_NEW', T_ref.n=T_ref.n, CallFunction.s='sMRFluxPartition', debug.l=debug.l)
, parsE0Regression))
if( sum(sTEMP$E_0_NEW==-111) != 0 )
return(invisible(-111)) # Abort flux partitioning if regression of E_0 failed
# Reanalyse R_ref with E_0 fixed
sTEMP$R_ref_NEW <<- sRegrRref('FP_VARnight', 'FP_Temp_NEW', 'E_0_NEW', T_ref.n=T_ref.n, CallFunction.s='sMRFluxPartition')
# Calculate the ecosystem respiration Reco
sTEMP$Reco_NEW <<- fLloydTaylor(sTEMP$R_ref_NEW, sTEMP$E_0_NEW, fConvertCtoK(ds[,TempVar.s]), T_ref.n=T_ref.n)
attr(sTEMP$Reco_NEW, 'varnames') <<- 'Reco'
attr(sTEMP$Reco_NEW, 'units') <<- attr(Var.V.n, 'units')
# Calculate the gross primary production GPP_f
sTEMP$GPP_NEW_f <<- -ds[,FluxVar.s] + sTEMP$Reco_NEW
sTEMP$GPP_NEW_fqc <<- ds[,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)
attr(sTEMP$GPP_NEW_f, 'varnames') <<- 'GPP_f'
attr(sTEMP$GPP_NEW_f, 'units') <<- attr(Var.V.n, 'units')
##details<< \describe{\item{Different processing setups on the same dataset}{
## Attention: When processing the same site data set with different setups for the gap filling or flux partitioning
## (e.g. due to different ustar filters),
## a string suffix is needed! This suffix is added to the result column names to distinguish the results of the different setups.
## If a Suffix.s is provided and if the defaults for FluxVar.s and QFFluxVar.s are used, the Suffix.s will be added to their variable names
## (e.g. 'NEE_f' will be renamed to 'NEE_WithUstar_f' and 'NEE_fqc' to 'NEE_WithUstar_fqc' for the Suffix.s='WithUstar').
## Currently, this works only with defaults of FluxVar.s='NEE_f' and QFFluxVar.s='NEE_fqc'.
## }}
# Rename new columns generated during flux partitioning:
# For nighttime NEE (FP_NEEnight or FP_NEEnight_Suffix)
colnames(sTEMP) <<- gsub('_VARnight', paste('_NEEnight', (if(fCheckValString(Suffix.s)) '_' else ''), Suffix.s, sep=''), colnames(sTEMP))
# For the results columns, the _NEW is dropped and the suffix added
colnames(sTEMP) <<- gsub('_NEW', paste((if(fCheckValString(Suffix.s)) '_' else ''), Suffix.s, sep=''), colnames(sTEMP))
# Check for duplicate columns (to detect if different processing setups were executed without different suffixes)
if( length(names(iDupl <- which(table(colnames(sTEMP)) > 1))) ) {
warning(paste0('sMRFluxPartition::: Duplicated columns found! (',
paste(names(iDupl),collapse=",")
,') Deleting each first of duplicate columns.'
,' Please use different Suffix.s when processing different setups on the same dataset!'))
for( cname in names(iDupl) ) sTEMP[cname] <<- NULL # need to remove columns else some tests fail
}
return(invisible(NULL))
##value<<
## Flux partitioning results (see variables in details) in sTEMP data frame (with renamed columns).
## On success, return value is NULL. On failure an integer scalar error code is returned:
## -111 if regression of E_0 failed due to insufficient relationship in the data.
})
Try the REddyProc package in your browser
Any scripts or data that you put into this service are public.
REddyProc documentation built on May 2, 2019, 5:19 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#+++ R script with sEddyProc methods for flux partitioning +++
#+++ Flux partitionig algorithm, adapted after the PV-Wave code and paper by Markus Reichstein +++
#+++ Dependencies: Eddy.R, DataFunctions.R
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#TEST: FluxVar.s <- 'NEE_f'; QFFluxVar.s <- 'NEE_fqc'; QFFluxValue.n <- 0; TempVar.s <- 'Tair_f'; QFTempVar.s <- 'Tair_fqc'; QFTempValue.n <- 0
#TEST: RadVar.s <- 'Rg'; Lat_deg.n <- 51.0; Long_deg.n <- 13.6; TimeZone_h.n <- 1.0; CallFunction.s='test'
#TEST: NightFlux.s='FP_VARnight'; TempVar.s='FP_Temp_NEW'; WinDays.i=7; DayStep.i=5; TempRange.n=5; NumE_0.n=3; Trim.n=5
#TEST: NightFlux.s='FP_VARnight'; TempVar.s='FP_Temp_NEW'; E_0.s='E_0_NEW'; WinDays.i=4; DayStep.i=4;
#TEST: sDATA <- EddyProc.C$sDATA; sTEMP <- EddyProc.C$sTEMP
fOptimSingleE0 <- function(
##title<<
## Estimate temperature sensitivity E_0 using a Newton type optimization
##description<<
## Estimate temperature sensitivity E_0 of \code{\link{fLloydTaylor}} for a single time series
## using a Newton type optimization.
NEEnight.V.n ##<< (Original) nighttime ecosystem carbon flux, i.e. respiration vector
,Temp_degK.V.n ##<< (Original) air or soil temperature vector (degC)
,T_ref.n #=273.15+15 ##<< Reference temperature in Kelvin (degK) used in \code{fLloydTaylor} for regressing Flux and Temperature
,Trim.n=5 ##<< Percentile to trim residual (%)
,recoverOnError=FALSE ##<< Set to TRUE to debug errors instead of catching them
,algorithm="default" ##<< optimization algorithm used (see \code{\link{nls}})
)
##author<<
## AMM, TW
{
# Original implementation by AMM
res <- tryCatch({
# Non-linear regression
NLS.L <- nls(formula=R_eco ~ fLloydTaylor(R_ref, E_0, Temp, T_ref.n=T_ref.n), trace=FALSE,
data=as.data.frame(cbind(R_eco=NEEnight.V.n,Temp=Temp_degK.V.n)), start=list(R_ref=2,E_0=200)
,algorithm=algorithm
)
# Remove points with residuals outside Trim.n quantiles
Residuals.V.n <- resid(NLS.L)
#Residuals.V.n <- fLloydTaylor(R_ref=coef(summary(NLS.L))['R_ref',1], E_0=coef(summary(NLS.L))['E_0',1],
# Temp_degK.V.n, T_ref.n=T_ref.n) - NEEnight.V.n
t.b <- Residuals.V.n >= quantile(Residuals.V.n, probs=c(Trim.n/100)) & Residuals.V.n <= quantile(Residuals.V.n, probs=c(1-Trim.n/100))
# Trimmed non-linear regression
NLS_trim.L <- nls(formula=R_eco ~ fLloydTaylor(R_ref, E_0, Temp, T_ref.n=T_ref.n), algorithm='default', trace=FALSE,
data=as.data.frame(cbind(R_eco=NEEnight.V.n[t.b], Temp=Temp_degK.V.n[t.b])), start=list(R_ref=2,E_0=200))
##value<< Numeric vector with following components:
res <- c(
R_ref=coef(summary(NLS.L))['R_ref',1] ##<< Estimate of respiration rate at reference temperature
, R_ref_SD=coef(summary(NLS.L))['R_ref',2] ##<< Standard deviation of R_ref
, E_0=coef(summary(NLS.L))['E_0',1] ##<< Estimate of temperature sensitivity ("activation energy") in Kelvin (degK) for untrimmed dataset
, E_0_SD=coef(summary(NLS.L))['E_0',2] ##<< Standard deviation of E_0
, E_0_trim=coef(summary(NLS_trim.L))['E_0',1] ##<< Estimate of temperature sensitivity ("activation energy") in Kelvin (degK) for trimmed dataset
, E_0_trim_SD=coef(summary(NLS_trim.L))['E_0',2] ##<< Standard deviation of E_0_trim
##end<<
)
# Note on other tested algorithms:
# Standard require(stats) nls with PORT algorithm and lower and upper bounds
# require(FME) for modFit and modCost, has PORT algorithm included (and other algorithms like MCMC)
# require(robustbase) for ltsReg but only linear regression
# require(nlme) for heteroscedastic and mixed NLR but no port algo with upper and lower bounds
# require(nlstools) for bootstrapping with nlsBoot(nls...)
}, error = function(e) {
if( isTRUE(recoverOnError) ) recover()
res <- c(R_ref=NA, R_ref_SD=NA, E_0=NA, E_0_SD=NA, E_0_trim=NA, E_0_trim_SD=NA)
} ) #Spaces between brackets required to avoid replacement on documentation generation
res
}
fOptimSingleE0_Lev <- function(
##title<<
## Estimate temperature sensitivity E_0 using Levenberg-Marquard optimization
##description<<
## Estimate temperature sensitivity E_0 of \code{\link{fLloydTaylor}} for a single time series
## using Levenberg-Marquard optimization.
NEEnight.V.n ##<< (Original) nighttime ecosystem carbon flux, i.e. respiration vector
,Temp_degK.V.n ##<< (Original) air or soil temperature vector (degC)
,T_ref.n #=273.15+15 ##<< Reference temperature in Kelvin (degK) used in \code{fLloydTaylor} for regressing Flux and Temperature
,Trim.n=5 ##<< Percentile to trim residual (%)
,recoverOnError=FALSE ##<< Set to TRUE to debug errors instead of catching them
,algorithm='LM' ##<< optimization algorithm used (see nlsLM from package minpack.lm)
)
##author<<
## TW
{
if( !requireNamespace('minpack.lm') ) stop("Need to install package minpack.lm before using LM optimization.")
res <- tryCatch({
# Non-linear regression
NLS.L <- minpack.lm::nlsLM(formula=R_eco ~ fLloydTaylor(R_ref, E_0, Temp, T_ref.n=T_ref.n), trace=FALSE,
data=as.data.frame(cbind(R_eco=NEEnight.V.n,Temp=Temp_degK.V.n)), start=list(R_ref=2,E_0=200)
,control=minpack.lm::nls.lm.control(maxiter = 20)
,algorithm=algorithm
)
# Remove points with residuals outside Trim.n quantiles
Residuals.V.n <- resid(NLS.L)
#plot( Residuals.V.n ~ Temp_degK.V.n )
t.b <- Residuals.V.n >= quantile(Residuals.V.n, probs=c(Trim.n/100)) & Residuals.V.n <= quantile(Residuals.V.n, probs=c(1-Trim.n/100))
#points( Residuals.V.n[!t.b] ~ Temp_degK.V.n[!t.b], col="red")
# Trimmed non-linear regression
NLS_trim.L <- minpack.lm::nlsLM(formula=R_eco ~ fLloydTaylor(R_ref, E_0, Temp, T_ref.n=T_ref.n), algorithm='default', trace=FALSE,
data=as.data.frame(cbind(R_eco=NEEnight.V.n[t.b], Temp=Temp_degK.V.n[t.b]))
, start=coef(NLS.L)
#, start=list(R_ref=2,E_0=200)
)
##value<< Numeric vector with following components:
res <- c(
R_ref=coef(summary(NLS.L))['R_ref',1] ##<< Estimate of espiration rate at reference temperature
, R_ref_SD=coef(summary(NLS.L))['R_ref',2] ##<< Standard deviation of R_ref
,E_0=coef(summary(NLS.L))['E_0',1] ##<< Estimate of temperature sensitivity ("activation energy") in Kelvin (degK) for untrimmed dataset
, E_0_SD=coef(summary(NLS.L))['E_0',2] ##<< Standard deviation of E_0
,E_0_trim=coef(summary(NLS_trim.L))['E_0',1] ##<< Estimate of temperature sensitivity ("activation energy") in Kelvin (degK) for trimmed dataset
, E_0_trim_SD=coef(summary(NLS_trim.L))['E_0',2] ##<< Standard deviation of E_0_trim
##end<<
)
}, error = function(e) {
if( isTRUE(recoverOnError) ) recover()
res <- c(R_ref=NA, R_ref_SD=NA, E_0=NA, E_0_SD=NA, E_0_trim=NA, E_0_trim_SD=NA)
} ) #Spaces between brackets required to avoid replacement on documentation generation
res
}
fRegrE0fromShortTerm = function(
##title<<
## Estimation of the temperature sensitivity E_0
##description<<
## Estimation of the temperature sensitivity E_0 from regression of \code{\link{fLloydTaylor}} for short periods
NightFlux.V.n ##<< (Original) nighttime ecosystem carbon flux, i.e. respiration vector
,TempVar.V.n ##<< (Original) air or soil temperature vector (degC)
,DayCounter.V.i ##<< Integer specifying the day of each record
,WinDays.i=7 ##<< Window size for \code{\link{fLloydTaylor}} regression in days
,DayStep.i=5 ##<< Window step for \code{\link{fLloydTaylor}} regression in days
,TempRange.n=5 ##<< Threshold temperature range to start regression (#! Could be larger for Tair)
,Trim.n=5 ##<< Percentile to trim residual (%)
,NumE_0.n=3 ##<< Number of best E_0's to average over
,MinE_0.n=30 ##<< Minimum E0 for validity check
,MaxE_0.n=450 ##<< Maximum E0 for validity check
,T_ref.n #=273.15+15 ##<< Reference temperature in Kelvin (degK) used in \code{fLloydTaylor} for regressing Flux and Temperature
,CallFunction.s='' ##<< Name of function called from
,optimAlgorithm='default' ##<< optimization algorithm used (see \code{\link{nls}} ) or 'LM' for Levenberg-Marquard (see nlsLM from package minpack.lm)
)
##author<<
## AMM
{
##details<<
## The coefficient E0 is estimated for windows with a length of \code{WinDays.i} days,
## for successive periods in steps of \code{DayStep.i} days.
## Only those windows with a sufficient number or records and with a sufficient temperature range \code{TempRange.n}
## are used for the \code{\link{fLloydTaylor}} regression of E0 using the optimization \code{\link{fOptimSingleE0}}.
## Unreasonable estimates are discarded (95% confidence interval inside \code{MinE_0.n} and \code{MaxE_0.n}) and
## the others are ordered by their standard deviations.
## The mean across the best (=lowest standard deviation) E0 estimates is reported
## with \code{NumE_0.n} defining the number of best estimates to use.
# Regression settings
#NLSRes.F <- data.frame(NULL) #Results of non-linear regression
#NLSRes_trim.F <- data.frame(NULL) #Results of non-linear regression
MinData.n <- 6 # Minimum number of data points
fOptim <- fOptimSingleE0
if( optimAlgorithm=='LM') fOptim <- fOptimSingleE0_Lev
#tw: better use rbind with a list instead of costly repeated extending a data.frame
NLSRes.F <- as.data.frame(do.call( rbind, NLSRes.l <- lapply( seq(WinDays.i+1, max(DayCounter.V.i), DayStep.i) ,function(DayMiddle.i){
#TEST: DayMiddle.i <- 8
DayStart.i <- DayMiddle.i-WinDays.i
DayEnd.i <- DayMiddle.i+WinDays.i
#! Window size of 7 days corresponds to full window length of 15 days as in paper, non-congruent with PV-Wave code of 14 days
#! New code: Last window has minimum width of WinDays.i
Subset.b <- DayCounter.V.i >= DayStart.i & DayCounter.V.i <= DayEnd.i & !is.na(NightFlux.V.n) & !is.na(TempVar.V.n)
NEEnight.V.n <- subset(NightFlux.V.n, Subset.b)
Temp.V.n <- subset(TempVar.V.n, Subset.b)
Temp_degK.V.n <- fConvertCtoK(Temp.V.n)
if( length(NEEnight.V.n) > MinData.n && diff(range(Temp_degK.V.n)) >= TempRange.n ) {
#CountRegr.i <- CountRegr.i+1
resOptim <- fOptim( NEEnight.V.n, Temp_degK.V.n, algorithm=optimAlgorithm, T_ref.n=T_ref.n)
NLSRes.F <- c(Start=DayStart.i, End=DayEnd.i, Num=length(NEEnight.V.n), TRange=diff(range(Temp_degK.V.n)),
resOptim)
} else NULL
} ) ))
Limits.b <- ( NLSRes.F$E_0_trim - NLSRes.F$E_0_trim_SD > MinE_0.n & NLSRes.F$E_0_trim + NLSRes.F$E_0_trim_SD < MaxE_0.n )
#! New code: Check validity with SD (standard deviation) limits, in PV-Wave without SD, in paper if E_0_SD < (E_0 * 50%)
NLSRes.F$E_0_trim_ok <- ifelse( Limits.b, NLSRes.F$E_0_trim, NA)
NLSRes.F$E_0_trim_SD_ok <- ifelse( Limits.b, NLSRes.F$E_0_trim_SD, NA)
#
# Sort data frame for smallest standard deviation
NLSsort.F <- NLSRes.F[order(NLSRes.F$E_0_trim_SD_ok),] # ordered data.frame
##details<<
## Take average of the three E_0 with lowest standard deviation
E_0_trim.n <- round(mean(NLSsort.F$E_0_trim_ok[1:NumE_0.n]), digits=2)
#
# Abort flux partitioning if regression of E_0 failed
if( is.na(E_0_trim.n) ) {
# twutz 150226: just warning and returning negative value gives problems later on, maybe better stop and catch exception
warning(CallFunction.s, ':::fRegrE0fromShortTerm::: Less than ', NumE_0.n, ' valid values for E_0 after regressing ',
nrow(NLSRes.F), ' periods! Aborting partitioning.')
return(-111)
}
E_0_trim.n
##value<<
## Estimated scalar temperature sensitivity (E_0, degK)
}
sEddyProc$methods(
sRegrE0fromShortTerm = function(
##title<<
## sEddyProc$sRegrE0fromShortTerm - Estimation of the temperature sensitivity E_0
##description<<
## Estimation of the temperature sensitivity E_0 from regression of \code{\link{fLloydTaylor}}
## for short periods by calling \code{\link{fRegrE0fromShortTerm}}
NightFlux.s ##<< Variable with (original) nighttime ecosystem carbon flux, i.e. respiration
,TempVar.s ##<< Variable with (original) air or soil temperature (degC)
,... ##<< Parameters passed to \code{\link{fRegrE0fromShortTerm}}
,CallFunction.s='' ##<< Name of function called from
,debug.l = list(fixedE0=NA) ##<< List with controls for debugging, see details
)
##author<<
## AMM, TW
{
'Estimation of the temperature sensitivity E_0 from regression of fLloydTaylor() for short periods by calling fRegrE0fromShortTerm()'
##details<< For further details see \code{\link{fRegrE0fromShortTerm}}.
#
# Check if specified columns are numeric
SubCallFunc.s <- paste(CallFunction.s,'sRegrE0fromShortTerm', sep=':::')
ds <- cbind(sDATA,sTEMP)
fCheckColNames(ds, c(NightFlux.s, TempVar.s), SubCallFunc.s)
fCheckColNum(ds, c(NightFlux.s, TempVar.s), SubCallFunc.s)
#
##details<< \describe{ \item{Debugging control}{
## When supplying a finite scalar value \code{debug.l$fixedE0}, then this value
## is used instead of the temperature sensitivity E_0 from short term data.
## }}
if( (length(debug.l$fixedE0) != 0) && is.finite(debug.l$fixedE0) ){
E_0_trim.n <- debug.l$fixedE0[1]
message('Using prescribed temperature sensitivity E_0 of: ', format(E_0_trim.n, digits=5), '.')
}else{
# Write into vectors since cbind of data frames is so slow (factor of ~20 (!))
NightFlux.V.n <- ds[,NightFlux.s]
TempVar.V.n <- ds[,TempVar.s]
DayCounter.V.i <- c(1:sINFO$DIMS) %/% sINFO$DTS
#
# Check for validity of E_0 regression results
if( grepl('Tair', TempVar.s) ) {
#Limits in PV-Wave code for Tair
MinE_0.n <- 30; MaxE_0.n <- 350
} else if( grepl('Tsoil', TempVar.s) ) {
#Limits in PV-Wave code for Tsoil
MinE_0.n <- 30; MaxE_0.n <- 550 # Higher values due to potentially high Q10 values
} else {
#Default limits taken from paper
MinE_0.n <- 30; MaxE_0.n <- 450
}
#
E_0_trim.n <- fRegrE0fromShortTerm( NightFlux.V.n, TempVar.V.n, DayCounter.V.i, ...,
MinE_0.n=MinE_0.n, MaxE_0.n=MaxE_0.n, CallFunction.s=SubCallFunc.s)
message('Estimate of the temperature sensitivity E_0 from short term data: ', format(E_0_trim.n, digits=5), '.')
}
# Add constant value of E_0 as column vector to sTEMP
E_0.V.n <- rep(E_0_trim.n, nrow(sTEMP))
attr(E_0.V.n, 'varnames') <- 'E_0'
attr(E_0.V.n, 'units') <- 'degK'
E_0.V.n
##value<<
## Data vector with (constant) temperature sensitivity (E_0, degK)
## On failure, all entries are set to -111
})
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
sEddyProc$methods(
sRegrRref = function(
##title<<
## sEddyProc$sRegrRref - Estimation of the time series of reference respiration Rref
##description<<
## Estimation of the reference respiration Rref of \code{\link{fLloydTaylor}} for successive periods
NightFlux.s ##<< Variable with (original) nighttime ecosystem carbon flux, i.e. respiration
,TempVar.s ##<< Variable with (original) air or soil temperature (degC)
,E_0.s ##<< Temperature sensitivity E_0 estimated with \code{\link{sRegrE0fromShortTerm}}
,T_ref.n #=273.15+15 ##<< Reference temperature in Kelvin (degK) used in \code{fLloydTaylor} for regressing Flux and Temperature
,WinDays.i=3 ##<< Window size for \code{\link{fLloydTaylor}} regression in days
,DayStep.i=4 ##<< Window step for \code{\link{fLloydTaylor}} regression in days
,CallFunction.s='' ##<< Name of function called from
)
##author<<
## AMM
##details<<
## While parameter E0 in the Temperature-Respiration relationship (\code{\link{fLloydTaylor}}) is kept konstant,
## parameter Rref is allowed to change with time.
## This method estimates Rref for a series of time windows of length 2*\code{WinDays.i}+1 days
## shifted by \code{DayStep.i} days.
##
## For some of the windows, it maybe not be possible to estimate Rref. These missing values are filled by linear
## interpolation by function \code{\link{fInterpolateGaps}}.
{
'Estimation of the reference respiration Rref of fLloydTaylor() for successive periods'
# Check if specified columns are numeric
SubCallFunc.s <- paste(CallFunction.s,'sRegrRref', sep=':::')
ds <- cbind(sDATA,sTEMP)
fCheckColNames(ds, c(NightFlux.s, TempVar.s, E_0.s), SubCallFunc.s)
fCheckColNum(ds, c(NightFlux.s, TempVar.s, E_0.s), SubCallFunc.s)
# Regression settings
LMRes.F <- data.frame(NULL) #Results of linear regression
MinData.n <- 2 # Minimum number of data points for regression
# Write into vectors since cbind of data frames is so slow (factor of ~20 (!))
NightFlux.V.n <- ds[,NightFlux.s]
TempVar.V.n <- ds[,TempVar.s]
# Loop regression periods
DayCounter.V.i <- c(1:sINFO$DIMS) %/% sINFO$DTS
CountRegr.i <- 0
for (DayMiddle.i in seq(WinDays.i+1, max(DayCounter.V.i), DayStep.i)) {
#TEST: DayMiddle.i <- 8
DayStart.i <- DayMiddle.i-WinDays.i
DayEnd.i <- DayMiddle.i+WinDays.i
#! Window size of 4 days corresponds to a full window length of 9 days, non-congruent with PV-Wave code of 8 days, in paper not mentioned
#! New code: Last window has minimum of window size
Subset.b <- DayCounter.V.i >= DayStart.i & DayCounter.V.i <= DayEnd.i & !is.na(NightFlux.V.n) & !is.na(TempVar.V.n)
MeanHour.i <- round(mean(which(Subset.b))) # Weighted middle of the time period
NEEnight.V.n <- subset(NightFlux.V.n, Subset.b)
Temp.V.n <- subset(TempVar.V.n, Subset.b)
Temp_degK.V.n <- fConvertCtoK(Temp.V.n)
E_0.V.n <- subset(sTEMP[,E_0.s], Subset.b) # (Constant value)
if( length(NEEnight.V.n) > MinData.n ) {
CountRegr.i <- CountRegr.i+1
tryCatch({
LM.L <- lm(R_eco ~ 0 + fLloydTaylor(R_ref, E_0, Temp_degK, T_ref.n=T_ref.n), data=as.data.frame(cbind(R_eco=NEEnight.V.n, R_ref=1, E_0=E_0.V.n, Temp_degK=Temp_degK.V.n)))
LMRes.F <- rbind(LMRes.F, cbind(Start=DayStart.i, End=DayEnd.i, Num=length(NEEnight.V.n), MeanH=MeanHour.i,
R_ref=coef(summary(LM.L))[1], R_ref_SD=coef(summary(LM.L))[2]))
#! Note on PV-Wave code: trimmed linear regression not used in the end, i.e. in online webtool
if( F ) { # Plot for testing
plot(NEEnight.V.n ~ fLloydTaylor(1, E_0.V.n, Temp_degK.V.n, T_ref.n=T_ref.n))
curve(coef(LM.L)[1] * x, add=T, col='green')
}
}, error = function(e) {
LMRes.F <- rbind(LMRes.F, cbind(Start=DayStart.i, End=DayEnd.i, Num=length(NEEnight.V.n), MeanH=MeanHour.i,
R_ref=NA, R_ref_SD=NA))
} ) #Spaces between brackets required to avoid replacement on documentation generation
}
}
# Check for validity of regression results
LMRes.F$R_ref_ok <- ifelse( LMRes.F$R_ref < 0 , NA, LMRes.F$R_ref )
#! New code: Omit regressions with R_ref <0, in PV-Wave smaller values are set to 0.000001, not mentioned in paper
#TODO later: Flag for long distances between R_refs, especially if long distance in the beginning
#TODO later: Provide some kind of uncertainty estimate from R_ref_SD
# Interpolate R_ref periods linearly between MeanHour.i and keep constant at beginning/end
Rref.V.n <- rep(NA, length(NightFlux.V.n))
Rref.V.n[LMRes.F$MeanH] <- LMRes.F$R_ref_ok
Rref.V.n <- fInterpolateGaps(Rref.V.n)
attr(Rref.V.n, 'varnames') <- 'R_ref'
attr(Rref.V.n, 'units') <- attr(sTEMP[,NightFlux.s], 'units')
message('Regression of reference temperature R_ref for ', sum(!is.na(LMRes.F$R_ref_ok)), ' periods.')
Rref.V.n
##value<<
## Data vector (length number of windows) with reference respiration (R_ref, flux units)
})
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
sEddyProc$methods(
sMRFluxPartition = function(
##title<<
## sEddyProc$sMRFluxPartition - Flux partitioning after Reichstein et al. (2005)
##description<<
## Nighttime-based partitioning of measured net ecosystem fluxes into gross primary production (GPP) and ecosystem respiration (Reco)
FluxVar.s='NEE_f' ##<< Variable name of column with original and filled net ecosystem fluxes (NEE)
,QFFluxVar.s='NEE_fqc' ##<< Quality flag of NEE variable
,QFFluxValue.n=0 ##<< Value of quality flag for _good_ (original) data
,TempVar.s='Tair_f' ##<< Filled air- or soil temperature variable (degC)
,QFTempVar.s='Tair_fqc'##<< Quality flag of filled temperature variable
,QFTempValue.n=0 ##<< Value of temperature quality flag for _good_ (original) data
,RadVar.s='Rg' ##<< Unfilled (original) radiation variable
,Lat_deg.n ##<< Latitude in (decimal) degrees
,Long_deg.n ##<< Longitude in (decimal) degrees
,TimeZone_h.n ##<< Time zone (in hours)
,T_ref.n=273.15+15 ##<< Reference temperature in Kelvin (degK) used in \code{fLloydTaylor} for regressing Flux and Temperature
,Suffix.s = '' ##<< String suffix needed for different processing setups on the same dataset (for explanations see below)
,debug.l=list( ##<< List with debugging control (passed also to \code{\link{sRegrE0fromShortTerm}}).
##describe<<
useLocaltime.b=FALSE ##<< see details on solar vs local time
##end<<
)
,parsE0Regression=list() ##<< list with further parameters passed down to \code{\link{sRegrE0fromShortTerm}} and \code{\link{fRegrE0fromShortTerm}}, such as \code{TempRange.n}
)
##author<<
## AMM,TW
##references<<
## Reichstein M, Falge E, Baldocchi D et al. (2005) On the separation of net ecosystem exchange
## into assimilation and ecosystem respiration: review and improved algorithm. Global Change Biology, 11, 1424-1439.
{
'Partitioning of measured net ecosystem fluxes into gross primary production (GPP) and ecosystem respiration (Reco)'
##details<< \describe{\item{Description of newly generated variables with partitioning results:}{
## \itemize{
## \item PotRad - Potential radiation \cr
## \item FP_NEEnight - Good (original) NEE nighttime fluxes used for flux partitioning \cr
## \item FP_Temp - Good (original) temperature measurements used for flux partitioning \cr
## \item E_0 - Estimated temperature sensitivity \cr
## \item R_ref - Estimated reference respiration \cr
## \item Reco - Estimated ecosystem respiration \cr
## \item GPP_f - Estimated gross primary production \cr
## }
## }}
##details<< \describe{\item{Background}{
## This partitioning is based on the regression of nighttime respiration with temperature
## using the Lloyd-Taylor-Function \code{\link{fLloydTaylor}}.
## First the temperature sensitivity E_0 is estimated from short term data, see \code{\link{sRegrE0fromShortTerm}}.
## Next the reference temperature R_ref is estimated for successive periods throughout the whole dataset (see \code{\link{sRegrRref}}).
## These estimates are then used to calculate the respiration during daytime and nighttime and with this GPP.
## Attention: Gap filling of the net ecosystem fluxes (NEE) and temperature measurements (Tair or Tsoil) is required
## prior to the partitioning!
## }}
# If the default variable names are used and a Suffix.s is provided, then the suffix is added to the variable name (see also comment below)
if( FluxVar.s=='NEE_f' && QFFluxVar.s=='NEE_fqc' && fCheckValString(Suffix.s) ) {
FluxVar.s <- paste('NEE_', Suffix.s, '_f', sep='')
QFFluxVar.s <- paste('NEE_', Suffix.s, '_fqc', sep='')
}
ds <- cbind(sDATA,sTEMP)
# Check if specified columns exist in sDATA or sTEMP and if numeric and plausible. Then apply quality flag
fCheckColNames(ds, c(FluxVar.s, QFFluxVar.s, TempVar.s, QFTempVar.s, RadVar.s), 'sMRFluxPartition')
fCheckColNum(ds, c(FluxVar.s, QFFluxVar.s, TempVar.s, QFTempVar.s, RadVar.s), 'sMRFluxPartition')
fCheckColPlausibility(ds, c(FluxVar.s, QFFluxVar.s, TempVar.s, QFTempVar.s, RadVar.s), 'sMRFluxPartition')
Var.V.n <- fSetQF(ds, FluxVar.s, QFFluxVar.s, QFFluxValue.n, 'sMRFluxPartition')
message('Start flux partitioning for variable ', FluxVar.s, ' with temperature ', TempVar.s, '.')
##details<< \describe{\item{Selection of daytime data based on solar time}{
## The respiration-temperature regression is very
## sensitive to the selection of night- and daytime data.
## Nighttime is selected by a combined threshold of current solar radiation and potential radiation.
## The current implementation calculates potential radiation based on exact solar time, based on latitude and longitude.
## (see \code{\link{fCalcPotRadiation}} )
## Therefore it might differ from implementations that use local winter clock time instead.
## }}
# Calculate potential radiation
#! New code: Local time and equation of time accounted for in potential radiation calculation
DoY.V.n <- as.numeric(format(sDATA$sDateTime, '%j'))
Hour.V.n <- as.numeric(format(sDATA$sDateTime, '%H')) + as.numeric(format(sDATA$sDateTime, '%M'))/60
sTEMP$PotRad_NEW <<- fCalcPotRadiation(DoY.V.n, Hour.V.n, Lat_deg.n, Long_deg.n, TimeZone_h.n
, useSolartime.b=!isTRUE(debug.l$useLocaltime.b) )
# Filter night time values only
#! Note: Rg <= 10 congruent with MR PV-Wave, in paper Rg <= 20
# Should be unfilled (original) radiation variable, therefore dataframe set to sDATA only
sTEMP$FP_VARnight <<- ifelse(sDATA[,RadVar.s] > 10 | sTEMP$PotRad_NEW > 0, NA, Var.V.n)
attr(sTEMP$FP_VARnight, 'varnames') <<- paste(attr(Var.V.n, 'varnames'), '_night', sep='')
attr(sTEMP$FP_VARnight, 'units') <<- attr(Var.V.n, 'units')
#! New code: Slightly different subset than PV-Wave due to time zone correction (avoids timezone offset between Rg and PotRad)
# Apply quality flag for temperature
sTEMP$FP_Temp_NEW <<- fSetQF(ds, TempVar.s, QFTempVar.s, QFTempValue.n, 'sMRFluxPartition')
# Estimate E_0 and R_ref (results are saved in sTEMP)
# twutz1508: changed to do.call in order to allow passing further parameters when calling sMRFluxPartition
sTEMP$E_0_NEW <<- do.call( .self$sRegrE0fromShortTerm, c(
list('FP_VARnight', 'FP_Temp_NEW', T_ref.n=T_ref.n, CallFunction.s='sMRFluxPartition', debug.l=debug.l)
, parsE0Regression))
if( sum(sTEMP$E_0_NEW==-111) != 0 )
return(invisible(-111)) # Abort flux partitioning if regression of E_0 failed
# Reanalyse R_ref with E_0 fixed
sTEMP$R_ref_NEW <<- sRegrRref('FP_VARnight', 'FP_Temp_NEW', 'E_0_NEW', T_ref.n=T_ref.n, CallFunction.s='sMRFluxPartition')
# Calculate the ecosystem respiration Reco
sTEMP$Reco_NEW <<- fLloydTaylor(sTEMP$R_ref_NEW, sTEMP$E_0_NEW, fConvertCtoK(ds[,TempVar.s]), T_ref.n=T_ref.n)
attr(sTEMP$Reco_NEW, 'varnames') <<- 'Reco'
attr(sTEMP$Reco_NEW, 'units') <<- attr(Var.V.n, 'units')
# Calculate the gross primary production GPP_f
sTEMP$GPP_NEW_f <<- -ds[,FluxVar.s] + sTEMP$Reco_NEW
sTEMP$GPP_NEW_fqc <<- ds[,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)
attr(sTEMP$GPP_NEW_f, 'varnames') <<- 'GPP_f'
attr(sTEMP$GPP_NEW_f, 'units') <<- attr(Var.V.n, 'units')
##details<< \describe{\item{Different processing setups on the same dataset}{
## Attention: When processing the same site data set with different setups for the gap filling or flux partitioning
## (e.g. due to different ustar filters),
## a string suffix is needed! This suffix is added to the result column names to distinguish the results of the different setups.
## If a Suffix.s is provided and if the defaults for FluxVar.s and QFFluxVar.s are used, the Suffix.s will be added to their variable names
## (e.g. 'NEE_f' will be renamed to 'NEE_WithUstar_f' and 'NEE_fqc' to 'NEE_WithUstar_fqc' for the Suffix.s='WithUstar').
## Currently, this works only with defaults of FluxVar.s='NEE_f' and QFFluxVar.s='NEE_fqc'.
## }}
# Rename new columns generated during flux partitioning:
# For nighttime NEE (FP_NEEnight or FP_NEEnight_Suffix)
colnames(sTEMP) <<- gsub('_VARnight', paste('_NEEnight', (if(fCheckValString(Suffix.s)) '_' else ''), Suffix.s, sep=''), colnames(sTEMP))
# For the results columns, the _NEW is dropped and the suffix added
colnames(sTEMP) <<- gsub('_NEW', paste((if(fCheckValString(Suffix.s)) '_' else ''), Suffix.s, sep=''), colnames(sTEMP))
# Check for duplicate columns (to detect if different processing setups were executed without different suffixes)
if( length(names(iDupl <- which(table(colnames(sTEMP)) > 1))) ) {
warning(paste0('sMRFluxPartition::: Duplicated columns found! (',
paste(names(iDupl),collapse=",")
,') Deleting each first of duplicate columns.'
,' Please use different Suffix.s when processing different setups on the same dataset!'))
for( cname in names(iDupl) ) sTEMP[cname] <<- NULL # need to remove columns else some tests fail
}
return(invisible(NULL))
##value<<
## Flux partitioning results (see variables in details) in sTEMP data frame (with renamed columns).
## On success, return value is NULL. On failure an integer scalar error code is returned:
## -111 if regression of E_0 failed due to insufficient relationship in the data.
})
Try the REddyProc package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.