R/window_idx.R

In efetac/PolarMetrics: Functions for analyzing seasonal patterns in environmental variables

#' @name window_idx
#' @title Calculate index values marking truncation endpoints for a
#'   window within a series of cumulative sums
#' @description Calculate index values marking the truncation endpoints
#'   for a window within a series of cumulative sums. In a phenological
#'   context these endpoints define the beginning and end of a
#'   quantitatively specified growing season.
#' @param s A vector of numeric values representing a uniformly
#'   sampled series of cumulative sums. \code{s} can be obtained using
#'   \code{\link{sum_cycle}}
#' @param c A numeric value indicating the number of cycles in series \code{s}.
#' @param cy A positive integer in the range [0,c] indicating the cycle
#'   (year) within which to calculate truncation endpoints.
#' @param lb A numeric value in the range [0,0.5] indicating the
#'   percent of cumulative total to truncate from the first half of
#'   cycle number \code{cy}. For example, 0.15 will result in remove
#'   the interval corresponding to [0\%,15\%] representing a window that
#'   begins after 15\% of cumulative annual total is crossed.
#' @param ub A numeric value in the range (0.5,1] indicating the percent
#'   of cumulative total to truncate from the last half of cycle
#'   \code{cy}. For example, 0.8 will result in removing (80\%,100\%]
#'   of cumulative NDVI representing a window that begins after the 85\%
#'   of cumulative annual total is crossed.
#' @details \code{window_idx} returns a vector of two values indicating
#'   endpoints for a window within series \eqn{S}. Cumulative sums within
#'   these endpoints are used to place threshold crossings in the data
#'   such as growing season start and end points. Since a growing season
#'   can be based on crossing points of cumulative sums. Make \eqn{S}
#'   a series of cumulative sums and \eqn{s} a subset of growing season
#'   values within \eqn{S}. \eqn{s} can be defined on the basis of
#'   thresholds crossings of cumulative sums as follows. For a given
#'   threshold value, \eqn{t}, let \eqn{j(t)} be the index \eqn{j}
#'   for \eqn{S} such that
#'   \deqn{\{i: t < S_i < (t-1)\}}{{i: t < S_i < (t-1)}}
#'   Then \eqn{\min j(t)}{min[j(t)]} is the index of \eqn{S} marking the
#'   first value larger than \eqn{t}, which is also the beginning of the
#'   growing season. \eqn{\max j(t)}{max[j(t)]} marks the ending index.
#' \cr\cr
#' Figure 2. This visualization shows how the subset of values
#'   representing the growing season are selected based on a cumulative
#'   NDVI threshold percentage of 15\%. The growing season is defined
#'   to start after cumulative NDVI reaches 15\% of the total for that
#'   phenology-centered year and ends at 85\%.\cr
#'   \if{html}{\figure{ex2.png}{options: alt="Figure: ex2.png"}}
#'   \if{latex}{\figure{ex2.pdf}{options: width=7cm}}
#' @return Returns a vector of the five indexes (from the input)
#'   marking the following timing metrics: early season (lb), early-mid
#'   season (midpoint of lb and 50% of cumulative total), middle season
#'   (50% of cumulative total), late-mid season (midpoint of 50% of
#'   cumulative total and ub), and late season (ub).
#' @examples
#' dpy <- 365                 # Days/yr
#' c <- 12                    # Num. of years/cycles
#' spc <- 46                  # Number of samples in one cycle (yr)
#' data(mndvi)                # Load data
#' t <- as.vector(mndvi$day)  # Days since January 1, 2000
#' r <- t2rad(t,dpy)          # Transform days of year to radians
#' v <- as.vector(mndvi$wc)   # MODIS NDVI for Willow Creek tower, WI
#' vx <- mean(vec.x(r,v), na.rm=TRUE) # Avg horizontal vector
#' vy <- mean(vec.y(r,v), na.rm=TRUE) # Avg vertical vector
#' rv_ang <- vec_ang(vx,vy)   # Angle of resultant vec (point of max activity)
#' av_ang <- avec_ang(rv_ang)  # Angle marking point of least activity
#' av_idx <- rad2idx(av_ang, spc=spc) # Index (1-spc) marking avg start of yr
#' ann_cum <- sum_cycle(v,av_idx,spc=spc)$cumsum # Accum. vals within each yr
#' # Find seasonal beg, end index for the 2nd yr using 15th pctile of cum NDVI
#' cy <- 2                    # The second yr of data (which is 2001 here)
#' es.idx <- window_idx(ann_cum,c-1,cy,0.15,0.8)[1]  # Idx of marking ES
#' ems.idx <- window_idx(ann_cum,c-1,cy,0.15,0.8)[2] # Idx for EMS
#' ms.idx <- window_idx(ann_cum,c-1,cy,0.15,0.8)[3]  # Idx for MS
#' lms.idx <- window_idx(ann_cum,c-1,cy,0.15,0.8)[4] # Idx for LMS
#' ls.idx <- window_idx(ann_cum,c-1,cy,0.15,0.8)[5]  # Idx for LS
#' es <- t[es.idx]          # Early growing season day of pheno yr
#' ems <- t[ems.idx]        # Early-mid growing season day of pheno yr
#' ms <- t[ms.idx]          # Mid (50th %tile) growing season day
#' lms <- t[lms.idx]        # Late-mid growing season day of pheno yr
#' ls <- t[ls.idx]          # Late growing season day of pheno yr
#' Sintv <- ls-es           # Days in the growing season
#' ann_cum[es.idx:ls.idx]   # Show cumulative NDVI vals for growing season
#' @author Bjorn J. Brooks, Danny C. Lee, William W. Hargrove, Lars Y. Pomara
#' @references Brooks, B.J., Lee, D.C., Desai, A.R., Pomara, L.Y.,
#'   Hargrove, W.W. (2017). Quantifying seasonal patterns in
#'   disparate environmental variables using the PolarMetrics R package.
#' @export

window_idx <- function(s,c,cy,lb,ub) {
  if (length(c) == 1 & length(cy) == 1 & length(lb) == 1 & length(ub) == 1 & lb >= 0 & lb <= 0.5 & ub > 0.5 & ub <= 1) {
    spc = length(s) / c # num of samples per cycle
    if (spc %% 1 != 0) {
      stop('Length of argument s is not evenly divisible by arg c')
    }
    idx.s <- (spc*(cy-1)+1):(spc*cy)               # Idx for this cycle's sbst
    if (sum(is.na(s[idx.s]))>0) {
      es <- NA; ems <- NA; ms <- NA; lms <- NA; ls <- NA
    } else {
      x <- s[idx.s]                                # Subset x
      pr.x <- x / max(x)                           # Convert to percentages
      es <- idx.s[min(which(pr.x > lb))]           # Lower threshold index
      ems <- idx.s[min(which(pr.x > ((lb+0.5)/2)))] # Midpoint of ES, MS index
      ms <- idx.s[min(which(pr.x > 0.5))]          # 50% of cum. total index
      lms <- idx.s[min(which(pr.x > ((ub+0.5)/2)))] # Midpoint of MS, LS index
      ls <- idx.s[min(which(pr.x >= ub))]          # Upper threshold index
    }

    return(c(es, ems, ms, lms, ls))
  } else {
    stop('Arguments lb (ub) should be 0<=lb<=0.5 (0.5<ub<=1)')
  }
}