LakeMetabolizer meta.bookkeep() example

LM_mbkkp_test.R

#Test case for calculating metabolism
library(LakeMetabolizer)

#set up test example where DO starts and ends at same level, is at saturation (= no flux)
DOobs = c(3, 2, 3, 4, 3, 2)
DOsat = c(3, 2, 3, 4, 3, 2)
kGas = c(.2, .2, .2, .2, .2, .2)
zMix = c(1, 1, 1, 1, 1, 1)
irr = c(0, 1, 1, 0, 0, 0)

metab.bookkeep(do.obs = DOobs, do.sat = DOsat, k.gas = kGas, z.mix = zMix, irr = irr)

do.obs = DOobs 
do.sat = DOsat
k.gas = kGas
z.mix = zMix
irr = irr

plot(DOobs, type="l", main="'Midnight to midnight'")
polygon(x = c(1, 2, 2, 1), y = c(2,2,4,4), col="light grey")
polygon(x = c(2, 4, 4, 2), y = c(2,2,4,4), col="light yellow")
polygon(x = c(4, 6, 6, 4), y = c(2,2,4,4), col="light grey")
lines(DOobs, type="l")


#Below is stripped down code from metab.bookkeep with fixes on line 29 and 42
nobs <- length(do.obs)
# freq <- nobs
freq <- nobs - 1 # SUBTRACT ONE B/C THAT'S THE NUMBER OF DELTA DO'S
dayI <- irr == 1L
nightI <- irr == 0L
delta.do <- diff(do.obs)
miss.delta <- sum(is.na(delta.do))
gas.flux <- (do.sat - do.obs) * (k.gas/freq)/z.mix
delta.do.metab <- delta.do - gas.flux[1:(length(gas.flux) - 1)]
nep.day <- delta.do.metab[dayI]
nep.night <- delta.do.metab[nightI]
R <- mean(nep.night, na.rm = TRUE) * freq
NEP <- mean(delta.do.metab, na.rm = TRUE) * freq
# LakeMetabolizer formula: 
# GPP <- mean(nep.day, na.rm = TRUE) * sum(dayI) - R
## Proposed fix
GPP <- mean(nep.day, na.rm = TRUE) * sum(dayI) - mean(nep.night, na.rm = TRUE) * sum(dayI)
metab <- data.frame(GPP = GPP, R = R, NEP = NEP)

print("Metabolism with changed GPP and freq formulas above:")
print(metab)

print("Metabolism from LakeMetabolizer")
metab.bookkeep(do.obs = DOobs, do.sat = DOsat, k.gas = kGas, z.mix = zMix, irr = irr)