R/GranierConversions.R

In bmcnellis/sapflux: Analyze Sapflux Traces

#' Performs Granier-relevant data transformation
#'
#' @description
#'
#' Conversion here converts raw/partially converted data to a per-tree flow
#' index, in three steps:
#'
#' 1. Convert raw voltage data to temperature differentials
#'
#' 2. Perform automatic dT zeroing
#'
#' 3. Apply the Granier calibration. Defaults to Granier (1987):
#' u = 0.000119 * K ^ 1.231 . Other calibrations supported.
#'
#'
#' @param flux         'flux' class object
#' @param stop_at      This is a hack that should be ignored.
#' @param zero_method  Method to use for temperature-zeroing routine
#' @param index_eqn    Equation to pass to generate the 'flux index'
#' @param refT         Reference probe temperature, calculate Seebeck coeff
#' @param multiplier   Alternative to refT - simple voltage multiplier
#' @param JIT          TRUE/FALSE, use JIT/cmpfun() from package 'compiler' to
#'                     try to speed up the zeroing function?
#' @param ...          Further arguments.
#'
#' @return
#'
#' Returns the 'flux' object with data represented as a transpiration 'density'
#' at the point of measurement, in units of centimeters of H2O per hour. Also
#' updates metadata, when appropriate, as well as the datatype and output logs.
#'
#' @export
#' @family zero
#' @examples
#' # Custom Granier calibration:
#' eqn <- function(x) {
#'   # 'x' is equivalent to 'K' in Granier (1987)
#'     a <- 0.000224
#'     b <- 1.55
#'     u <- a * (K ^ b)
#'     return(u)
#' }
#' # Added a refT to account for lower temperatures at e.g. boreal study site
#' data <- GranierConversions(flux = flux, zero.method = "Omit.NTE",
#'                            index.eqn = eqn, refT = 7)
GranierConversions <- function(flux, stop_at = "flux_density",
                               multiplier = 25, refT = NULL,
                               zero_method = "omit_NTE",
                               index_eqn = NULL, JIT = FALSE,
                               window = c(2300, 300)) {
  # Input checks
  validObject(flux)
  params <- LoadDefaults(flux)
  stopifnot(
    is.numeric(multiplier), length(multiplier) == 1,
    zero_method %in% c("omit_NTE"),
    length(zero_method) == 1,
    is.null(refT) | all(is.numeric(refT), length(refT) == 1),
    is.null(index_eqn) | inherits(index_eqn, "function"),
    length(window) == 2, is.numeric(window)
  )
  zero_log <- FALSE
  # Pull slots
  data     <- slot(object = flux, name = "data")
  NAprev <- sum(is.na(data), na.rm = TRUE)
  datatype <- slot(object = flux, name = "datatype")
  time     <- slot(object = flux, name = "time")
  #####
  # 1. Raw voltage to dT
  if (datatype == "voltages") {
    message("Converting voltages...")
    #if (!exists("multiplier", frame = topframe, inherits = FALSE)) {
    if (length(refT) < 1) {
    #if (!exists("refT", frame = topframe, inherits = FALSE)) {
      #data <- apply(data, 2, function(x) x * multiplier)
      data <- data * multiplier
    } else {
      coeff <- CalcSeebeck(refT)
      data <- data / coeff
      #data <- apply(data, 2, function(x) x / coeff)
    }
    datatype <- "temperatures"
  }
  #####
  # 2. dT to 'dimensionless flux index'
  if (datatype == "temperatures") {
    message("Zeroing temperatures...")
    # Set defaults
    if (zero_method == "omit_NTE") {
      if (JIT == TRUE) {
        compiler::enableJIT(3)
        zero_compiled <- compiler::cmpfun(ByColumnOmitNTE)
      } else {
        zero_compiled <- ByColumnOmitNTE
      }
      time_lt <- asPOSIXlt(time)
      for (i in 1:ncol(data)) {
        data_out <- tryCatch(expr = {
          data[, i] <- zero_compiled(data = data[, i],
                                       time = time_lt, window = window)
          perc <- 100 * (i / ncol(data))
          cat(round(perc, 2), "% complete.\n")
        }, error = function(cond) {
          print(cond)
          stop(paste("Error in column", i))
        }, finally = {
          invisible()
        })
      }
      if (JIT == TRUE) {
        compiler::enableJIT(0)
      }
    } else {
      stop("Can't zero data - invalid method given")
    }
    datatype <- "flux index"
  }
  #####
  # 3. 'flux index' to E 'density'
  if (datatype == "flux index") {
    message("Scaling index...")
    if (length(index_eqn) < 1) {
      eqn <- GranierEqn
      data <- apply(data, c(1, 2), eqn)
      zero_log <- TRUE
      # Raw returns are in m/s
      #data <- apply(data, c(1, 2), function(x) x * 360000)
      data <- data * 360000
      # This converts to cm/hr
    }
    datatype <- "flux density"
  }
  # Logging, clean-up and return ####
  NApost <- sum(is.na(data))
  NAdiff <- round(abs(NApost - NAprev) / length(data) * 100, 2)
  log <- slot(flux, "log")
  if (zero.method == "Omit_NTE" &&
      zero.log == TRUE) {
    log_message <- paste(
      "Data was zeroed using function 'Omit NTE', which zeroes each 24-hour",
      "window of measurements by the following night's pseudo-temperature",
      "maximum. This minimum is looked for in a window that excludes any",
      "potential nighttime transpiration and may cause a concomittant",
      "underestimation of total fluxes by between 0% and >65%, depending",
      "on species, ecosystem, and environmental variables."
      )
    log <- c(slot(flux, "log"), log_message)
  }
  if ((NApost - NAprev) > 0) {
    log_message <- paste(
      "Granier conversions dropped", NApost - NAprev,
      "data points, corresponding to", NAdiff,
      "% of total data."
    )
    log <- c(slot(flux, "log"), log_message)
  }
  slot(flux, "log") <- log
  slot(flux, "data") <- data
  slot(flux, "datatype") <- datatype
  return(flux)
}