R/ByColumnOmitNTE.R

In bmcnellis/sapflux: Analyze Sapflux Traces

#' @title Zero method that deliberately excludes NTE.
#'
#' @description Finds a local maximum for each night in the dataset
#'              and zeroes by that - the default window avoids anomalous
#'              local maximums that are *probably* the true maximums, but
#'              can complicate data interpretation.
#'
#' @param flux    'flux' class object
#' @param time    'time' vector, should be same length as data
#' @param window  Window to search for nightly maximums
#'
#'
#' @family zero
#' @examples
#' ByColumnOmitNTE(data = data, days = time)
ByColumnOmitNTE <- function(data, time, window = c(2300, 300)) {
  # Check input validity:
  stopifnot(
    is.vector(data),
    inherits(time, "POSIXlt"),
    length(window) == 2,
    is.numeric(window),
    length(time) == length(data)
  )
  ##### Create 'days' object
  # Window before midnight, as a difftime object:
  diffprev <- as.difftime(tim = 24 - window[1] / 100, units = "hours" )
  diffpost <- as.difftime(tim = window[2] / 100, units = "hours")
  diff24hr <- as.difftime(tim = 24, units = "hours")
  # Find unique days in 'time' vector
  days <- time
  lubridate::hour(days) <- 0
  lubridate::minute(days) <- 0
  lubridate::second(days) <- 0
  days <- unique(days) # 'days' is a vector of type integer
  day_length <- length(days)
  # Fix the 'beginning night' issue:
  ##### The 'days' vector has to be updated with information
  # from the current 'i' loop
  # Check to see if data is totally missing for this midnight:
  for (j in 1:day_length) {
    midnight <- days[j]
    night.sub <- time[which(time > (days[j] - diffprev))]
    night.sub <- night.sub[which(night.sub < (days[j] + diffpost))]
    night.sub <- data[which(time %in% night.sub)]
    if (all(is.na(night.sub))) {
      days[j] <- NA # Preallocated
    }
  }
  # Adenostoma column 1 time, to this from last note: 2.70 minutes
  ##### Midnight processing - this is the bulk of processing
  # time, so the only non-NA days should have real data
  night.count <- 0
  dtm <- vector(mode = "numeric", length = day_length)
  for (j in 1:day_length) {
    night.count <- night.count + 1
    if (all(is.na(days[j:day_length]))) {
      break
    }
    midnight <- days[j]
    if (is.na(midnight)) {
      next
    }
    # Find all four initial boundaries:
    nbeg <- which(time == (midnight - diffprev))
    nend <- which(time == (midnight + diffpost))
    if ((night.count == day_length) &
        length(midnight + diff24hr) < 1) {
      dbeg <- which(time == (midnight + diffpost + diff24hr))
    } else {
      dbeg <- which(time == (midnight + diffpost - diff24hr))
    }
    # dend is identical to nend
    # If any of the boundaries are undefined, look for them:
    rmdr.tstep <- difftime(time[2], time[1])
    rmdr.count <- 0
    while(any(length(dbeg) < 1,
              length(nbeg) < 1,
              length(nend) < 1)) {
      rmdr.count <- rmdr.count + 1
      if (length(dbeg) < 1) {
        dbeg <- which(time == (midnight + diffpost - diff24hr
                               + (rmdr.count * rmdr.tstep)))
      }
      if (length(nbeg) < 1) {
        nbeg <- which(time == (midnight - diffprev +
                                 (rmdr.tstep * rmdr.count)))
      }
      if (length(nend) < 1) {
        nend <- which(time == (midnight + diffpost -
                                 (rmdr.tstep * rmdr.count)))
      }
    }
    if (all(is.na(data[dbeg:nend]))) {
      next
    }
    rmdr.count <- 0
    repeat {
      if (all(is.na(data[nbeg]),
              nbeg < nend)) {
        rmdr.count <- rmdr.count + 1
        nbeg <- which(time == (midnight - diffprev +
                                 (rmdr.tstep * rmdr.count)))
      } else {
        break
      }
    }
    rmdr.count <- 0
    repeat {
      if (all(is.na(data[nend]),
              nend > nbeg)) {
        rmdr.count <- rmdr.count + 1
        nend <- which(time == (midnight + diffpost -
                                 (rmdr.tstep * rmdr.count)))
      } else {
        break
      }
    }
    rmdr.count <- 0
    repeat {
      if (all(is.na(data[dbeg]),
              dbeg < nbeg)) {
        rmdr.count <- rmdr.count + 1
        dbeg <- which(time == (midnight + diffpost - diff24hr
                               + (rmdr.count * rmdr.tstep)))
      } else {
        break
      }
    }
    if (any(
      dbeg == nbeg,
      dbeg == nend,
      nbeg == nend
    )) {
      n <- 0
      while (is.na(dtm[j])) {
        n <- n + 1
        if (n == j) {
          if (dbeg != nend) {
            dtm[j] <- max(data[dbeg:nend])
          } else {
            stop("Need 1st nights data at a minimum")
          }
        } else {
          dtm[j] <- dtm[j - n]
        }
      }
    } else {
      dtm[j] <- max(data[nbeg:nend], na.rm = TRUE)
    }
    # With appropriate boundaries, find DTM:
    # This is the nightly MAXIMUM temperature differential
    # Now Calculate Granier's K:
    if ((night.count == day_length) &
        (length(midnight + diff24hr) < 1)) {
      # For the tail-end remainder
      data[nend:dbeg] <- (dtm[j] - data[nend:dbeg]) / data[nend:dbeg]
    } else {
      # For everything else
      data[dbeg:nend] <- ((dtm[j] - data[dbeg:nend]) / data[dbeg:nend])
    }
  } # end 'j' loop
  # Adenostoma column 1 time, to this from last note: 2.77 minutes
  # Drop unreasonable outliers
  data <- ifelse(data > 1, NA, data)
  return(data)
}