metab {LakeMetabolizer}

Calculate metabolism

Description

Returns daily time series of gross primary production (GPP), respiration (R), and net ecosystem production (NEP). Depending on the method used, other information may be returned as well. Calculations are made using one of 5 statistical methods.

Usage

metab(data, method, wtr.name="wtr", irr.name="irr", do.obs.name="do.obs", ...)

Arguments

data

a data.frame whose columns are

"datetime" = class POSIXct vector

"do.obs" = numeric vector of oxygen concentration in mg/L

"do.sat" = numeric vector of saturated oxygen concentration in mg/L

"k.gas" = numeric vector of gas exchange coefficient values in m/day, should be 0 when depth of do.obs is deeper than z.mix

"z.mix" = numeric vector of mixing depth values in meters

"irr" = numeric vector of PAR values, arbitrary units

"wtr" = numeric vector of water temperature values, arbitrary units

Columns that are not used by a particular statistical method do not need to be supplied.

method

a character string specifying one of the 5 statistical methods (bayesian, bookkeep, kalman, ols, mle)

wtr.name

the name of the column containing temperature at the depth of do.obs (predictor variable for R)

irr.name

the name of the column containing irradiance (predictor variable for GPP)

do.obs.name

the name of the column in data containing the DO observations (in mg/L) to be used as the response variable

...

arguments to be passed on to the metabolism model specified by method

Value

A data.frame containing columns for year, doy (day of year, julian day plus fraction of day), GPP, R, and NEP

year

integer year

doy

numeric, day of year + fraction of day, where the day is the julian day, and a fraction of 0.5 corresponds to noon

GPP

numeric, gross primary production, in units of mg O2 per liter per day. By convention, this value is positive.

R

numeric, respiration, in units of mg O2 per liter per day. By convention, this value is negative

NEP

numeric, net ecosystem production, in units of mg O2 per liter per day. For most methods this equal GPP+R, but this is not necessarily the case for "method"="bookkeep"

Note that different models will have different attributes attached to them. See examples.

Author(s)

Ryan D. Batt

See Also

Metabolism models: metab.bookkeep, metab.ols, metab.mle, metab.kalman, metab.bayesian

For smoothing noisy temperature: temp.kalman

To calculate do.sat: o2.at.sat

To calculate k.gas: k600.2.kGAS

To calculate k600 values for k.gas: k.cole, k.crusius, k.macIntyre, k.read

Examples

# fake data
datetime <- seq(as.POSIXct("2014-06-16 00:00:00", tz="GMT"),
     as.POSIXct("2014-06-17 23:55:00", tz="GMT"), length.out=288*2)
do.obs <- 2*sin(2*pi*(1/288)*(1:(288*2))+1.1*pi) + 8 + rnorm(288*2, 0, 0.5)
wtr <- 3*sin(2*pi*(1/288)*(1:(288*2))+pi) + 17 + rnorm(288*2, 0, 0.15)
do.sat <- LakeMetabolizer::o2.at.sat.base(wtr, 960)
irr <- (1500*sin(2*pi*(1/288)*(1:(288*2))+1.5*pi) +650 + rnorm(288*2, 0, 0.25)) *
     ifelse(is.day(datetime, 42.3), 1, 0)
k.gas <- 0.4
z.mix <- 1

# plot time series
plot(wtr, type="l", xaxt="n", yaxt="n", xlab="", ylab="")
par(new=TRUE); plot(do.obs, type="l", col="blue", xaxt="n", yaxt="n", xlab="", ylab="")
par(new=TRUE); plot(irr, type="l", col="orange", xaxt="n", yaxt="n", xlab="", ylab="")
abline(v=144, lty="dotted")
abline(v=288)
legend("topleft", legend=c("wtr", "do.obs", "irr"), lty=1,
     col=c("black", "blue", "orange"), inset=c(0.08, 0.01))

# put data in a data.frame
data <- data.frame(datetime=datetime, do.obs=do.obs, do.sat=do.sat, k.gas=k.gas,
     z.mix=z.mix, irr=irr, wtr=wtr)

# run each metabolism model
m.bk <- metab(data, "bookkeep", lake.lat=42.6)
 m.bk <- metab(data, lake.lat=42.6) # no method defaults to "bookeep"
m.ols <- metab(data, "ols", lake.lat=42.6)
m.mle <- metab(data, "mle", lake.lat=42.6)
m.kal <- metab(data, "kalman", lake.lat=42.6)
## Not run: m.bay <- metab(data, "bayesian", lake.lat=42.6)

# example attributes
names(attributes(m.ols))
attr(m.ols, "mod")

# To get full JAGS model
# including posterior draws:
## Not run: names(attributes(m.bay))
## Not run: attr(m.bay, "model")