Add ppm to umol conversion utility; update vignette to use it (#29)

NAMESPACE

@@ -4,6 +4,7 @@ export(wtf_compute_fluxes)
 export(wtf_fit_models)
 export(wtf_metadata_match)
 export(wtf_normalize_time)
+export(wtf_ppm_to_umol)
 export(wtf_read_LGR915)
 export(wtf_read_LI7810)
 export(wtf_read_LI7820)

R/models.R

@@ -3,7 +3,10 @@
 #' Fit various models to gas concentration data
 #'
 #' @param time Relative time of observation (typically seconds), numeric
-#' @param conc Greenhouse gas concentration, numeric
+#' @param conc Greenhouse gas concentration (typically ppm or ppb), numeric
+#' @param area Area covered by the measurement chamber (typically cm2), numeric
+#' @param volume Volume of the system
+#' (chamber + tubing + analyzer, typically cm3), numeric
 #' @return A wide-form \code{\link{data.frame}} with fit statistics for linear,
 #' robust linear, and polynomial models. By default, extensive details are
 #' provided only for the linear fi; for robust linear and polynomial, only
@@ -24,7 +27,7 @@
 #' dat$SECONDS <- dat$SECONDS - min(dat$SECONDS) # normalize time to start at 0
 #' plot(dat$SECONDS, dat$CO2)
 #' wtf_fit_models(dat$SECONDS, dat$CO2)
-wtf_fit_models <- function(time, conc) {
+wtf_fit_models <- function(time, conc, area, volume) {
   # Basic linear model
   try(mod <- lm(conc ~ time))
   if(!exists("mod")) {
@@ -157,3 +160,38 @@ wtf_compute_fluxes <- function(data,
   onleft <- c(group_column, time_column)
   return(z[c(onleft, setdiff(names(z), onleft))])
 }
+
+
+#' Convert ppm to micromoles using the Ideal Gas Law
+#'
+#' @param ppm Gas concentration (ppmv), numeric
+#' @param volume System volume (chamber + tubing + analyzer, m3), numeric
+#' @param temp Optional chamber temperature (degrees C), numeric
+#' @param atm Optional atmospheric pressure (Pa), numeric
+#' @return The value(s) in micromoles.
+#' @export
+#' @references Steduto et al.:
+#' Automated closed-system canopy-chamber for continuous field-crop monitoring
+#' of CO2 and H2O fluxes, Agric. For. Meteorol., 111:171-186, 2002.
+#' \url{http://dx.doi.org/10.1016/S0168-1923(02)00023-0}
+#' @note If \code{temp} and/or \code{atm} are not provided, the defaults
+#' are NIST normal temperature and pressure.
+wtf_ppm_to_umol <- function(ppm, volume, temp, atm) {
+
+  if(missing(temp)) {
+    temp <- 20
+    wtf_message("Assuming temp = ", temp, " C")
+  }
+  if(missing(atm)) {
+    atm <- 101325
+    wtf_message("Assuming atm = ", atm, " Pa")
+  }
+
+  # Gas constant, from https://en.wikipedia.org/wiki/Gas_constant
+  R <- 8.31446261815324	 # m3 Pa K−1 mol−1
+  wtf_message("Using R = ", R, " m3 Pa K-1 mol-1")
+  TEMP_K <- temp + 273.15
+
+  # Use ideal gas law to calculate micromoles: n = pV/RT
+  return(ppm * atm * volume / (R * TEMP_K))
+}

tests/testthat/test-models.R

@@ -56,3 +56,21 @@ test_that("wtf_compute_fluxes works", {
   expect_identical(out$time_min, rep(min(times), nrow(out))) # min of raw times
   expect_identical(out$time_max, rep(max(times), nrow(out))) # max of raw times
 })
+
+test_that("wtf_ppm_to_umol works", {
+  withr::local_options(whattheflux.quiet = TRUE)
+
+  # Fixed-value test at STP, assuming our ideal gas law implementation is correct
+  # 0.18 is the slope (ppm CO2/s) of the TG10-01087 example data
+  x <- wtf_ppm_to_umol(0.18, volume = 0.1)
+  expect_equal(round(x, 4), 0.7483)
+
+  # Colder temperature means larger flux for a given concentration rise
+  expect_gt(wtf_ppm_to_umol(0.18, volume = 0.1, temp = 5), x)
+  # Lower pressure means smaller flux for a given concentration rise
+  expect_lt(wtf_ppm_to_umol(0.18, volume = 0.1, atm = 75000), x)
+  # Larger volume means larger flux for a given concentration rise
+  expect_gt(wtf_ppm_to_umol(0.18, volume = 1), x)
+
+  # Not sure what else to test
+})

vignettes/intro-to-whattheflux.Rmd

@@ -112,6 +112,26 @@ p %+% dat
 
 That looks better!
 
+## Unit conversion
+
+We'd like our final units to be in µmol/m2/s, and so need to do some unit
+conversion. (This can happen either before or after flux computation, below.)
+The package provides `wtf_ppm_to_umol()` that does this conversion using the
+[Ideal Gas Law](https://en.wikipedia.org/wiki/Ideal_gas_law).
+
+```{r units}
+dat$CO2_umol <- wtf_ppm_to_umol(dat$CO2, 
+                            volume = 0.1, # m3
+                            temp = 24)    # degrees C
+
+# Note that because we didn't provide the 'atm' parameter, 
+# wtf_ppm_to_umol assumed standard pressure.
+
+# Also normalize by ground area (0.16 m2)
+dat$CO2_umol_m2 <- dat$CO2_umol / 0.16
+```
+
+
 ## Compute fluxes
 
 The `wtf_compute_fluxes` function provides a general-purpose tool for
@@ -122,7 +142,7 @@ QA/QC information.
 fluxes <- wtf_compute_fluxes(dat,
                            group_column = "Plot", 
                            time_column = "TIMESTAMP", 
-                           conc_column = "CO2",
+                           conc_column = "CO2_umol_m2",
                            dead_band = 10)
 
 # By default, wtf_compute_fluxes returns a data.frame with one row per
@@ -141,7 +161,8 @@ fluxes[c(1:2, 4, 14, 15, 20)]
 ggplot(fluxes, aes(Plot, slope_estimate, color = adj.r.squared)) +
   geom_point() +
   geom_linerange(aes(ymin = slope_estimate - slope_std.error,
-                     ymax = slope_estimate + slope_std.error))
+                     ymax = slope_estimate + slope_std.error)) +
+  ylab("CO2 flux (µmol/m2/s)")
 ```
 
 We might want to check whether the robust-linear slope diverges from the