sapflux: Analyze Sapflux Traces

#' @title Zero method that uses a threshold vapor-pressure deficit approach.
#'
#' @description This function uses a threshold minimum vapor pressure deficit
#'              that serves as an indicator of probable 'zero-flow' conditions.
#'
#'              Not implemented yet.
#'
#' @param flux       'flux' class object
#' @param time       'time' vector, must be same length as data.
#' @param vpd        Vapor-pressure deficit vector, must be same length as time.
#' @param threshold  Zero-flow VPD threshold, in kPA
#' @param limit      Time limit, in days, in which to find a VPD threshold.
#'
#' @details
#'
#' The time limit serves to limit how far back or forward a dT maximum is looked
#' for. This can solve computation errors as well as limit the effect that
#' mean daily temperature has on dT maximums. For example, a dT max of 40C might
#' be valid in California during the summer, but the same probes at zero flow
#' may max out at 20C in boreal Canada during March.
#'
#' The time limit has units of full days and extends forward and backward of the
#' measurement.
#'
#' @references
#' Kavanagh, K., Pangle, R., D Schotzko, A., 2007.
#'     Nocturnal transpiration causing disequilibrium between soil and stem
#'     predawn water potential in mixed conifer forests of Idaho.
#'     Tree physiology 27, 621–9. https://doi.org/10.1093/treephys/27.4.621
#'
#' @family zero
#' @examples
#' VPDThresholdNTE(data = data, time = time)
VPDThresholdNTE <- function(flux, time, vpd,
                            threshold = 0.01, limit = 15
                            ) {
  stop("Not implemented yet.")
  # Check input validity:
  stopifnot(
    is.vector(data),
    is.vector(vpd),
    class(vpd) == "numeric",
    inherits(time, "POSIXct"),
    length(threshold) == 1,
    class(threshold) == "numeric",
    length(time) == length(data),
    length(time) == length(vpd),
    length(limit) == 1,
    class(limit) == "numeric"
  )
  # Filter data by VPD threshold
  vpd_data <- ifelse(vpd < threshold, data, NA)
  limit <- limit * 86400 # Convert limit to POSIXct interval
  time_secs <- as.numeric(time)
  dT_max <- vector(mode = "numeric", length = length(data))
  for (i in 1:length(vpd_data)) {
    # Find the timing bounds, as an index interval.
    if ((time_secs[i] - limit) < time_secs[1]) {
      dT_start <- 1
    } else {
      dT_start <- which(time_secs == (time_secs[i] - limit))
    }
    if ((time_secs[i] + limit) > time_secs[length(time_secs)]) {
      dT_end <- length(time_secs)
    } else {
      dT_end <- which(time_secs == (time_secs[i] + limit))
    }
    dT_max[i] <- max(vpd_data[dT_start:dT_end], na.rm = TRUE)
  }
  stop("Unfinished!")
  # Drop unreasonable outliers.
  data <- ifelse(data > 1, NA, data)
  return(data)
}