From efaa75a48aeeff92df24b05b646b3377c0b69ffb Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 1 Apr 2026 10:19:35 -0400 Subject: [PATCH 01/20] Constrain c1 and c2 positive --- src/mleLoop.c | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/mleLoop.c b/src/mleLoop.c index 9c93d20..6784904 100644 --- a/src/mleLoop.c +++ b/src/mleLoop.c @@ -6,17 +6,21 @@ void mleLoopCoe(double *alpha, double *doobs, double *c1, double *c2, double *be double a1, z=0; // z has to be a double b/c kz is passed in as a double (for cases where it isn't 0) // double alpha=*alpha, doobs=*doobs, c1=*c1, c2=*c2, beta=*beta, irr=*irr, wtr=*wtr, kz=*kz, dosat=*dosat; int ni=*nobs; + double c1_pos = exp(*c1); + double c2_pos = exp(*c2); for(i=1; i Date: Wed, 1 Apr 2026 10:21:32 -0400 Subject: [PATCH 02/20] Correct C calls --- R/metab.mle.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/metab.mle.R b/R/metab.mle.R index 4df7de1..8bb4c8d 100644 --- a/R/metab.mle.R +++ b/R/metab.mle.R @@ -167,7 +167,7 @@ metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", . # ============================================ mleLoopOE <- function(alpha, doobs, c1, c2, beta, irr, wtr, kz, dosat){ nobs <- length(doobs) - a.loop <- .C("mleLoopCoe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer") + a.loop <- .C("mleLoopCoe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer2") return(a.loop[["alpha"]]) } @@ -176,7 +176,7 @@ mleLoopOE <- function(alpha, doobs, c1, c2, beta, irr, wtr, kz, dosat){ # ============================================ mleLoopPE <- function(alpha, doobs, c1, c2, beta, irr, wtr, kz, dosat){ nobs <- length(doobs) - a.loop <- .C("mleLoopCpe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer") + a.loop <- .C("mleLoopCpe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer2") return(a.loop[["alpha"]]) } From a6a34ca5bc4dd54d47b59682e141fc0b17b27586 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 1 Apr 2026 10:35:24 -0400 Subject: [PATCH 03/20] Undo C call updates --- R/metab.mle.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/metab.mle.R b/R/metab.mle.R index 8bb4c8d..4df7de1 100644 --- a/R/metab.mle.R +++ b/R/metab.mle.R @@ -167,7 +167,7 @@ metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", . # ============================================ mleLoopOE <- function(alpha, doobs, c1, c2, beta, irr, wtr, kz, dosat){ nobs <- length(doobs) - a.loop <- .C("mleLoopCoe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer2") + a.loop <- .C("mleLoopCoe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer") return(a.loop[["alpha"]]) } @@ -176,7 +176,7 @@ mleLoopOE <- function(alpha, doobs, c1, c2, beta, irr, wtr, kz, dosat){ # ============================================ mleLoopPE <- function(alpha, doobs, c1, c2, beta, irr, wtr, kz, dosat){ nobs <- length(doobs) - a.loop <- .C("mleLoopCpe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer2") + a.loop <- .C("mleLoopCpe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer") return(a.loop[["alpha"]]) } From e2d219a35437f5a0bc2d88ead33fdf7f028a9f23 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Mon, 8 Jun 2026 15:20:51 -0400 Subject: [PATCH 04/20] Added k.klaus and k.klaus.base --- NAMESPACE | 2 ++ R/k.klaus.R | 62 ++++++++++++++++++++++++++++++++++++++++++++++ R/k.read.R | 31 ++++++++++++++++++++--- man/k.read.Rd | 18 ++++++++++++-- man/k.read.base.Rd | 18 ++++++++++++-- 5 files changed, 123 insertions(+), 8 deletions(-) create mode 100644 R/k.klaus.R diff --git a/NAMESPACE b/NAMESPACE index 0298434..9c78e66 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -15,6 +15,8 @@ export(k.crusius) export(k.crusius.base) export(k.heiskanen) export(k.heiskanen.base) +export(k.klaus) +export(k.klaus.base) export(k.macIntyre) export(k.macIntyre.base) export(k.read) diff --git a/R/k.klaus.R b/R/k.klaus.R new file mode 100644 index 0000000..8571f50 --- /dev/null +++ b/R/k.klaus.R @@ -0,0 +1,62 @@ +# ---Author: Bennett McAfee, 2026-06-08 --- +# Last update: 2026-06-08 + +#'@export +k.klaus = function(ts.data, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", "power")){ + + if(!has.vars(ts.data, 'wnd')){ + stop('k.klaus requires a "wnd" (wind speed) column in the supplied data') + } + + wind = get.vars(ts.data, 'wnd') + + k600 = k.klaus.base(wind[,2], wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", "power")) + + return(data.frame(datetime=ts.data$datetime, k600=k600)) +} + + +#'@export +k.klaus.base <- function(wnd, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", "power")) { + + method <- match.arg(method) + + # converting m2 to km2 + lake.area <- lake.area / 1e6 + + # Converting uz to u10 + if (wnd.z != 10){ + wnd <- LakeMetabolizer::wind.scale.base(wnd = wnd, wnd.z = wnd.z) + } + + # helper logit function + logit_custom <- function(x){ + return(log(x / (1 - x))) + } + + # sanity checks + if (any(sin <= 0 | sin >= 1)) { + stop("sin must be between 0 and 1 (exclusive)") + } + + if (method == "linear") { + k600 <- (0.328 * log10(lake.area) + 1.581) * wnd - 0.066 * logit_custom(sin) + 1.266 + } else if (method == "power") { + k600 <- (0.281 * log10(lake.area) + 1.361) * (wnd ^ 1.097) - 0.072 * logit_custom(sin) + 1.401 + } else if (method == "exp") { + if (is.null(sdi)) { + stop("sdi must be provided for exponential model") + } + k600 <- (-0.057 * logit_custom(sin) + 2.366) * exp(wnd * (0.144 * log10(sdi) + 0.156)) + } + + k600 <- k600 * (24/100) # convert cm/hr to m/day + + return(k600) +} + + +# -- References +# Klaus, Marcus, and Dominic Vachon. 2020. +# Challenges of Predicting Gas Transfer Velocity from Wind Measurements over Global Lakes. +# Aquatic Sciences 82 (3): 53. https://doi.org/10.1007/s00027-020-00729-9. diff --git a/R/k.read.R b/R/k.read.R index 49a1cb1..589c9ac 100644 --- a/R/k.read.R +++ b/R/k.read.R @@ -7,6 +7,7 @@ #'k.crusius #'k.vachon #'k.heiskanen +#'k.klaus #'@title Returns a timeseries of gas exchange velocity #'@description #'Returns the gas exchange velocity based on the chosen model in units of m/day @@ -24,13 +25,17 @@ #'k.vachon(ts.data, lake.area, params=c(2.51,1.48,0.39)) #' #'k.heiskanen(ts.data, wnd.z, Kd, atm.press) +#' +#'k.klaus(ts.data, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", "power")) #'@param ts.data vector of datetime in POSIXct format -#'@param method Only for \link{k.crusius}. String of valid method . Either "linear", "bilinear", or "power" +#'@param method Only for \link{k.crusius} and \link{k.klaus}. String of valid method . For k.crusius either "linear", "bilinear", or "power". For k.klaus either "linear", "exp", or "power" #'@param wnd.z height of wind measurement, m #'@param Kd Light attenuation coefficient (Units:m^-1) #'@param atm.press atmospheric pressure in mb #'@param lat Latitude, degrees north #'@param lake.area Lake area, m^2 +#'@param sdi Numeric value of shoreline development index. Only for \link{k.klaus}'s "exp" method. +#'@param sin Numeric value >0 and ≤1 of the scale of spatial integration. Only for \link{k.klaus}. #'@param params Only for \link{k.vachon.base} and \link{k.macIntyre}. See details. #' #'@details Can change default parameters of MacIntyre and Vachon models. Default for Vachon is @@ -70,14 +75,18 @@ #'\emph{An approach to estimation of near-surface turbulence and CO2 transfer #'velocity from remote sensing data}. Journal of Marine Systems 66, (2007): 182-194. #' +#'Marcus Klaus and Dominic Vachon. \emph{Challenges of Predicting Gas Transfer Velocity from Wind Measurements over Global Lakes}. +#'Aquatic Sciences 82 (3): 53. (2020). +#' #'@author -#'Hilary Dugan, Jake Zwart, Luke Winslow, R. Iestyn. Woolway, Jordan S. Read +#'Hilary Dugan, Jake Zwart, Luke Winslow, R. Iestyn. Woolway, Jordan S. Read, Bennett McAfee #'@seealso #'\link{k.cole} #'\link{k.crusius} #'\link{k.macIntyre} #'\link{k.vachon} #'\link{k.heiskanen} +#'\link{k.klaus} #'@examples #'data.path = system.file('extdata', package="LakeMetabolizer") #' @@ -101,6 +110,7 @@ #'atm.press = 1018 #'lat = tb.data$metadata$latitude #'lake.area = tb.data$metadata$lakearea +#'sin = 1/lake.area #' #'#for k.read and k.macIntyre, we need LW_net. #'#Calculate from the observations we have available. @@ -116,6 +126,8 @@ #' #'k600_macIntyre = k.macIntyre(ts.data, wnd.z=wnd.z, Kd=kd, atm.press=atm.press) #' +#'k600_klaus = k.klaus(ts.data, wnd.z=wnd.z, lake.area=lake.area, sin=sin) +#' #'@export k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ @@ -178,6 +190,7 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #'k.crusius.base #'k.vachon.base #'k.heiskanen.base +#'k.klaus.base #'@title Returns a timeseries of gas exchange velocity #'@description #'Returns the gas exchange velocity based on the chosen model in units of m/day @@ -199,8 +212,9 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #' #'k.heiskanen.base(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet) #' +#'k.klaus.base(wnd, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", "power")) #'@param wnd Numeric value of wind speed, (Units:m/s) -#'@param method Only for \link{k.crusius.base}. String of valid method . Either "constant", "bilinear", or "power" +#'@param method Only for \link{k.crusius.base} and \link{k.klaus.base}. String of valid method . For k.crusius.base either "constant", "bilinear", or "power". For k.klaus.base either "linear", "exp", or "power" #'@param wnd.z Height of wind measurement, (Units: m) #'@param Kd Light attenuation coefficient (Units: m^-1) #'@param lat Latitude, degrees north @@ -213,6 +227,8 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #'@param RH Numeric value of relative humidity, \% #'@param sw Numeric value of short wave radiation, W m^-2 #'@param lwnet Numeric value net long wave radiation, W m^-2 +#'@param sdi Numeric value of shoreline development index. Only for \link{k.klaus.base}'s "exp" method. +#'@param sin Numeric value >0 and ≤1 of the scale of spatial integration. Only for \link{k.klaus.base}. #'@param params Optional parameter input, only for \link{k.vachon.base} and \link{k.macIntyre.base}. See details. #'@details Can change default parameters of MacIntyre and Vachon models. Default for Vachon is #'c(2.51,1.48,0.39). Default for MacIntyre is c(1.2,0.4872,1.4784). Heiskanen et al. (2014) uses MacIntyre @@ -248,8 +264,11 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #'\emph{An approach to estimation of near-surface turbulence and CO2 transfer #'velocity from remote sensing data}. Journal of Marine Systems 66, (2007): 182-194. #' +#'Marcus Klaus and Dominic Vachon. \emph{Challenges of Predicting Gas Transfer Velocity from Wind Measurements over Global Lakes}. +#'Aquatic Sciences 82 (3): 53. (2020). +#' #'@author -#'R. Iestyn. Woolway, Hilary Dugan, Luke Winslow, Jordan S Read, GLEON fellows +#'R. Iestyn. Woolway, Hilary Dugan, Luke Winslow, Jordan S Read, Bennett McAfee, GLEON fellows #'@seealso #'\link{k.cole} #'\link{k.read} @@ -257,6 +276,7 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #'\link{k.macIntyre} #'\link{k.vachon} #'\link{k.heiskanen} +#'\link{k.klaus} #'@examples #'wnd.z <- 2 #'Kd <- 2 @@ -271,6 +291,7 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #'RH <- 90 #'sw <- 800 #'lwnet <- -55 +#'sin <- 1/lake.area #'timeStep <- 30 #' #'U10 <- wind.scale.base(wnd, wnd.z) @@ -288,6 +309,8 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #'k600_macInytre <- k.macIntyre.base(wnd.z, Kd, atm.press, #'dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet) #' +#'k600_klaus <- k.klaus.base(wnd, wnd.z, lake.area, sin) +#' #'@export k.read.base <- function(wnd.z, Kd, lat, lake.area, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet){ diff --git a/man/k.read.Rd b/man/k.read.Rd index 5ec6ba1..e34ad21 100644 --- a/man/k.read.Rd +++ b/man/k.read.Rd @@ -8,6 +8,7 @@ \alias{k.crusius} \alias{k.vachon} \alias{k.heiskanen} +\alias{k.klaus} \title{Returns a timeseries of gas exchange velocity} \usage{ k.cole(ts.data) @@ -23,6 +24,8 @@ k.macIntyre(ts.data, wnd.z, Kd, atm.press,params=c(1.2,0.4872,1.4784)) k.vachon(ts.data, lake.area, params=c(2.51,1.48,0.39)) k.heiskanen(ts.data, wnd.z, Kd, atm.press) + +k.klaus(ts.data, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", "power")) } \arguments{ \item{ts.data}{vector of datetime in POSIXct format} @@ -37,7 +40,11 @@ k.heiskanen(ts.data, wnd.z, Kd, atm.press) \item{lake.area}{Lake area, m^2} -\item{method}{Only for \link{k.crusius}. String of valid method . Either "linear", "bilinear", or "power"} +\item{method}{Only for \link{k.crusius} and \link{k.klaus}. String of valid method . For k.crusius either "linear", "bilinear", or "power". For k.klaus either "linear", "exp", or "power"} + +\item{sdi}{Numeric value of shoreline development index. Only for \link{k.klaus}'s "exp" method.} + +\item{sin}{Numeric value >0 and ≤1 of the scale of spatial integration. Only for \link{k.klaus}.} \item{params}{Only for \link{k.vachon.base} and \link{k.macIntyre}. See details.} } @@ -75,6 +82,7 @@ wnd.z = 10 #because we converted to u10 atm.press = 1018 lat = tb.data$metadata$latitude lake.area = tb.data$metadata$lakearea +sin = 1/lake.area #for k.read and k.macIntyre, we need LW_net. #Calculate from the observations we have available. @@ -90,6 +98,8 @@ atm.press=atm.press, lat=lat, lake.area=lake.area) k600_macIntyre = k.macIntyre(ts.data, wnd.z=wnd.z, Kd=kd, atm.press=atm.press) +k600_klaus = k.klaus(ts.data, wnd.z=wnd.z, lake.area=lake.area, sin=sin) + } \references{ Cole, J., J. Nina, and F. Caraco. \emph{Atmospheric exchange of carbon dioxide @@ -118,6 +128,9 @@ transfer coefficients in a boreal lake}. Tellus B 66, no.22827 (2014) Alexander Soloviev, Mark Donelan, Hans Graber, Brian Haus, Peter Schlussel. \emph{An approach to estimation of near-surface turbulence and CO2 transfer velocity from remote sensing data}. Journal of Marine Systems 66, (2007): 182-194. + +Marcus Klaus and Dominic Vachon. \emph{Challenges of Predicting Gas Transfer Velocity from Wind Measurements over Global Lakes}. +Aquatic Sciences 82 (3): 53. (2020). } \seealso{ \link{k.cole} @@ -125,9 +138,10 @@ velocity from remote sensing data}. Journal of Marine Systems 66, (2007): 182-19 \link{k.macIntyre} \link{k.vachon} \link{k.heiskanen} +\link{k.klaus} } \author{ -Hilary Dugan, Jake Zwart, Luke Winslow, R. Iestyn. Woolway, Jordan S. Read +Hilary Dugan, Jake Zwart, Luke Winslow, R. Iestyn. Woolway, Jordan S. Read, Bennett McAfee } \keyword{math} \keyword{methods} diff --git a/man/k.read.base.Rd b/man/k.read.base.Rd index d0ec805..7a0cfa2 100644 --- a/man/k.read.base.Rd +++ b/man/k.read.base.Rd @@ -8,6 +8,7 @@ \alias{k.crusius.base} \alias{k.vachon.base} \alias{k.heiskanen.base} +\alias{k.klaus.base} \title{Returns a timeseries of gas exchange velocity} \usage{ k.cole.base(wnd) @@ -26,6 +27,8 @@ lwnet, params=c(1.2,0.4872,1.4784)) k.vachon.base(wnd, lake.area, params=c(2.51,1.48,0.39)) k.heiskanen.base(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet) + +k.klaus.base(wnd, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", "power")) } \arguments{ \item{wnd.z}{Height of wind measurement, (Units: m)} @@ -54,7 +57,11 @@ k.heiskanen.base(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, l \item{lwnet}{Numeric value net long wave radiation, W m^-2} -\item{method}{Only for \link{k.crusius.base}. String of valid method . Either "constant", "bilinear", or "power"} +\item{method}{Only for \link{k.crusius.base} and \link{k.klaus.base}. String of valid method . For k.crusius.base either "constant", "bilinear", or "power". For k.klaus.base either "linear", "exp", or "power"} + +\item{sdi}{Numeric value of shoreline development index. Only for \link{k.klaus.base}'s "exp" method.} + +\item{sin}{Numeric value >0 and ≤1 of the scale of spatial integration. Only for \link{k.klaus.base}.} \item{params}{Optional parameter input, only for \link{k.vachon.base} and \link{k.macIntyre.base}. See details.} } @@ -84,6 +91,7 @@ wnd <- 6 RH <- 90 sw <- 800 lwnet <- -55 +sin <- 1/lake.area timeStep <- 30 U10 <- wind.scale.base(wnd, wnd.z) @@ -101,6 +109,8 @@ atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet) k600_macInytre <- k.macIntyre.base(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet) +k600_klaus <- k.klaus.base(wnd, wnd.z, lake.area, sin) + } \references{ Cole, J., J. Nina, and F. Caraco. \emph{Atmospheric exchange of carbon dioxide @@ -129,6 +139,9 @@ transfer coefficients in a boreal lake}. Tellus B 66, no.22827 (2014) Alexander Soloviev, Mark Donelan, Hans Graber, Brian Haus, Peter Schlussel. \emph{An approach to estimation of near-surface turbulence and CO2 transfer velocity from remote sensing data}. Journal of Marine Systems 66, (2007): 182-194. + +Marcus Klaus and Dominic Vachon. \emph{Challenges of Predicting Gas Transfer Velocity from Wind Measurements over Global Lakes}. +Aquatic Sciences 82 (3): 53. (2020). } \seealso{ \link{k.cole} @@ -137,9 +150,10 @@ velocity from remote sensing data}. Journal of Marine Systems 66, (2007): 182-19 \link{k.macIntyre} \link{k.vachon} \link{k.heiskanen} +\link{k.klaus} } \author{ -R. Iestyn. Woolway, Hilary Dugan, Luke Winslow, Jordan S Read, GLEON fellows +R. Iestyn. Woolway, Hilary Dugan, Luke Winslow, Jordan S Read, Bennett McAfee, GLEON fellows } \keyword{math} \keyword{methods} From 62ef8c51e12141ac3b6fe1a97630c4bd0b61074c Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Tue, 9 Jun 2026 09:23:25 -0400 Subject: [PATCH 05/20] k.klaus documentation quick fix --- R/k.read.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/k.read.R b/R/k.read.R index 589c9ac..af33e70 100644 --- a/R/k.read.R +++ b/R/k.read.R @@ -35,7 +35,7 @@ #'@param lat Latitude, degrees north #'@param lake.area Lake area, m^2 #'@param sdi Numeric value of shoreline development index. Only for \link{k.klaus}'s "exp" method. -#'@param sin Numeric value >0 and ≤1 of the scale of spatial integration. Only for \link{k.klaus}. +#'@param sin Numeric value >0 and <1 of the scale of spatial integration. Only for \link{k.klaus}. #'@param params Only for \link{k.vachon.base} and \link{k.macIntyre}. See details. #' #'@details Can change default parameters of MacIntyre and Vachon models. Default for Vachon is @@ -228,7 +228,7 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #'@param sw Numeric value of short wave radiation, W m^-2 #'@param lwnet Numeric value net long wave radiation, W m^-2 #'@param sdi Numeric value of shoreline development index. Only for \link{k.klaus.base}'s "exp" method. -#'@param sin Numeric value >0 and ≤1 of the scale of spatial integration. Only for \link{k.klaus.base}. +#'@param sin Numeric value >0 and <1 of the scale of spatial integration. Only for \link{k.klaus.base}. #'@param params Optional parameter input, only for \link{k.vachon.base} and \link{k.macIntyre.base}. See details. #'@details Can change default parameters of MacIntyre and Vachon models. Default for Vachon is #'c(2.51,1.48,0.39). Default for MacIntyre is c(1.2,0.4872,1.4784). Heiskanen et al. (2014) uses MacIntyre From aa7ede025d0c5a98e613d48ed7965bde90536b46 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Tue, 9 Jun 2026 10:51:47 -0400 Subject: [PATCH 06/20] Revert all MLE changes --- src/mleLoop.c | 17 ++++------------- 1 file changed, 4 insertions(+), 13 deletions(-) diff --git a/src/mleLoop.c b/src/mleLoop.c index 6784904..9c93d20 100644 --- a/src/mleLoop.c +++ b/src/mleLoop.c @@ -6,21 +6,17 @@ void mleLoopCoe(double *alpha, double *doobs, double *c1, double *c2, double *be double a1, z=0; // z has to be a double b/c kz is passed in as a double (for cases where it isn't 0) // double alpha=*alpha, doobs=*doobs, c1=*c1, c2=*c2, beta=*beta, irr=*irr, wtr=*wtr, kz=*kz, dosat=*dosat; int ni=*nobs; - double c1_pos = exp(*c1); - double c2_pos = exp(*c2); for(i=1; i Date: Tue, 9 Jun 2026 14:10:54 -0400 Subject: [PATCH 07/20] Constrain MLE estimates to correct sign --- R/metab.mle.R | 76 +++++++++++++++++++++++++-------------------------- 1 file changed, 38 insertions(+), 38 deletions(-) diff --git a/R/metab.mle.R b/R/metab.mle.R index 4df7de1..24779b6 100644 --- a/R/metab.mle.R +++ b/R/metab.mle.R @@ -130,7 +130,7 @@ metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", . Q0 <- ((diff(range(do.obs,na.rm=TRUE)) - mean(do.obs,na.rm=TRUE))^2 / length(do.obs)) - guesses <- c(1E-4, 1E-4, log(Q0)) + guesses <- c(log(1E-4), log(1E-4), log(Q0)) #We have a different number of fitted parameters depending on error type of the model if(error.type=='OE'){ @@ -139,7 +139,7 @@ metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", . fit <- optim(guesses, fn=mleNllOE, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) pars0 <- fit$par - pars <- c("gppCoeff"=pars0[1], "rCoeff"=pars0[2], "Q"=exp(pars0[3]), "nll"=fit$value, "doInit"=pars0[4]) + pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "nll"=fit$value, "doInit"=pars0[4]) }else if(error.type=='PE'){ guesses <- c(guesses) @@ -147,7 +147,7 @@ metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", . fit <- optim(guesses, fn=mleNllPE, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) pars0 <- fit$par - pars <- c("gppCoeff"=pars0[1], "rCoeff"=pars0[2], "Q"=exp(pars0[3]), "nll"=fit$value) + pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "nll"=fit$value) }else{ stop("error.type must be either 'OE' or 'PE', Observation Error or Process Error respectively.") @@ -184,45 +184,45 @@ mleLoopPE <- function(alpha, doobs, c1, c2, beta, irr, wtr, kz, dosat){ # = mle NLL function = # ==================== mleNllPE <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr){ - c1 <- Params[1] #PAR coeff - c2 <- Params[2] #log(Temp) coeff - Q <- exp(Params[3]) # Variance of the process error - - # See KalmanDO_smooth.R comments for explanation of beta - kz <- k.gas/z.mix # K and Zmix are both vector of length nobs - beta <- exp(-kz) # This beta is for using the differential equation form - - # Set first true value equal to first observation - alpha <- rep(0, length(do.obs)) - alpha[1] <- do.obs[1]#Let's give this model some starting values - - #R version of C loop - #for(i in 2:length(do.obs)){ - # a1 <- c1*irr[i-1] + c2*log(wtr[i-1]) + kz[i-1]*do.sat[i-1] - # alpha[i] <- a1/kz[i-1] + -exp(-kz[i-1])*a1/kz[i-1] + beta[i-1]*alpha[i-1] # NOTE: beta==exp(-kz); kz=K/Zmix - #} - alpha <- mleLoopPE(alpha=alpha, doobs=do.obs, c1=c1, c2=c2, beta=beta, irr=irr, wtr=wtr, kz=kz, dosat=do.sat) - - return(-sum(dnorm(do.obs, alpha, sd=sqrt(Q), log=TRUE), na.rm=TRUE)) + c1 <- exp(Params[1]) #PAR coeff + c2 <- -exp(Params[2]) #log(Temp) coeff + Q <- exp(Params[3]) # Variance of the process error + + # See KalmanDO_smooth.R comments for explanation of beta + kz <- k.gas/z.mix # K and Zmix are both vector of length nobs + beta <- exp(-kz) # This beta is for using the differential equation form + + # Set first true value equal to first observation + alpha <- rep(0, length(do.obs)) + alpha[1] <- do.obs[1]#Let's give this model some starting values + + #R version of C loop + #for(i in 2:length(do.obs)){ + # a1 <- c1*irr[i-1] + c2*log(wtr[i-1]) + kz[i-1]*do.sat[i-1] + # alpha[i] <- a1/kz[i-1] + -exp(-kz[i-1])*a1/kz[i-1] + beta[i-1]*alpha[i-1] # NOTE: beta==exp(-kz); kz=K/Zmix + #} + alpha <- mleLoopPE(alpha=alpha, doobs=do.obs, c1=c1, c2=c2, beta=beta, irr=irr, wtr=wtr, kz=kz, dosat=do.sat) + + return(-sum(dnorm(do.obs, alpha, sd=sqrt(Q), log=TRUE), na.rm=TRUE)) }#End function # ==================== # = mle NLL function = # ==================== mleNllOE <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type){ - c1 <- Params[1] #PAR coeff - c2 <- Params[2] #log(Temp) coeff - Q <- exp(Params[3]) # Variance of the process error - - # See KalmanDO_smooth.R comments for explanation of beta - kz <- k.gas/z.mix # K and Zmix are both vector of length nobs - beta <- exp(-kz) # This beta is for using the differential equation form - - # Set first true value equal to first observation - alpha <- rep(0, length(do.obs)) - alpha[1] <- Params[4] #Free varying initial DO value - - alpha <- mleLoopOE(alpha=alpha, doobs=do.obs, c1=c1, c2=c2, beta=beta, irr=irr, wtr=wtr, kz=kz, dosat=do.sat) - - return(-sum(dnorm(do.obs, alpha, sd=sqrt(Q), log=TRUE), na.rm=TRUE)) + c1 <- exp(Params[1]) #PAR coeff + c2 <- -exp(Params[2]) #log(Temp) coeff + Q <- exp(Params[3]) # Variance of the process error + + # See KalmanDO_smooth.R comments for explanation of beta + kz <- k.gas/z.mix # K and Zmix are both vector of length nobs + beta <- exp(-kz) # This beta is for using the differential equation form + + # Set first true value equal to first observation + alpha <- rep(0, length(do.obs)) + alpha[1] <- Params[4] #Free varying initial DO value + + alpha <- mleLoopOE(alpha=alpha, doobs=do.obs, c1=c1, c2=c2, beta=beta, irr=irr, wtr=wtr, kz=kz, dosat=do.sat) + + return(-sum(dnorm(do.obs, alpha, sd=sqrt(Q), log=TRUE), na.rm=TRUE)) }#End function From 52e08e01c1e3f5c04a9d272bc5ba8d8b23dc5c57 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Tue, 9 Jun 2026 14:11:59 -0400 Subject: [PATCH 08/20] Update indentation --- R/metab.mle.R | 152 +++++++++++++++++++++++++------------------------- 1 file changed, 76 insertions(+), 76 deletions(-) diff --git a/R/metab.mle.R b/R/metab.mle.R index 24779b6..754befa 100644 --- a/R/metab.mle.R +++ b/R/metab.mle.R @@ -12,9 +12,9 @@ #'@return #'A data.frame with columns corresponding to components of metabolism #'\describe{ - #'\item{GPP}{numeric estimate of Gross Primary Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} - #'\item{R}{numeric estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} - #'\item{NEP}{numeric estimate of Net Ecosystem production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{GPP}{numeric estimate of Gross Primary Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{R}{numeric estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} +#'\item{NEP}{numeric estimate of Net Ecosystem production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} #'} #' The maximum likelihood estimates of model parameters can be accessed via \code{attributes(metab.mle(...))[["params"]]} #' @@ -88,96 +88,96 @@ #'metab.mle(doobs[,2], do.sat, k.gas, z.mix[,2], irr[,2], wtr[,3]) #'@export metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", ...){ - + complete.inputs(do.obs=do.obs, do.sat=do.sat, k.gas=k.gas, z.mix=z.mix, irr=irr, wtr=wtr, error=TRUE) - - match.arg(error.type, choices=c('OE', 'PE')) - - nobs <- length(do.obs) - - mm.args <- list(...) - + + match.arg(error.type, choices=c('OE', 'PE')) + + nobs <- length(do.obs) + + mm.args <- list(...) + if(any(z.mix <= 0)){ stop("z.mix must be greater than zero.") } - if(any(wtr <= 0)){ - stop("all wtr must be positive.") - } - - if("datetime"%in%names(mm.args)){ # check to see if datetime is in the ... args - datetime <- mm.args$datetime # extract datetime - freq <- calc.freq(datetime) # calculate sampling frequency from datetime - if(nobs!=freq){ # nobs and freq should agree, if they don't issue a warning - bad.date <- format.Date(datetime[1], format="%Y-%m-%d") - warning("number of observations on ", bad.date, " (", nobs, ") ", "does not equal estimated sampling frequency", " (", freq, ")", sep="") - } - }else{ # if datetime is *not* in the ... args - warning("datetime not found, inferring sampling frequency from # of observations") # issue a warning (note checks in addNAs) - # NOTE: because of the checks in addNA's, it is unlikely a user would receive this warning via metab() - # warning will only be seen through direct use of metab.bookkeep when datettime is not supplied - freq <- nobs - } - - chk.list <- list(do.obs, irr, do.sat, z.mix, k.gas, wtr) - if(!all(sapply(chk.list, is.numeric)) || !all(sapply(chk.list, is.vector))){ - stop('All metab.mle inputs must be numeric vectors.') - } - - if(!all(nobs==sapply(chk.list, length))){ - stop('All input data to metab.mle must be the same length') - } - - Q0 <- ((diff(range(do.obs,na.rm=TRUE)) - mean(do.obs,na.rm=TRUE))^2 / length(do.obs)) - - guesses <- c(log(1E-4), log(1E-4), log(Q0)) - - #We have a different number of fitted parameters depending on error type of the model - if(error.type=='OE'){ - guesses <- c(guesses, do.obs[1]) - - fit <- optim(guesses, fn=mleNllOE, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) - - pars0 <- fit$par - pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "nll"=fit$value, "doInit"=pars0[4]) - - }else if(error.type=='PE'){ - guesses <- c(guesses) - - fit <- optim(guesses, fn=mleNllPE, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) - - pars0 <- fit$par - pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "nll"=fit$value) - - }else{ - stop("error.type must be either 'OE' or 'PE', Observation Error or Process Error respectively.") - } - - # ==================================== - # = Use fits to calculate metabolism = - # ==================================== - GPP <- mean(pars[1]*irr, na.rm=TRUE) * freq - R <- mean(pars[2]*log(wtr), na.rm=TRUE) * freq - - return(list("params"=pars, "metab"=c("GPP"=GPP,"R"=R,"NEP"=GPP+R))) + if(any(wtr <= 0)){ + stop("all wtr must be positive.") + } + + if("datetime"%in%names(mm.args)){ # check to see if datetime is in the ... args + datetime <- mm.args$datetime # extract datetime + freq <- calc.freq(datetime) # calculate sampling frequency from datetime + if(nobs!=freq){ # nobs and freq should agree, if they don't issue a warning + bad.date <- format.Date(datetime[1], format="%Y-%m-%d") + warning("number of observations on ", bad.date, " (", nobs, ") ", "does not equal estimated sampling frequency", " (", freq, ")", sep="") + } + }else{ # if datetime is *not* in the ... args + warning("datetime not found, inferring sampling frequency from # of observations") # issue a warning (note checks in addNAs) + # NOTE: because of the checks in addNA's, it is unlikely a user would receive this warning via metab() + # warning will only be seen through direct use of metab.bookkeep when datettime is not supplied + freq <- nobs + } + + chk.list <- list(do.obs, irr, do.sat, z.mix, k.gas, wtr) + if(!all(sapply(chk.list, is.numeric)) || !all(sapply(chk.list, is.vector))){ + stop('All metab.mle inputs must be numeric vectors.') + } + + if(!all(nobs==sapply(chk.list, length))){ + stop('All input data to metab.mle must be the same length') + } + + Q0 <- ((diff(range(do.obs,na.rm=TRUE)) - mean(do.obs,na.rm=TRUE))^2 / length(do.obs)) + + guesses <- c(log(1E-4), log(1E-4), log(Q0)) + + #We have a different number of fitted parameters depending on error type of the model + if(error.type=='OE'){ + guesses <- c(guesses, do.obs[1]) + + fit <- optim(guesses, fn=mleNllOE, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) + + pars0 <- fit$par + pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "nll"=fit$value, "doInit"=pars0[4]) + + }else if(error.type=='PE'){ + guesses <- c(guesses) + + fit <- optim(guesses, fn=mleNllPE, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) + + pars0 <- fit$par + pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "nll"=fit$value) + + }else{ + stop("error.type must be either 'OE' or 'PE', Observation Error or Process Error respectively.") + } + + # ==================================== + # = Use fits to calculate metabolism = + # ==================================== + GPP <- mean(pars[1]*irr, na.rm=TRUE) * freq + R <- mean(pars[2]*log(wtr), na.rm=TRUE) * freq + + return(list("params"=pars, "metab"=c("GPP"=GPP,"R"=R,"NEP"=GPP+R))) } # ============================================ # = The R loop for Observation Error mle NLL = # ============================================ mleLoopOE <- function(alpha, doobs, c1, c2, beta, irr, wtr, kz, dosat){ - nobs <- length(doobs) - a.loop <- .C("mleLoopCoe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer") - return(a.loop[["alpha"]]) + nobs <- length(doobs) + a.loop <- .C("mleLoopCoe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer") + return(a.loop[["alpha"]]) } # ============================================ # = The R loop for Process Error mle NLL = # ============================================ mleLoopPE <- function(alpha, doobs, c1, c2, beta, irr, wtr, kz, dosat){ - nobs <- length(doobs) - a.loop <- .C("mleLoopCpe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer") - return(a.loop[["alpha"]]) + nobs <- length(doobs) + a.loop <- .C("mleLoopCpe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer") + return(a.loop[["alpha"]]) } # ==================== From f924e5505aff2f6e84b47b9fa0497ceb81b06cd1 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Tue, 9 Jun 2026 15:40:55 -0400 Subject: [PATCH 09/20] Update k.klaus documentation --- R/k.klaus.R | 2 +- man/k.read.Rd | 2 +- man/k.read.base.Rd | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/R/k.klaus.R b/R/k.klaus.R index 8571f50..7aacb0a 100644 --- a/R/k.klaus.R +++ b/R/k.klaus.R @@ -59,4 +59,4 @@ k.klaus.base <- function(wnd, wnd.z, lake.area, sin, sdi = NULL, method = c("lin # -- References # Klaus, Marcus, and Dominic Vachon. 2020. # Challenges of Predicting Gas Transfer Velocity from Wind Measurements over Global Lakes. -# Aquatic Sciences 82 (3): 53. https://doi.org/10.1007/s00027-020-00729-9. +# Aquatic Sciences 82 (3): 53. https://doi.org/10.1007/s00027-020-00729-9 diff --git a/man/k.read.Rd b/man/k.read.Rd index e34ad21..34af434 100644 --- a/man/k.read.Rd +++ b/man/k.read.Rd @@ -44,7 +44,7 @@ k.klaus(ts.data, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", \item{sdi}{Numeric value of shoreline development index. Only for \link{k.klaus}'s "exp" method.} -\item{sin}{Numeric value >0 and ≤1 of the scale of spatial integration. Only for \link{k.klaus}.} +\item{sin}{Numeric value >0 and <1 of the scale of spatial integration. Only for \link{k.klaus}.} \item{params}{Only for \link{k.vachon.base} and \link{k.macIntyre}. See details.} } diff --git a/man/k.read.base.Rd b/man/k.read.base.Rd index 7a0cfa2..4efebed 100644 --- a/man/k.read.base.Rd +++ b/man/k.read.base.Rd @@ -61,7 +61,7 @@ k.klaus.base(wnd, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", \item{sdi}{Numeric value of shoreline development index. Only for \link{k.klaus.base}'s "exp" method.} -\item{sin}{Numeric value >0 and ≤1 of the scale of spatial integration. Only for \link{k.klaus.base}.} +\item{sin}{Numeric value >0 and <1 of the scale of spatial integration. Only for \link{k.klaus.base}.} \item{params}{Optional parameter input, only for \link{k.vachon.base} and \link{k.macIntyre.base}. See details.} } From ef0471aa8a73d251f9f8ccbe614f0bc6a54bafee Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Tue, 7 Jul 2026 18:03:36 -0400 Subject: [PATCH 10/20] Revert metab.mle changes; moving to separate PR --- R/metab.mle.R | 222 +++++++++++++++++++++++++------------------------- 1 file changed, 111 insertions(+), 111 deletions(-) diff --git a/R/metab.mle.R b/R/metab.mle.R index 754befa..4df7de1 100644 --- a/R/metab.mle.R +++ b/R/metab.mle.R @@ -12,9 +12,9 @@ #'@return #'A data.frame with columns corresponding to components of metabolism #'\describe{ -#'\item{GPP}{numeric estimate of Gross Primary Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} -#'\item{R}{numeric estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} -#'\item{NEP}{numeric estimate of Net Ecosystem production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} + #'\item{GPP}{numeric estimate of Gross Primary Production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} + #'\item{R}{numeric estimate of Respiration, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} + #'\item{NEP}{numeric estimate of Net Ecosystem production, \eqn{mg O_2 L^{-1} d^{-1}}{mg O2 / L / d}} #'} #' The maximum likelihood estimates of model parameters can be accessed via \code{attributes(metab.mle(...))[["params"]]} #' @@ -88,141 +88,141 @@ #'metab.mle(doobs[,2], do.sat, k.gas, z.mix[,2], irr[,2], wtr[,3]) #'@export metab.mle <- function(do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type="OE", ...){ - + complete.inputs(do.obs=do.obs, do.sat=do.sat, k.gas=k.gas, z.mix=z.mix, irr=irr, wtr=wtr, error=TRUE) - - match.arg(error.type, choices=c('OE', 'PE')) - - nobs <- length(do.obs) - - mm.args <- list(...) - + + match.arg(error.type, choices=c('OE', 'PE')) + + nobs <- length(do.obs) + + mm.args <- list(...) + if(any(z.mix <= 0)){ stop("z.mix must be greater than zero.") } - if(any(wtr <= 0)){ - stop("all wtr must be positive.") - } - - if("datetime"%in%names(mm.args)){ # check to see if datetime is in the ... args - datetime <- mm.args$datetime # extract datetime - freq <- calc.freq(datetime) # calculate sampling frequency from datetime - if(nobs!=freq){ # nobs and freq should agree, if they don't issue a warning - bad.date <- format.Date(datetime[1], format="%Y-%m-%d") - warning("number of observations on ", bad.date, " (", nobs, ") ", "does not equal estimated sampling frequency", " (", freq, ")", sep="") - } - }else{ # if datetime is *not* in the ... args - warning("datetime not found, inferring sampling frequency from # of observations") # issue a warning (note checks in addNAs) - # NOTE: because of the checks in addNA's, it is unlikely a user would receive this warning via metab() - # warning will only be seen through direct use of metab.bookkeep when datettime is not supplied - freq <- nobs - } - - chk.list <- list(do.obs, irr, do.sat, z.mix, k.gas, wtr) - if(!all(sapply(chk.list, is.numeric)) || !all(sapply(chk.list, is.vector))){ - stop('All metab.mle inputs must be numeric vectors.') - } - - if(!all(nobs==sapply(chk.list, length))){ - stop('All input data to metab.mle must be the same length') - } - - Q0 <- ((diff(range(do.obs,na.rm=TRUE)) - mean(do.obs,na.rm=TRUE))^2 / length(do.obs)) - - guesses <- c(log(1E-4), log(1E-4), log(Q0)) - - #We have a different number of fitted parameters depending on error type of the model - if(error.type=='OE'){ - guesses <- c(guesses, do.obs[1]) - - fit <- optim(guesses, fn=mleNllOE, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) - - pars0 <- fit$par - pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "nll"=fit$value, "doInit"=pars0[4]) - - }else if(error.type=='PE'){ - guesses <- c(guesses) - - fit <- optim(guesses, fn=mleNllPE, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) - - pars0 <- fit$par - pars <- c("gppCoeff"=exp(pars0[1]), "rCoeff"=-exp(pars0[2]), "Q"=exp(pars0[3]), "nll"=fit$value) - - }else{ - stop("error.type must be either 'OE' or 'PE', Observation Error or Process Error respectively.") - } - - # ==================================== - # = Use fits to calculate metabolism = - # ==================================== - GPP <- mean(pars[1]*irr, na.rm=TRUE) * freq - R <- mean(pars[2]*log(wtr), na.rm=TRUE) * freq - - return(list("params"=pars, "metab"=c("GPP"=GPP,"R"=R,"NEP"=GPP+R))) + if(any(wtr <= 0)){ + stop("all wtr must be positive.") + } + + if("datetime"%in%names(mm.args)){ # check to see if datetime is in the ... args + datetime <- mm.args$datetime # extract datetime + freq <- calc.freq(datetime) # calculate sampling frequency from datetime + if(nobs!=freq){ # nobs and freq should agree, if they don't issue a warning + bad.date <- format.Date(datetime[1], format="%Y-%m-%d") + warning("number of observations on ", bad.date, " (", nobs, ") ", "does not equal estimated sampling frequency", " (", freq, ")", sep="") + } + }else{ # if datetime is *not* in the ... args + warning("datetime not found, inferring sampling frequency from # of observations") # issue a warning (note checks in addNAs) + # NOTE: because of the checks in addNA's, it is unlikely a user would receive this warning via metab() + # warning will only be seen through direct use of metab.bookkeep when datettime is not supplied + freq <- nobs + } + + chk.list <- list(do.obs, irr, do.sat, z.mix, k.gas, wtr) + if(!all(sapply(chk.list, is.numeric)) || !all(sapply(chk.list, is.vector))){ + stop('All metab.mle inputs must be numeric vectors.') + } + + if(!all(nobs==sapply(chk.list, length))){ + stop('All input data to metab.mle must be the same length') + } + + Q0 <- ((diff(range(do.obs,na.rm=TRUE)) - mean(do.obs,na.rm=TRUE))^2 / length(do.obs)) + + guesses <- c(1E-4, 1E-4, log(Q0)) + + #We have a different number of fitted parameters depending on error type of the model + if(error.type=='OE'){ + guesses <- c(guesses, do.obs[1]) + + fit <- optim(guesses, fn=mleNllOE, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) + + pars0 <- fit$par + pars <- c("gppCoeff"=pars0[1], "rCoeff"=pars0[2], "Q"=exp(pars0[3]), "nll"=fit$value, "doInit"=pars0[4]) + + }else if(error.type=='PE'){ + guesses <- c(guesses) + + fit <- optim(guesses, fn=mleNllPE, do.obs=do.obs, do.sat=do.sat, k.gas=(k.gas/freq), z.mix=z.mix, irr=irr, wtr=wtr) + + pars0 <- fit$par + pars <- c("gppCoeff"=pars0[1], "rCoeff"=pars0[2], "Q"=exp(pars0[3]), "nll"=fit$value) + + }else{ + stop("error.type must be either 'OE' or 'PE', Observation Error or Process Error respectively.") + } + + # ==================================== + # = Use fits to calculate metabolism = + # ==================================== + GPP <- mean(pars[1]*irr, na.rm=TRUE) * freq + R <- mean(pars[2]*log(wtr), na.rm=TRUE) * freq + + return(list("params"=pars, "metab"=c("GPP"=GPP,"R"=R,"NEP"=GPP+R))) } # ============================================ # = The R loop for Observation Error mle NLL = # ============================================ mleLoopOE <- function(alpha, doobs, c1, c2, beta, irr, wtr, kz, dosat){ - nobs <- length(doobs) - a.loop <- .C("mleLoopCoe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer") - return(a.loop[["alpha"]]) + nobs <- length(doobs) + a.loop <- .C("mleLoopCoe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer") + return(a.loop[["alpha"]]) } # ============================================ # = The R loop for Process Error mle NLL = # ============================================ mleLoopPE <- function(alpha, doobs, c1, c2, beta, irr, wtr, kz, dosat){ - nobs <- length(doobs) - a.loop <- .C("mleLoopCpe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer") - return(a.loop[["alpha"]]) + nobs <- length(doobs) + a.loop <- .C("mleLoopCpe", alpha=as.double(alpha), as.double(doobs), as.double(c1), as.double(c2), as.double(beta), as.double(irr), as.double(wtr), as.double(kz), as.double(dosat), as.integer(nobs), PACKAGE="LakeMetabolizer") + return(a.loop[["alpha"]]) } # ==================== # = mle NLL function = # ==================== mleNllPE <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr){ - c1 <- exp(Params[1]) #PAR coeff - c2 <- -exp(Params[2]) #log(Temp) coeff - Q <- exp(Params[3]) # Variance of the process error - - # See KalmanDO_smooth.R comments for explanation of beta - kz <- k.gas/z.mix # K and Zmix are both vector of length nobs - beta <- exp(-kz) # This beta is for using the differential equation form - - # Set first true value equal to first observation - alpha <- rep(0, length(do.obs)) - alpha[1] <- do.obs[1]#Let's give this model some starting values - - #R version of C loop - #for(i in 2:length(do.obs)){ - # a1 <- c1*irr[i-1] + c2*log(wtr[i-1]) + kz[i-1]*do.sat[i-1] - # alpha[i] <- a1/kz[i-1] + -exp(-kz[i-1])*a1/kz[i-1] + beta[i-1]*alpha[i-1] # NOTE: beta==exp(-kz); kz=K/Zmix - #} - alpha <- mleLoopPE(alpha=alpha, doobs=do.obs, c1=c1, c2=c2, beta=beta, irr=irr, wtr=wtr, kz=kz, dosat=do.sat) - - return(-sum(dnorm(do.obs, alpha, sd=sqrt(Q), log=TRUE), na.rm=TRUE)) + c1 <- Params[1] #PAR coeff + c2 <- Params[2] #log(Temp) coeff + Q <- exp(Params[3]) # Variance of the process error + + # See KalmanDO_smooth.R comments for explanation of beta + kz <- k.gas/z.mix # K and Zmix are both vector of length nobs + beta <- exp(-kz) # This beta is for using the differential equation form + + # Set first true value equal to first observation + alpha <- rep(0, length(do.obs)) + alpha[1] <- do.obs[1]#Let's give this model some starting values + + #R version of C loop + #for(i in 2:length(do.obs)){ + # a1 <- c1*irr[i-1] + c2*log(wtr[i-1]) + kz[i-1]*do.sat[i-1] + # alpha[i] <- a1/kz[i-1] + -exp(-kz[i-1])*a1/kz[i-1] + beta[i-1]*alpha[i-1] # NOTE: beta==exp(-kz); kz=K/Zmix + #} + alpha <- mleLoopPE(alpha=alpha, doobs=do.obs, c1=c1, c2=c2, beta=beta, irr=irr, wtr=wtr, kz=kz, dosat=do.sat) + + return(-sum(dnorm(do.obs, alpha, sd=sqrt(Q), log=TRUE), na.rm=TRUE)) }#End function # ==================== # = mle NLL function = # ==================== mleNllOE <- function(Params, do.obs, do.sat, k.gas, z.mix, irr, wtr, error.type){ - c1 <- exp(Params[1]) #PAR coeff - c2 <- -exp(Params[2]) #log(Temp) coeff - Q <- exp(Params[3]) # Variance of the process error - - # See KalmanDO_smooth.R comments for explanation of beta - kz <- k.gas/z.mix # K and Zmix are both vector of length nobs - beta <- exp(-kz) # This beta is for using the differential equation form - - # Set first true value equal to first observation - alpha <- rep(0, length(do.obs)) - alpha[1] <- Params[4] #Free varying initial DO value - - alpha <- mleLoopOE(alpha=alpha, doobs=do.obs, c1=c1, c2=c2, beta=beta, irr=irr, wtr=wtr, kz=kz, dosat=do.sat) - - return(-sum(dnorm(do.obs, alpha, sd=sqrt(Q), log=TRUE), na.rm=TRUE)) + c1 <- Params[1] #PAR coeff + c2 <- Params[2] #log(Temp) coeff + Q <- exp(Params[3]) # Variance of the process error + + # See KalmanDO_smooth.R comments for explanation of beta + kz <- k.gas/z.mix # K and Zmix are both vector of length nobs + beta <- exp(-kz) # This beta is for using the differential equation form + + # Set first true value equal to first observation + alpha <- rep(0, length(do.obs)) + alpha[1] <- Params[4] #Free varying initial DO value + + alpha <- mleLoopOE(alpha=alpha, doobs=do.obs, c1=c1, c2=c2, beta=beta, irr=irr, wtr=wtr, kz=kz, dosat=do.sat) + + return(-sum(dnorm(do.obs, alpha, sd=sqrt(Q), log=TRUE), na.rm=TRUE)) }#End function From 22b8909a07ffd955ade83a4883bda47f85b53add Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 8 Jul 2026 08:53:19 -0400 Subject: [PATCH 11/20] Update k.klaus.base arguments in k.klaus Co-authored-by: Jake Zwart --- R/k.klaus.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/k.klaus.R b/R/k.klaus.R index 7aacb0a..e969d24 100644 --- a/R/k.klaus.R +++ b/R/k.klaus.R @@ -10,7 +10,7 @@ k.klaus = function(ts.data, wnd.z, lake.area, sin, sdi = NULL, method = c("linea wind = get.vars(ts.data, 'wnd') - k600 = k.klaus.base(wind[,2], wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", "power")) + k600 = k.klaus.base(wind[,2], wnd.z, lake.area, sin, sdi = sdi, method = method) return(data.frame(datetime=ts.data$datetime, k600=k600)) } From 6da9a69ffb150ebf4133b8f794146cbb8f8d0b2c Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 8 Jul 2026 08:54:13 -0400 Subject: [PATCH 12/20] Remove redundant qualifier from k.klaus Co-authored-by: Jake Zwart --- R/k.klaus.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/k.klaus.R b/R/k.klaus.R index e969d24..da5910b 100644 --- a/R/k.klaus.R +++ b/R/k.klaus.R @@ -26,7 +26,7 @@ k.klaus.base <- function(wnd, wnd.z, lake.area, sin, sdi = NULL, method = c("lin # Converting uz to u10 if (wnd.z != 10){ - wnd <- LakeMetabolizer::wind.scale.base(wnd = wnd, wnd.z = wnd.z) + wnd <- wind.scale.base(wnd = wnd, wnd.z = wnd.z) } # helper logit function From b2486162a0d4e90ca28248a199b81d923086ded9 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 8 Jul 2026 08:55:48 -0400 Subject: [PATCH 13/20] Rename logit_custom to logit.custom Co-authored-by: Jake Zwart --- R/k.klaus.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/k.klaus.R b/R/k.klaus.R index da5910b..3471beb 100644 --- a/R/k.klaus.R +++ b/R/k.klaus.R @@ -30,7 +30,7 @@ k.klaus.base <- function(wnd, wnd.z, lake.area, sin, sdi = NULL, method = c("lin } # helper logit function - logit_custom <- function(x){ + logit.custom <- function(x){ return(log(x / (1 - x))) } From fdf93ffccb408ea5c3a1eb3e4882fa59a305cb19 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 8 Jul 2026 08:56:20 -0400 Subject: [PATCH 14/20] Apply new logit.custom name Co-authored-by: Jake Zwart --- R/k.klaus.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/R/k.klaus.R b/R/k.klaus.R index 3471beb..a4814f8 100644 --- a/R/k.klaus.R +++ b/R/k.klaus.R @@ -40,14 +40,14 @@ k.klaus.base <- function(wnd, wnd.z, lake.area, sin, sdi = NULL, method = c("lin } if (method == "linear") { - k600 <- (0.328 * log10(lake.area) + 1.581) * wnd - 0.066 * logit_custom(sin) + 1.266 + k600 <- (0.328 * log10(lake.area) + 1.581) * wnd - 0.066 * logit.custom(sin) + 1.266 } else if (method == "power") { - k600 <- (0.281 * log10(lake.area) + 1.361) * (wnd ^ 1.097) - 0.072 * logit_custom(sin) + 1.401 + k600 <- (0.281 * log10(lake.area) + 1.361) * (wnd ^ 1.097) - 0.072 * logit.custom(sin) + 1.401 } else if (method == "exp") { if (is.null(sdi)) { stop("sdi must be provided for exponential model") } - k600 <- (-0.057 * logit_custom(sin) + 2.366) * exp(wnd * (0.144 * log10(sdi) + 0.156)) + k600 <- (-0.057 * logit.custom(sin) + 2.366) * exp(wnd * (0.144 * log10(sdi) + 0.156)) } k600 <- k600 * (24/100) # convert cm/hr to m/day From aa559e36fcb9342bb3382beeff9272b6092bb2f4 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 8 Jul 2026 08:57:44 -0400 Subject: [PATCH 15/20] Update k.klaus documentation in k.read.R Co-authored-by: Jake Zwart --- R/k.read.R | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/R/k.read.R b/R/k.read.R index af33e70..ef28b59 100644 --- a/R/k.read.R +++ b/R/k.read.R @@ -34,8 +34,13 @@ #'@param atm.press atmospheric pressure in mb #'@param lat Latitude, degrees north #'@param lake.area Lake area, m^2 -#'@param sdi Numeric value of shoreline development index. Only for \link{k.klaus}'s "exp" method. -#'@param sin Numeric value >0 and <1 of the scale of spatial integration. Only for \link{k.klaus}. +#'@param sdi Only for \link{k.klaus}'s "exp" method. Numeric shoreline development index +#'(unitless, >= 1): the ratio of shoreline length to the circumference of a circle of equal +#'lake area, SDI = P / (2 * sqrt(pi * lake.area)). +#'@param sin Only for \link{k.klaus}. Numeric scale of spatial integration, expressed relative +#'to total lake surface area: from near 0 (a point-scale measurement, e.g. 1 m^2 / lake.area) +#'up to but not including 1 (whole-lake integration). Must satisfy 0 < sin < 1; the value 1 is +#'not supported because the model uses logit(sin). #'@param params Only for \link{k.vachon.base} and \link{k.macIntyre}. See details. #' #'@details Can change default parameters of MacIntyre and Vachon models. Default for Vachon is From 558a48c440a6f0941cb0cab8a478965144b78c85 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 8 Jul 2026 08:58:16 -0400 Subject: [PATCH 16/20] Update k.klaus documentation in k.read.R Co-authored-by: Jake Zwart --- R/k.read.R | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/R/k.read.R b/R/k.read.R index ef28b59..99d608e 100644 --- a/R/k.read.R +++ b/R/k.read.R @@ -232,8 +232,13 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #'@param RH Numeric value of relative humidity, \% #'@param sw Numeric value of short wave radiation, W m^-2 #'@param lwnet Numeric value net long wave radiation, W m^-2 -#'@param sdi Numeric value of shoreline development index. Only for \link{k.klaus.base}'s "exp" method. -#'@param sin Numeric value >0 and <1 of the scale of spatial integration. Only for \link{k.klaus.base}. +#'@param sdi Only for \link{k.klaus}'s "exp" method. Numeric shoreline development index +#'(unitless, >= 1): the ratio of shoreline length to the circumference of a circle of equal +#'lake area, SDI = P / (2 * sqrt(pi * lake.area)). +#'@param sin Only for \link{k.klaus}. Numeric scale of spatial integration, expressed relative +#'to total lake surface area: from near 0 (a point-scale measurement, e.g. 1 m^2 / lake.area) +#'up to but not including 1 (whole-lake integration). Must satisfy 0 < sin < 1; the value 1 is +#'not supported because the model uses logit(sin). #'@param params Optional parameter input, only for \link{k.vachon.base} and \link{k.macIntyre.base}. See details. #'@details Can change default parameters of MacIntyre and Vachon models. Default for Vachon is #'c(2.51,1.48,0.39). Default for MacIntyre is c(1.2,0.4872,1.4784). Heiskanen et al. (2014) uses MacIntyre From 80affb8c0d037c74cac51aa2592c98483e9d1336 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 8 Jul 2026 09:16:51 -0400 Subject: [PATCH 17/20] Rename sin to spatial.int and update k.klaus documentation --- R/k.klaus.R | 16 ++++++++-------- R/k.read.R | 46 ++++++++++++++++++++++++++++++++-------------- man/k.read.Rd | 26 ++++++++++++++++++++------ man/k.read.base.Rd | 26 ++++++++++++++++++++------ 4 files changed, 80 insertions(+), 34 deletions(-) diff --git a/R/k.klaus.R b/R/k.klaus.R index a4814f8..d42f5c6 100644 --- a/R/k.klaus.R +++ b/R/k.klaus.R @@ -2,7 +2,7 @@ # Last update: 2026-06-08 #'@export -k.klaus = function(ts.data, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", "power")){ +k.klaus = function(ts.data, wnd.z, lake.area, spatial.int, sdi = NULL, method = c("linear", "exp", "power")){ if(!has.vars(ts.data, 'wnd')){ stop('k.klaus requires a "wnd" (wind speed) column in the supplied data') @@ -10,14 +10,14 @@ k.klaus = function(ts.data, wnd.z, lake.area, sin, sdi = NULL, method = c("linea wind = get.vars(ts.data, 'wnd') - k600 = k.klaus.base(wind[,2], wnd.z, lake.area, sin, sdi = sdi, method = method) + k600 = k.klaus.base(wind[,2], wnd.z, lake.area, spatial.int, sdi = sdi, method = method) return(data.frame(datetime=ts.data$datetime, k600=k600)) } #'@export -k.klaus.base <- function(wnd, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", "power")) { +k.klaus.base <- function(wnd, wnd.z, lake.area, spatial.int, sdi = NULL, method = c("linear", "exp", "power")) { method <- match.arg(method) @@ -35,19 +35,19 @@ k.klaus.base <- function(wnd, wnd.z, lake.area, sin, sdi = NULL, method = c("lin } # sanity checks - if (any(sin <= 0 | sin >= 1)) { - stop("sin must be between 0 and 1 (exclusive)") + if (any(spatial.int <= 0 | spatial.int >= 1)) { + stop("spatial.int must be between 0 and 1 (exclusive)") } if (method == "linear") { - k600 <- (0.328 * log10(lake.area) + 1.581) * wnd - 0.066 * logit.custom(sin) + 1.266 + k600 <- (0.328 * log10(lake.area) + 1.581) * wnd - 0.066 * logit.custom(spatial.int) + 1.266 } else if (method == "power") { - k600 <- (0.281 * log10(lake.area) + 1.361) * (wnd ^ 1.097) - 0.072 * logit.custom(sin) + 1.401 + k600 <- (0.281 * log10(lake.area) + 1.361) * (wnd ^ 1.097) - 0.072 * logit.custom(spatial.int) + 1.401 } else if (method == "exp") { if (is.null(sdi)) { stop("sdi must be provided for exponential model") } - k600 <- (-0.057 * logit.custom(sin) + 2.366) * exp(wnd * (0.144 * log10(sdi) + 0.156)) + k600 <- (-0.057 * logit.custom(spatial.int) + 2.366) * exp(wnd * (0.144 * log10(sdi) + 0.156)) } k600 <- k600 * (24/100) # convert cm/hr to m/day diff --git a/R/k.read.R b/R/k.read.R index 99d608e..c6ee9d8 100644 --- a/R/k.read.R +++ b/R/k.read.R @@ -26,9 +26,9 @@ #' #'k.heiskanen(ts.data, wnd.z, Kd, atm.press) #' -#'k.klaus(ts.data, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", "power")) +#'k.klaus(ts.data, wnd.z, lake.area, spatial.int, sdi = NULL, method = c("linear", "exp", "power")) #'@param ts.data vector of datetime in POSIXct format -#'@param method Only for \link{k.crusius} and \link{k.klaus}. String of valid method . For k.crusius either "linear", "bilinear", or "power". For k.klaus either "linear", "exp", or "power" +#'@param method Only for \link{k.crusius} and \link{k.klaus}. String of valid method. For k.crusius either "linear", "bilinear", or "power". For k.klaus either "linear", "exp", or "power" #'@param wnd.z height of wind measurement, m #'@param Kd Light attenuation coefficient (Units:m^-1) #'@param atm.press atmospheric pressure in mb @@ -37,16 +37,25 @@ #'@param sdi Only for \link{k.klaus}'s "exp" method. Numeric shoreline development index #'(unitless, >= 1): the ratio of shoreline length to the circumference of a circle of equal #'lake area, SDI = P / (2 * sqrt(pi * lake.area)). -#'@param sin Only for \link{k.klaus}. Numeric scale of spatial integration, expressed relative +#'@param spatial.int Only for \link{k.klaus}. Numeric scale of spatial integration, expressed relative #'to total lake surface area: from near 0 (a point-scale measurement, e.g. 1 m^2 / lake.area) -#'up to but not including 1 (whole-lake integration). Must satisfy 0 < sin < 1; the value 1 is -#'not supported because the model uses logit(sin). +#'up to but not including 1 (whole-lake integration). Must satisfy 0 < spatial.int < 1; the value 1 is +#'not supported because the model uses logit(spatial.int). #'@param params Only for \link{k.vachon.base} and \link{k.macIntyre}. See details. #' #'@details Can change default parameters of MacIntyre and Vachon models. Default for Vachon is #'c(2.51,1.48,0.39). Default for MacIntyre is c(1.2,0.4872,1.4784). Heiskanen 2014 uses MacIntyre #'model with c(0.5,0.77,0.3) and z.aml constant at 0.15. #' +#'\code{k.klaus} implements the three wind-based models of Klaus and Vachon (2020), +#'each fit to 46 globally distributed lakes (their Table 3, lowest-AIC parameterizations). +#'The \code{method} argument selects the functional form of the k600-wind relationship: +#'"linear" and "power" scale k600 with wind speed, lake area, and the scale of spatial +#'integration (\code{spatial.int}); "exp" is exponential in wind speed and additionally depends on +#'the shoreline development index (\code{sdi}). Note: Klaus and Vachon (2020) caution that +#'the SDI-based ("exp") model is not mechanistically meaningful and should not be used for +#'predictive purposes. +#' #'@return Returns a data.frame with a datetime column and a k600 column. k600 is in units of meters per day (m/d). #'@import rLakeAnalyzer #'@useDynLib LakeMetabolizer, .registration = TRUE @@ -115,7 +124,7 @@ #'atm.press = 1018 #'lat = tb.data$metadata$latitude #'lake.area = tb.data$metadata$lakearea -#'sin = 1/lake.area +#'spatial.int = 1/lake.area #' #'#for k.read and k.macIntyre, we need LW_net. #'#Calculate from the observations we have available. @@ -131,7 +140,7 @@ #' #'k600_macIntyre = k.macIntyre(ts.data, wnd.z=wnd.z, Kd=kd, atm.press=atm.press) #' -#'k600_klaus = k.klaus(ts.data, wnd.z=wnd.z, lake.area=lake.area, sin=sin) +#'k600_klaus = k.klaus(ts.data, wnd.z=wnd.z, lake.area=lake.area, spatial.int=spatial.int) #' #'@export k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ @@ -217,9 +226,9 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #' #'k.heiskanen.base(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet) #' -#'k.klaus.base(wnd, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", "power")) +#'k.klaus.base(wnd, wnd.z, lake.area, spatial.int, sdi = NULL, method = c("linear", "exp", "power")) #'@param wnd Numeric value of wind speed, (Units:m/s) -#'@param method Only for \link{k.crusius.base} and \link{k.klaus.base}. String of valid method . For k.crusius.base either "constant", "bilinear", or "power". For k.klaus.base either "linear", "exp", or "power" +#'@param method Only for \link{k.crusius.base} and \link{k.klaus.base}. String of valid method. For k.crusius.base either "constant", "bilinear", or "power". For k.klaus.base either "linear", "exp", or "power" #'@param wnd.z Height of wind measurement, (Units: m) #'@param Kd Light attenuation coefficient (Units: m^-1) #'@param lat Latitude, degrees north @@ -235,14 +244,23 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #'@param sdi Only for \link{k.klaus}'s "exp" method. Numeric shoreline development index #'(unitless, >= 1): the ratio of shoreline length to the circumference of a circle of equal #'lake area, SDI = P / (2 * sqrt(pi * lake.area)). -#'@param sin Only for \link{k.klaus}. Numeric scale of spatial integration, expressed relative +#'@param spatial.int Only for \link{k.klaus}. Numeric scale of spatial integration, expressed relative #'to total lake surface area: from near 0 (a point-scale measurement, e.g. 1 m^2 / lake.area) -#'up to but not including 1 (whole-lake integration). Must satisfy 0 < sin < 1; the value 1 is -#'not supported because the model uses logit(sin). +#'up to but not including 1 (whole-lake integration). Must satisfy 0 < spatial.int < 1; the value 1 is +#'not supported because the model uses logit(spatial.int). #'@param params Optional parameter input, only for \link{k.vachon.base} and \link{k.macIntyre.base}. See details. #'@details Can change default parameters of MacIntyre and Vachon models. Default for Vachon is #'c(2.51,1.48,0.39). Default for MacIntyre is c(1.2,0.4872,1.4784). Heiskanen et al. (2014) uses MacIntyre #'model with c(0.5,0.77,0.3) and z.aml constant at 0.15. +#' +#'\code{k.klaus} implements the three wind-based models of Klaus and Vachon (2020), +#'each fit to 46 globally distributed lakes (their Table 3, lowest-AIC parameterizations). +#'The \code{method} argument selects the functional form of the k600-wind relationship: +#'"linear" and "power" scale k600 with wind speed, lake area, and the scale of spatial +#'integration (\code{spatial.int}); "exp" is exponential in wind speed and additionally depends on +#'the shoreline development index (\code{sdi}). Note: Klaus and Vachon (2020) caution that +#'the SDI-based ("exp") model is not mechanistically meaningful and should not be used for +#'predictive purposes #'@return Numeric value of gas exchange velocity (k600) in units of m/day. Before use, #'should be converted to appropriate gas using \link{k600.2.kGAS}. #'@keywords methods math @@ -301,7 +319,7 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #'RH <- 90 #'sw <- 800 #'lwnet <- -55 -#'sin <- 1/lake.area +#'spatial.int <- 1/lake.area #'timeStep <- 30 #' #'U10 <- wind.scale.base(wnd, wnd.z) @@ -319,7 +337,7 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #'k600_macInytre <- k.macIntyre.base(wnd.z, Kd, atm.press, #'dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet) #' -#'k600_klaus <- k.klaus.base(wnd, wnd.z, lake.area, sin) +#'k600_klaus <- k.klaus.base(wnd, wnd.z, lake.area, spatial.int) #' #'@export k.read.base <- function(wnd.z, Kd, lat, lake.area, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet){ diff --git a/man/k.read.Rd b/man/k.read.Rd index 34af434..1fe4ebd 100644 --- a/man/k.read.Rd +++ b/man/k.read.Rd @@ -25,7 +25,7 @@ k.vachon(ts.data, lake.area, params=c(2.51,1.48,0.39)) k.heiskanen(ts.data, wnd.z, Kd, atm.press) -k.klaus(ts.data, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", "power")) +k.klaus(ts.data, wnd.z, lake.area, spatial.int, sdi = NULL, method = c("linear", "exp", "power")) } \arguments{ \item{ts.data}{vector of datetime in POSIXct format} @@ -40,11 +40,16 @@ k.klaus(ts.data, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", \item{lake.area}{Lake area, m^2} -\item{method}{Only for \link{k.crusius} and \link{k.klaus}. String of valid method . For k.crusius either "linear", "bilinear", or "power". For k.klaus either "linear", "exp", or "power"} +\item{method}{Only for \link{k.crusius} and \link{k.klaus}. String of valid method. For k.crusius either "linear", "bilinear", or "power". For k.klaus either "linear", "exp", or "power"} -\item{sdi}{Numeric value of shoreline development index. Only for \link{k.klaus}'s "exp" method.} +\item{sdi}{Only for \link{k.klaus}'s "exp" method. Numeric shoreline development index +(unitless, >= 1): the ratio of shoreline length to the circumference of a circle of equal +lake area, SDI = P / (2 * sqrt(pi * lake.area)).} -\item{sin}{Numeric value >0 and <1 of the scale of spatial integration. Only for \link{k.klaus}.} +\item{spatial.int}{Only for \link{k.klaus}. Numeric scale of spatial integration, expressed relative +to total lake surface area: from near 0 (a point-scale measurement, e.g. 1 m^2 / lake.area) +up to but not including 1 (whole-lake integration). Must satisfy 0 < spatial.int < 1; the value 1 is +not supported because the model uses logit(spatial.int).} \item{params}{Only for \link{k.vachon.base} and \link{k.macIntyre}. See details.} } @@ -58,6 +63,15 @@ Returns the gas exchange velocity based on the chosen model in units of m/day Can change default parameters of MacIntyre and Vachon models. Default for Vachon is c(2.51,1.48,0.39). Default for MacIntyre is c(1.2,0.4872,1.4784). Heiskanen 2014 uses MacIntyre model with c(0.5,0.77,0.3) and z.aml constant at 0.15. + +\code{k.klaus} implements the three wind-based models of Klaus and Vachon (2020), +each fit to 46 globally distributed lakes (their Table 3, lowest-AIC parameterizations). +The \code{method} argument selects the functional form of the k600-wind relationship: +"linear" and "power" scale k600 with wind speed, lake area, and the scale of spatial +integration (\code{spatial.int}); "exp" is exponential in wind speed and additionally depends on +the shoreline development index (\code{sdi}). Note: Klaus and Vachon (2020) caution that +the SDI-based ("exp") model is not mechanistically meaningful and should not be used for +predictive purposes. } \examples{ data.path = system.file('extdata', package="LakeMetabolizer") @@ -82,7 +96,7 @@ wnd.z = 10 #because we converted to u10 atm.press = 1018 lat = tb.data$metadata$latitude lake.area = tb.data$metadata$lakearea -sin = 1/lake.area +spatial.int = 1/lake.area #for k.read and k.macIntyre, we need LW_net. #Calculate from the observations we have available. @@ -98,7 +112,7 @@ atm.press=atm.press, lat=lat, lake.area=lake.area) k600_macIntyre = k.macIntyre(ts.data, wnd.z=wnd.z, Kd=kd, atm.press=atm.press) -k600_klaus = k.klaus(ts.data, wnd.z=wnd.z, lake.area=lake.area, sin=sin) +k600_klaus = k.klaus(ts.data, wnd.z=wnd.z, lake.area=lake.area, spatial.int=spatial.int) } \references{ diff --git a/man/k.read.base.Rd b/man/k.read.base.Rd index 4efebed..1c595a3 100644 --- a/man/k.read.base.Rd +++ b/man/k.read.base.Rd @@ -28,7 +28,7 @@ k.vachon.base(wnd, lake.area, params=c(2.51,1.48,0.39)) k.heiskanen.base(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet) -k.klaus.base(wnd, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", "power")) +k.klaus.base(wnd, wnd.z, lake.area, spatial.int, sdi = NULL, method = c("linear", "exp", "power")) } \arguments{ \item{wnd.z}{Height of wind measurement, (Units: m)} @@ -57,11 +57,16 @@ k.klaus.base(wnd, wnd.z, lake.area, sin, sdi = NULL, method = c("linear", "exp", \item{lwnet}{Numeric value net long wave radiation, W m^-2} -\item{method}{Only for \link{k.crusius.base} and \link{k.klaus.base}. String of valid method . For k.crusius.base either "constant", "bilinear", or "power". For k.klaus.base either "linear", "exp", or "power"} +\item{method}{Only for \link{k.crusius.base} and \link{k.klaus.base}. String of valid method. For k.crusius.base either "constant", "bilinear", or "power". For k.klaus.base either "linear", "exp", or "power"} -\item{sdi}{Numeric value of shoreline development index. Only for \link{k.klaus.base}'s "exp" method.} +\item{sdi}{Only for \link{k.klaus}'s "exp" method. Numeric shoreline development index +(unitless, >= 1): the ratio of shoreline length to the circumference of a circle of equal +lake area, SDI = P / (2 * sqrt(pi * lake.area)).} -\item{sin}{Numeric value >0 and <1 of the scale of spatial integration. Only for \link{k.klaus.base}.} +\item{spatial.int}{Only for \link{k.klaus}. Numeric scale of spatial integration, expressed relative +to total lake surface area: from near 0 (a point-scale measurement, e.g. 1 m^2 / lake.area) +up to but not including 1 (whole-lake integration). Must satisfy 0 < spatial.int < 1; the value 1 is +not supported because the model uses logit(spatial.int).} \item{params}{Optional parameter input, only for \link{k.vachon.base} and \link{k.macIntyre.base}. See details.} } @@ -76,6 +81,15 @@ Returns the gas exchange velocity based on the chosen model in units of m/day Can change default parameters of MacIntyre and Vachon models. Default for Vachon is c(2.51,1.48,0.39). Default for MacIntyre is c(1.2,0.4872,1.4784). Heiskanen et al. (2014) uses MacIntyre model with c(0.5,0.77,0.3) and z.aml constant at 0.15. + +\code{k.klaus} implements the three wind-based models of Klaus and Vachon (2020), +each fit to 46 globally distributed lakes (their Table 3, lowest-AIC parameterizations). +The \code{method} argument selects the functional form of the k600-wind relationship: +"linear" and "power" scale k600 with wind speed, lake area, and the scale of spatial +integration (\code{spatial.int}); "exp" is exponential in wind speed and additionally depends on +the shoreline development index (\code{sdi}). Note: Klaus and Vachon (2020) caution that +the SDI-based ("exp") model is not mechanistically meaningful and should not be used for +predictive purposes } \examples{ wnd.z <- 2 @@ -91,7 +105,7 @@ wnd <- 6 RH <- 90 sw <- 800 lwnet <- -55 -sin <- 1/lake.area +spatial.int <- 1/lake.area timeStep <- 30 U10 <- wind.scale.base(wnd, wnd.z) @@ -109,7 +123,7 @@ atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet) k600_macInytre <- k.macIntyre.base(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet) -k600_klaus <- k.klaus.base(wnd, wnd.z, lake.area, sin) +k600_klaus <- k.klaus.base(wnd, wnd.z, lake.area, spatial.int) } \references{ From 1a577838a1e6a1e1ea72909b4705d6e27d408d21 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 8 Jul 2026 09:41:31 -0400 Subject: [PATCH 18/20] Removed exp model from k.klaus --- R/k.klaus.R | 16 +++++++++------- R/k.read.R | 29 +++-------------------------- 2 files changed, 12 insertions(+), 33 deletions(-) diff --git a/R/k.klaus.R b/R/k.klaus.R index d42f5c6..0a6aa46 100644 --- a/R/k.klaus.R +++ b/R/k.klaus.R @@ -2,7 +2,7 @@ # Last update: 2026-06-08 #'@export -k.klaus = function(ts.data, wnd.z, lake.area, spatial.int, sdi = NULL, method = c("linear", "exp", "power")){ +k.klaus = function(ts.data, wnd.z, lake.area, spatial.int, method = c("linear", "power")){ if(!has.vars(ts.data, 'wnd')){ stop('k.klaus requires a "wnd" (wind speed) column in the supplied data') @@ -10,14 +10,14 @@ k.klaus = function(ts.data, wnd.z, lake.area, spatial.int, sdi = NULL, method = wind = get.vars(ts.data, 'wnd') - k600 = k.klaus.base(wind[,2], wnd.z, lake.area, spatial.int, sdi = sdi, method = method) + k600 = k.klaus.base(wind[,2], wnd.z, lake.area, spatial.int, method = method) return(data.frame(datetime=ts.data$datetime, k600=k600)) } #'@export -k.klaus.base <- function(wnd, wnd.z, lake.area, spatial.int, sdi = NULL, method = c("linear", "exp", "power")) { +k.klaus.base <- function(wnd, wnd.z, lake.area, spatial.int, method = c("linear", "power")) { method <- match.arg(method) @@ -44,10 +44,12 @@ k.klaus.base <- function(wnd, wnd.z, lake.area, spatial.int, sdi = NULL, method } else if (method == "power") { k600 <- (0.281 * log10(lake.area) + 1.361) * (wnd ^ 1.097) - 0.072 * logit.custom(spatial.int) + 1.401 } else if (method == "exp") { - if (is.null(sdi)) { - stop("sdi must be provided for exponential model") - } - k600 <- (-0.057 * logit.custom(spatial.int) + 2.366) * exp(wnd * (0.144 * log10(sdi) + 0.156)) + # if (is.null(sdi)) { + # stop("sdi must be provided for exponential model") + # } + # k600 <- (-0.057 * logit.custom(spatial.int) + 2.366) * exp(wnd * (0.144 * log10(sdi) + 0.156)) + stop('The exponential model was removed per Klaus and Vachon (2020): "...[We] do not recommend using our +models that include SDI."') } k600 <- k600 * (24/100) # convert cm/hr to m/day diff --git a/R/k.read.R b/R/k.read.R index c6ee9d8..c321af7 100644 --- a/R/k.read.R +++ b/R/k.read.R @@ -26,17 +26,14 @@ #' #'k.heiskanen(ts.data, wnd.z, Kd, atm.press) #' -#'k.klaus(ts.data, wnd.z, lake.area, spatial.int, sdi = NULL, method = c("linear", "exp", "power")) +#'k.klaus(ts.data, wnd.z, lake.area, spatial.int, method = c("linear", "exp", "power")) #'@param ts.data vector of datetime in POSIXct format -#'@param method Only for \link{k.crusius} and \link{k.klaus}. String of valid method. For k.crusius either "linear", "bilinear", or "power". For k.klaus either "linear", "exp", or "power" +#'@param method Only for \link{k.crusius} and \link{k.klaus}. String of valid method. For k.crusius either "linear", "bilinear", or "power". For k.klaus either "linear" or "power" #'@param wnd.z height of wind measurement, m #'@param Kd Light attenuation coefficient (Units:m^-1) #'@param atm.press atmospheric pressure in mb #'@param lat Latitude, degrees north #'@param lake.area Lake area, m^2 -#'@param sdi Only for \link{k.klaus}'s "exp" method. Numeric shoreline development index -#'(unitless, >= 1): the ratio of shoreline length to the circumference of a circle of equal -#'lake area, SDI = P / (2 * sqrt(pi * lake.area)). #'@param spatial.int Only for \link{k.klaus}. Numeric scale of spatial integration, expressed relative #'to total lake surface area: from near 0 (a point-scale measurement, e.g. 1 m^2 / lake.area) #'up to but not including 1 (whole-lake integration). Must satisfy 0 < spatial.int < 1; the value 1 is @@ -47,15 +44,6 @@ #'c(2.51,1.48,0.39). Default for MacIntyre is c(1.2,0.4872,1.4784). Heiskanen 2014 uses MacIntyre #'model with c(0.5,0.77,0.3) and z.aml constant at 0.15. #' -#'\code{k.klaus} implements the three wind-based models of Klaus and Vachon (2020), -#'each fit to 46 globally distributed lakes (their Table 3, lowest-AIC parameterizations). -#'The \code{method} argument selects the functional form of the k600-wind relationship: -#'"linear" and "power" scale k600 with wind speed, lake area, and the scale of spatial -#'integration (\code{spatial.int}); "exp" is exponential in wind speed and additionally depends on -#'the shoreline development index (\code{sdi}). Note: Klaus and Vachon (2020) caution that -#'the SDI-based ("exp") model is not mechanistically meaningful and should not be used for -#'predictive purposes. -#' #'@return Returns a data.frame with a datetime column and a k600 column. k600 is in units of meters per day (m/d). #'@import rLakeAnalyzer #'@useDynLib LakeMetabolizer, .registration = TRUE @@ -226,7 +214,7 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #' #'k.heiskanen.base(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet) #' -#'k.klaus.base(wnd, wnd.z, lake.area, spatial.int, sdi = NULL, method = c("linear", "exp", "power")) +#'k.klaus.base(wnd, wnd.z, lake.area, spatial.int, method = c("linear", "exp", "power")) #'@param wnd Numeric value of wind speed, (Units:m/s) #'@param method Only for \link{k.crusius.base} and \link{k.klaus.base}. String of valid method. For k.crusius.base either "constant", "bilinear", or "power". For k.klaus.base either "linear", "exp", or "power" #'@param wnd.z Height of wind measurement, (Units: m) @@ -241,9 +229,6 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #'@param RH Numeric value of relative humidity, \% #'@param sw Numeric value of short wave radiation, W m^-2 #'@param lwnet Numeric value net long wave radiation, W m^-2 -#'@param sdi Only for \link{k.klaus}'s "exp" method. Numeric shoreline development index -#'(unitless, >= 1): the ratio of shoreline length to the circumference of a circle of equal -#'lake area, SDI = P / (2 * sqrt(pi * lake.area)). #'@param spatial.int Only for \link{k.klaus}. Numeric scale of spatial integration, expressed relative #'to total lake surface area: from near 0 (a point-scale measurement, e.g. 1 m^2 / lake.area) #'up to but not including 1 (whole-lake integration). Must satisfy 0 < spatial.int < 1; the value 1 is @@ -253,14 +238,6 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #'c(2.51,1.48,0.39). Default for MacIntyre is c(1.2,0.4872,1.4784). Heiskanen et al. (2014) uses MacIntyre #'model with c(0.5,0.77,0.3) and z.aml constant at 0.15. #' -#'\code{k.klaus} implements the three wind-based models of Klaus and Vachon (2020), -#'each fit to 46 globally distributed lakes (their Table 3, lowest-AIC parameterizations). -#'The \code{method} argument selects the functional form of the k600-wind relationship: -#'"linear" and "power" scale k600 with wind speed, lake area, and the scale of spatial -#'integration (\code{spatial.int}); "exp" is exponential in wind speed and additionally depends on -#'the shoreline development index (\code{sdi}). Note: Klaus and Vachon (2020) caution that -#'the SDI-based ("exp") model is not mechanistically meaningful and should not be used for -#'predictive purposes #'@return Numeric value of gas exchange velocity (k600) in units of m/day. Before use, #'should be converted to appropriate gas using \link{k600.2.kGAS}. #'@keywords methods math From 007e1d9eb1f1f7fd9ceacc8d321372857ed902b0 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 8 Jul 2026 09:45:51 -0400 Subject: [PATCH 19/20] Rebuild documentation --- man/k.read.Rd | 17 ++--------------- man/k.read.base.Rd | 15 +-------------- 2 files changed, 3 insertions(+), 29 deletions(-) diff --git a/man/k.read.Rd b/man/k.read.Rd index 1fe4ebd..266e709 100644 --- a/man/k.read.Rd +++ b/man/k.read.Rd @@ -25,7 +25,7 @@ k.vachon(ts.data, lake.area, params=c(2.51,1.48,0.39)) k.heiskanen(ts.data, wnd.z, Kd, atm.press) -k.klaus(ts.data, wnd.z, lake.area, spatial.int, sdi = NULL, method = c("linear", "exp", "power")) +k.klaus(ts.data, wnd.z, lake.area, spatial.int, method = c("linear", "exp", "power")) } \arguments{ \item{ts.data}{vector of datetime in POSIXct format} @@ -40,11 +40,7 @@ k.klaus(ts.data, wnd.z, lake.area, spatial.int, sdi = NULL, method = c("linear", \item{lake.area}{Lake area, m^2} -\item{method}{Only for \link{k.crusius} and \link{k.klaus}. String of valid method. For k.crusius either "linear", "bilinear", or "power". For k.klaus either "linear", "exp", or "power"} - -\item{sdi}{Only for \link{k.klaus}'s "exp" method. Numeric shoreline development index -(unitless, >= 1): the ratio of shoreline length to the circumference of a circle of equal -lake area, SDI = P / (2 * sqrt(pi * lake.area)).} +\item{method}{Only for \link{k.crusius} and \link{k.klaus}. String of valid method. For k.crusius either "linear", "bilinear", or "power". For k.klaus either "linear" or "power"} \item{spatial.int}{Only for \link{k.klaus}. Numeric scale of spatial integration, expressed relative to total lake surface area: from near 0 (a point-scale measurement, e.g. 1 m^2 / lake.area) @@ -63,15 +59,6 @@ Returns the gas exchange velocity based on the chosen model in units of m/day Can change default parameters of MacIntyre and Vachon models. Default for Vachon is c(2.51,1.48,0.39). Default for MacIntyre is c(1.2,0.4872,1.4784). Heiskanen 2014 uses MacIntyre model with c(0.5,0.77,0.3) and z.aml constant at 0.15. - -\code{k.klaus} implements the three wind-based models of Klaus and Vachon (2020), -each fit to 46 globally distributed lakes (their Table 3, lowest-AIC parameterizations). -The \code{method} argument selects the functional form of the k600-wind relationship: -"linear" and "power" scale k600 with wind speed, lake area, and the scale of spatial -integration (\code{spatial.int}); "exp" is exponential in wind speed and additionally depends on -the shoreline development index (\code{sdi}). Note: Klaus and Vachon (2020) caution that -the SDI-based ("exp") model is not mechanistically meaningful and should not be used for -predictive purposes. } \examples{ data.path = system.file('extdata', package="LakeMetabolizer") diff --git a/man/k.read.base.Rd b/man/k.read.base.Rd index 1c595a3..81c9864 100644 --- a/man/k.read.base.Rd +++ b/man/k.read.base.Rd @@ -28,7 +28,7 @@ k.vachon.base(wnd, lake.area, params=c(2.51,1.48,0.39)) k.heiskanen.base(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet) -k.klaus.base(wnd, wnd.z, lake.area, spatial.int, sdi = NULL, method = c("linear", "exp", "power")) +k.klaus.base(wnd, wnd.z, lake.area, spatial.int, method = c("linear", "exp", "power")) } \arguments{ \item{wnd.z}{Height of wind measurement, (Units: m)} @@ -59,10 +59,6 @@ k.klaus.base(wnd, wnd.z, lake.area, spatial.int, sdi = NULL, method = c("linear" \item{method}{Only for \link{k.crusius.base} and \link{k.klaus.base}. String of valid method. For k.crusius.base either "constant", "bilinear", or "power". For k.klaus.base either "linear", "exp", or "power"} -\item{sdi}{Only for \link{k.klaus}'s "exp" method. Numeric shoreline development index -(unitless, >= 1): the ratio of shoreline length to the circumference of a circle of equal -lake area, SDI = P / (2 * sqrt(pi * lake.area)).} - \item{spatial.int}{Only for \link{k.klaus}. Numeric scale of spatial integration, expressed relative to total lake surface area: from near 0 (a point-scale measurement, e.g. 1 m^2 / lake.area) up to but not including 1 (whole-lake integration). Must satisfy 0 < spatial.int < 1; the value 1 is @@ -81,15 +77,6 @@ Returns the gas exchange velocity based on the chosen model in units of m/day Can change default parameters of MacIntyre and Vachon models. Default for Vachon is c(2.51,1.48,0.39). Default for MacIntyre is c(1.2,0.4872,1.4784). Heiskanen et al. (2014) uses MacIntyre model with c(0.5,0.77,0.3) and z.aml constant at 0.15. - -\code{k.klaus} implements the three wind-based models of Klaus and Vachon (2020), -each fit to 46 globally distributed lakes (their Table 3, lowest-AIC parameterizations). -The \code{method} argument selects the functional form of the k600-wind relationship: -"linear" and "power" scale k600 with wind speed, lake area, and the scale of spatial -integration (\code{spatial.int}); "exp" is exponential in wind speed and additionally depends on -the shoreline development index (\code{sdi}). Note: Klaus and Vachon (2020) caution that -the SDI-based ("exp") model is not mechanistically meaningful and should not be used for -predictive purposes } \examples{ wnd.z <- 2 From 4d88791c43d8b3a37d2a4e19bbb89543a5fc0f43 Mon Sep 17 00:00:00 2001 From: Bennett McAfee <105509324+bmcafee@users.noreply.github.com> Date: Wed, 8 Jul 2026 09:54:38 -0400 Subject: [PATCH 20/20] Fix errors in k.klaus documentation --- R/k.read.R | 6 +++--- man/k.read.Rd | 2 +- man/k.read.base.Rd | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/R/k.read.R b/R/k.read.R index c321af7..8fc0371 100644 --- a/R/k.read.R +++ b/R/k.read.R @@ -26,7 +26,7 @@ #' #'k.heiskanen(ts.data, wnd.z, Kd, atm.press) #' -#'k.klaus(ts.data, wnd.z, lake.area, spatial.int, method = c("linear", "exp", "power")) +#'k.klaus(ts.data, wnd.z, lake.area, spatial.int, method = c("linear", "power")) #'@param ts.data vector of datetime in POSIXct format #'@param method Only for \link{k.crusius} and \link{k.klaus}. String of valid method. For k.crusius either "linear", "bilinear", or "power". For k.klaus either "linear" or "power" #'@param wnd.z height of wind measurement, m @@ -214,9 +214,9 @@ k.read = function(ts.data, wnd.z, Kd, atm.press, lat, lake.area){ #' #'k.heiskanen.base(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet) #' -#'k.klaus.base(wnd, wnd.z, lake.area, spatial.int, method = c("linear", "exp", "power")) +#'k.klaus.base(wnd, wnd.z, lake.area, spatial.int, method = c("linear", "power")) #'@param wnd Numeric value of wind speed, (Units:m/s) -#'@param method Only for \link{k.crusius.base} and \link{k.klaus.base}. String of valid method. For k.crusius.base either "constant", "bilinear", or "power". For k.klaus.base either "linear", "exp", or "power" +#'@param method Only for \link{k.crusius.base} and \link{k.klaus.base}. String of valid method. For k.crusius.base either "constant", "bilinear", or "power". For k.klaus.base either "linear" or "power" #'@param wnd.z Height of wind measurement, (Units: m) #'@param Kd Light attenuation coefficient (Units: m^-1) #'@param lat Latitude, degrees north diff --git a/man/k.read.Rd b/man/k.read.Rd index 266e709..8bf9b84 100644 --- a/man/k.read.Rd +++ b/man/k.read.Rd @@ -25,7 +25,7 @@ k.vachon(ts.data, lake.area, params=c(2.51,1.48,0.39)) k.heiskanen(ts.data, wnd.z, Kd, atm.press) -k.klaus(ts.data, wnd.z, lake.area, spatial.int, method = c("linear", "exp", "power")) +k.klaus(ts.data, wnd.z, lake.area, spatial.int, method = c("linear", "power")) } \arguments{ \item{ts.data}{vector of datetime in POSIXct format} diff --git a/man/k.read.base.Rd b/man/k.read.base.Rd index 81c9864..7af7312 100644 --- a/man/k.read.base.Rd +++ b/man/k.read.base.Rd @@ -28,7 +28,7 @@ k.vachon.base(wnd, lake.area, params=c(2.51,1.48,0.39)) k.heiskanen.base(wnd.z, Kd, atm.press, dateTime, Ts, z.aml, airT, wnd, RH, sw, lwnet) -k.klaus.base(wnd, wnd.z, lake.area, spatial.int, method = c("linear", "exp", "power")) +k.klaus.base(wnd, wnd.z, lake.area, spatial.int, method = c("linear", "power")) } \arguments{ \item{wnd.z}{Height of wind measurement, (Units: m)} @@ -57,7 +57,7 @@ k.klaus.base(wnd, wnd.z, lake.area, spatial.int, method = c("linear", "exp", "po \item{lwnet}{Numeric value net long wave radiation, W m^-2} -\item{method}{Only for \link{k.crusius.base} and \link{k.klaus.base}. String of valid method. For k.crusius.base either "constant", "bilinear", or "power". For k.klaus.base either "linear", "exp", or "power"} +\item{method}{Only for \link{k.crusius.base} and \link{k.klaus.base}. String of valid method. For k.crusius.base either "constant", "bilinear", or "power". For k.klaus.base either "linear" or "power"} \item{spatial.int}{Only for \link{k.klaus}. Numeric scale of spatial integration, expressed relative to total lake surface area: from near 0 (a point-scale measurement, e.g. 1 m^2 / lake.area)