#' @name calc_metrics
#' @title Easy function for calculating all polar meatrics.
#' @description A basic, but fast function for calculating all polar
#' metrics from arrays or xts objects.
#' @importFrom stats sd
#' @importFrom zoo index
#' @param input Either a vector of values or a 1-column xts object
#' indexed with time stamps. If \code{input} is an xts object then omit
#' argument \code{t}. Else if \code{input} is a vector of values then
#' \code{t} should be set to a corresponding vector of time values.
#' @param t An optional vector of time values (e.g., days) corresponding
#' to the input vector sampling points. Do not use this argument if
#' \code{input} is an xts object.
#' @param yr_type Argument specifying either 'cal_yr' for output (of
#' timing variables) given in days starting from Jan. 1, or 'rot_yr'
#' for output in days starting from the average seasonal minimum.
#' @param timing Argument specifying 'from_vecs_integer' or
#' 'from_vecs_dpc' or 'from_percentiles_dpc'. Default is 'from_vecs_integer',
#' which means the timing variable EMS, MS, & LMS will be calculated
#' from the mean of VECTORS within their respective milestones with a
#' precision to the nearest INTEGER (not the original dpc sampling of
#' the data). 'from_vecs_dpc' means the timing variables EMS, MS & LMS
#' will again be calculated from the mean of VECTORS, but only at one
#' of the 1->dpc sampling points. 'from_percentiles_dpc' means all timing
#' variables will be calculated from their percentile completion
#' milestones at one of dpc sampling points. Note that ES & LS
#' variables are always calculated from percentiles as is required by
#' this methodology of defining the range of the growing season.
#' @param spc Integer value specifying the number of samples per cycle
#' (measurements per year) in input.
#' @param lcut Numeric value in the range [0,0.5] passed to
#' \code{window_idx} function. Indicates the percentile to truncate
#' from the first half of each cycle. For example, 0.15 will result
#' in remove the interval corresponding to [0\%,15\%] representing a
#' window that begins after the 15th-percentile is crossed.
#' @param hcut Numeric value in the range (0.5,1] passed to
#' \code{window_idx} function. Indicates the percentile to truncate
#' from the last half of each cycle. For example, 0.85 will result in
#' remove the interval corresponding to [85\%,100\%] representing a
#' window that begins after the 85th-percentile is crossed.
#' @param return_vecs logical argument specifying whether or not to
#' include all of the horizontal and vertical component vectors in output.
#' @param sin_cos logical argument. If TRUE then each timing
#' metric (es, ms, etc.) is returned as its sine and cosine components,
#' es is returned as es_sin and es_cos.
#' @details \code{calc_metrics} runs through the entire polar
#' transformation process and conveniently outputs the final polar
#' metrics for all years included in the input.
#' @return Returns a list with all of the derived polar metrics (e.g.,
#' early season, mid season, etc.). Timing variables in output can
#' be returned relative to the standard calendar year or rotated to
#' a relative year using \code{yr_type} argument.
#' Sine and cosine components can also be returned instead of days since
#' start using \code{sin_cos} argument.
#' The optional argument \code{return_vecs} can be used to add
#' the actual vecturs used in deriving the polar metrics to the returned
#' list. This will also return the overall average (all years) statistics
#' of the resultant vector (rv) and its opposite the anti-vector (av).
#' Table below indicates the variable components of each object within
#' the returned list.
#' metrics (data frame): year,
#' es (or es_sin, es_cos), # Early ssn DOY
#' ems (or ems_sin, ems_cos), # Early-mid ssn DOY
#' ms (or ms_sin, ms_cos), # Mid ssn DOY
#' lms (or lms_sin, lms_cos), # Late-mid ssn DOY
#' ls (or ls_sin, ls_cos), # Late ssn DOY
#' s_intv, # Season length(days)
#' s_avg, # Avg data val in ssn
#' s_sd, # StDev of data dur ssn
#' s_mag, # Mag of avg vec in ssn
#' ems_mag, # Mag early-mid ssn vec
#' lms_mag, # Mag late-mid ssn vec
#' s_mag_std, # s_mag stdzd by GS NDVI
#' a_avg # Avg data val of yr
#' component_vectors (data frame): VX, VY # Hrz, vert vec comp.
#' average_vectors (data frame): rv_idx, rv_ang, # Resultant vector
#' rv_doy, rv_mag, # attributes
#' av_idx, av_ang, av_doy # Anti-vector attr.
#' @examples
#' library(PolarMetrics)
#' library(xts)
#' input <- xts(mndvi$fef, as.Date(mndvi$date))
#' ### Calculate polar measures relative to calendar year
#' calc_metrics(input, yr_type='cal_yr', spc=46, lcut=0.15, hcut=0.8,
#' return_vecs=FALSE, sin_cos=FALSE)
#' ### Calculate as above and return sine, cosine components of timing metrics
#' calc_metrics(input, yr_type='cal_yr', spc=46, lcut=0.15, hcut=0.8,
#' return_vecs=FALSE, sin_cos=TRUE)
#' ### Calculate & return the average vectors for the entire time series
#' calc_metrics(input, yr_type='cal_yr', spc=46, lcut=0.15, hcut=0.8,
#' return_vecs=TRUE, sin_cos=FALSE)$avg_vectors
#' ### Calculate & return the horizontal and vertical vector components
#' head(calc_metrics(input, yr_type='cal_yr', spc=46, lcut=0.15, hcut=0.8,
#' return_vecs=TRUE, sin_cos=FALSE)$vectors)
#' @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
calc_metrics <- function(input, t=NULL, timing='from_vecs_integer', yr_type, spc, lcut, hcut, return_vecs, sin_cos) {
dpy <- 365 # Days per year
if (is.null(t)) { # If no values supplied for t then obtain from xts index
t <- as.integer(format(index(input), '%j')) # Time coordinates
}
t <- as.vector(t)
input <- as.vector(input)
if (length(input) != length(t)) {
stop('input and t should be have the same number of values')
}
if (all(timing != c('from_vecs_integer', 'from_vecs_dpc', 'from_percentiles_dpc'))) {
stop('Acceptable arguments for timing are: from_vecs_integer, from_vecs_dpc, or from_percentiles_dpc')
}
nyr <- length(input)/spc # No. of years in input
r <- t2rad(t, dpc=dpy) # Transform days of year to radians
v <- as.vector(input) # Input values (e.g., NDVI)
VX <- vec.x(r,v) # Horizontal vectors
VY <- vec.y(r,v) # Vertical vectors
vx <- mean(VX, na.rm=T) # Avg horizontal vector
vy <- mean(VY, na.rm=T) # Avg vertical vector
rv_ang <- vec_ang(vx,vy) # Angle of resultant vec (point of max activity)
rv_idx <- rad2idx(rv_ang, spc=spc) # Index (1-spc) marking avg vector
rv_doy <- rad2d(rv_ang, dpc=dpy) # Day of year equivalent of rv_ang
rv_mag <- vec_mag(vx,vy) # Magnitude (length) of resultant vector
av_ang <- avec_ang(rv_ang) # Angle marking point of least activity
av_doy <- rad2d(av_ang, dpc=dpy) # Day of year equivalent of av_ang
av_idx <- rad2idx(av_ang, spc=spc) # Idx (1-spc) marking avg start of yr
ann_cum <- sum_cycle(v,av_idx,spc=spc)$cumsum # Accum vals within rot. yrs
npy <- length(ann_cum)/spc # No. of complete years in output
# Note, output is re-centered, and accumulates starting from av_idx
# (offset of variable yr from calendar yr) and has nyr-1 yrs of data
# due to re-centering.
if (!isTRUE(sin_cos)) { # if false
output <- data.frame(yr=rep(NA,npy),
es=rep(NA,npy), ems=rep(NA,npy),
ms=rep(NA,npy), lms=rep(NA,npy), ls=rep(NA,npy),
s_intv=rep(NA,npy), s_avg=rep(NA,npy),
s_sd=rep(NA,npy), s_mag=rep(NA,npy),
s_mag_std=rep(NA,npy),
ems_mag=rep(NA,npy), lms_mag=rep(NA,npy),
a_avg=rep(NA,npy))
} else if (isTRUE(sin_cos)) { # if true
output <- data.frame(yr=rep(NA,npy),
es=rep(NA,npy), ems=rep(NA,npy),
ms=rep(NA,npy), lms=rep(NA,npy), ls=rep(NA,npy),
s_intv=rep(NA,npy), s_avg=rep(NA,npy),
s_sd=rep(NA,npy), s_mag=rep(NA,npy),
s_mag_std=rep(NA,npy),
ems_mag=rep(NA,npy), lms_mag=rep(NA,npy),
a_avg=rep(NA,npy),
es_sin=rep(NA,npy), es_cos=rep(NA,npy),
ems_sin=rep(NA,npy), ems_cos=rep(NA,npy),
ms_sin=rep(NA,npy), ms_cos=rep(NA,npy),
lms_sin=rep(NA,npy), lms_cos=rep(NA,npy),
ls_sin=rep(NA,npy), ls_cos=rep(NA,npy))
}
for (J in 1:npy) { # Calculate pheno param's for each yr
output$yr[J] <- J
wi <- window_idx(ann_cum,npy,J,lcut,hcut) # early, mid, late-S ann_cum idx
es_idx <- wi[1]
ls_idx <- wi[5]
be_idx <- c((spc*(J-1)+1),(spc*J)) # This cycle's beg/end indices
es2ls <- es_idx:(ls_idx-1) # ann_cum idx from es2ls
# Calculate from VECTOR ANGLES & with DPY (365) divisions
if (timing == 'from_vecs_integer') {
ms_ang <- vec_ang(mean(VX[es2ls], na.rm=TRUE), # Avg GS vector angle
mean(VY[es2ls], na.rm=TRUE))
ms_idx <- which.max(r[es2ls] > ms_ang) + es_idx # Index of MS milestone
es2ms <- es_idx:(ms_idx-1) # ann_cum idx erly to mid
ms2ls <- ms_idx:(ls_idx-1) # ann_cum idx mid to late
# Angle of early-to-mid season average vector
ems_ang <- vec_ang(mean(VX[es2ms], na.rm=TRUE),
mean(VY[es2ms], na.rm=TRUE))
# Angle of mid-to-late season average vector
lms_ang <- vec_ang(mean(VX[ms2ls], na.rm=TRUE),
mean(VY[ms2ls], na.rm=TRUE))
ms <- rad2d(ms_ang, dpc=dpy) # DOY for mean vector between ES & LS
ems <- rad2d(ems_ang, dpc=dpy) # DOY for mean vector between ES & MS
lms <- rad2d(lms_ang, dpc=dpy) # DOY for mean vector between MS & LS
# Calculate from VECTOR ANGLES & with DPC (divisions
# per cycle) divisions. In the case of 8-day MODIS
# NDVI, this means that the output timing dates will
# only fall on one of the 46 8-day time points.
} else if (timing == 'from_vecs_dpc') {
ms_ang <- vec_ang(mean(VX[es2ls], na.rm=TRUE), # Avg GS vector angle
mean(VY[es2ls], na.rm=TRUE))
ms_idx <- which.max(r[es2ls] > ms_ang) + es_idx # Index of MS milestone
es2ms <- es_idx:(ms_idx-1) # ann_cum idx erly to mid
ms2ls <- ms_idx:(ls_idx-1) # ann_cum idx mid to late
# Angle of early-to-mid season average vector
ems_ang <- vec_ang(mean(VX[es2ms], na.rm=TRUE),
mean(VY[es2ms], na.rm=TRUE))
# Angle of mid-to-late season average vector
lms_ang <- vec_ang(mean(VX[ms2ls], na.rm=TRUE),
mean(VY[ms2ls], na.rm=TRUE))
ems_idx <- which.max(r[es2ms] > ems_ang) + es_idx # Idx of EMS mlestne
lms_idx <- which.max(r[ms2ls] > lms_ang) + ms_idx # Idx of LMS mlestne
ems <- t[ems_idx] # DOY for lcut+50/2
ms <- t[ms_idx] # DOY for 50 %tile
lms <- t[lms_idx] # DOY for hcut+50/2
# Else get indices from PERCENTILES and use only
# (dpc divisions per cycle) divisions, such that the
# timing indices can only fall on one of those time
# points. For example if dpc=46, then any timing
# output would only fall on one of those 46 time
# points (spaced roughly 8 days apart) during the
# year, as opposed to any of the 365 days.
} else if (timing == 'from_percentiles_dpc') {
ems_idx <- wi[2]
ms_idx <- wi[3]
lms_idx <- wi[4]
es2ms <- es_idx:(ms_idx-1) # ann_cum idx erly to mid
ms2ls <- ms_idx:(ls_idx-1) # ann_cum idx mid to late
ems <- t[ems_idx] # DOY for lcut+50/2
ms <- t[ms_idx] # DOY for 50 %tile
lms <- t[lms_idx] # DOY for hcut+50/2
}
es <- t[es_idx] # DOY lcut, 15 %tile
ls <- t[ls_idx] # DOY for hcut
if (yr_type == 'cal_yr') {
# Timing variables relative to the calendar year
output$es[J] <- as.integer(ceiling((es + # DOY for ES milestone
av_doy) %% dpy))
output$ems[J] <- as.integer(ceiling((ems + # DOY for EMS milestone
av_doy) %% dpy))
output$ms[J] <- as.integer(ceiling((ms + # DOY for MS milestone
av_doy) %% dpy))
output$lms[J] <- as.integer(ceiling((lms + # DOY for MLS milestone
av_doy) %% dpy))
output$ls[J] <- as.integer(ceiling((ls + # DOY for LS milestone
av_doy) %% dpy))
if (isTRUE(sin_cos)) {
# Sine, cosine timing variables relative to the calendar year
output$es_sin[J] <- sin(r[es_idx] + av_ang) # sin rad.ang. for ES
output$es_cos[J] <- cos(r[es_idx] + av_ang) # cos rad.ang. for ES
if (all(timing == 'from_vecs_integer')) {
output$ems_sin[J] <- sin(ems_ang + av_ang) # sin rad.ang. for EMS
output$ems_cos[J] <- cos(ems_ang + av_ang) # cos rad.ang. for EMS
output$ms_sin[J] <- sin(ms_ang + av_ang) # sin rad.ang. for MS
output$ms_cos[J] <- cos(ms_ang + av_ang) # cos rad.ang. for MS
output$lms_sin[J] <- sin(lms_ang + av_ang) # sin rad.ang. for LMS
output$lms_cos[J] <- cos(lms_ang + av_ang) # cos rad.ang. for LMS
} else {
output$ems_sin[J] <- sin(r[ems_idx] + av_ang) # sin rad.ang. for EMS
output$ems_cos[J] <- cos(r[ems_idx] + av_ang) # cos rad.ang. for EMS
output$ms_sin[J] <- sin(r[ms_idx] + av_ang) # sin rad.ang. for MS
output$ms_cos[J] <- cos(r[ms_idx] + av_ang) # cos rad.ang. for MS
output$lms_sin[J] <- sin(r[lms_idx] + av_ang) # sin rad.ang. for LMS
output$lms_cos[J] <- cos(r[lms_idx] + av_ang) # cos rad.ang. for LMS
}
output$ls_sin[J] <- sin(r[ls_idx] + av_ang) # sin rad.ang. for LS
output$ls_cos[J] <- cos(r[ls_idx] + av_ang) # cos rad.ang. for LS
}
} else if (yr_type == 'rot_yr') {
# Timing variables relative to the rotated year
output$es[J] <- es %% dpy # DOY for ES milestone
output$ems[J] <- ems %% dpy # DOY for EMS milestone
output$ms[J] <- ms %% dpy # DOY for MS milestone
output$lms[J] <- lms %% dpy # DOY for MLS milestone
output$ls[J] <- ls %% dpy # DOY for LS milestone
if (isTRUE(sin_cos)) {
# Sine, cosine timing vars relative to the rotated yr
output$es_sin[J] <- sin(r[es_idx]) # sin rad.ang. for ES
output$es_cos[J] <- cos(r[es_idx]) # cos rad.ang. for ES
if (all(timing == 'from_vecs_integer')) {
output$ems_sin[J] <- sin(ems_ang) # sin rad.ang. for EMS
output$ems_cos[J] <- cos(ems_ang) # cos rad.ang. for EMS
output$ms_sin[J] <- sin(ms_ang) # sin rad.ang. for MS
output$ms_cos[J] <- cos(ms_ang) # cos rad.ang. for MS
output$lms_sin[J] <- sin(lms_ang) # sin rad.ang. for LMS
output$lms_cos[J] <- cos(lms_ang) # cos rad.ang. for LMS
} else {
output$ems_sin[J] <- sin(r[ems_idx]) # sin rad.ang. for EMS
output$ems_cos[J] <- cos(r[ems_idx]) # cos rad.ang. for EMS
output$ms_sin[J] <- sin(r[ms_idx]) # sin rad.ang. for MS
output$ms_cos[J] <- cos(r[ms_idx]) # cos rad.ang. for MS
output$lms_sin[J] <- sin(r[lms_idx]) # sin rad.ang. for LMS
output$lms_cos[J] <- cos(r[lms_idx]) # cos rad.ang. for LMS
}
output$ls_sin[J] <- sin(r[ls_idx]) # sin rad.ang. for LS
output$ls_cos[J] <- cos(r[ls_idx]) # cos rad.ang. for LS
}
}
output$s_intv[J] <- ls - es # Days in the grw season
v_mu <- mean(v[es2ls], na.rm=TRUE) # Mean for Seasn.
output$s_avg[J] <- v_mu
output$s_sd[J] <- sd(v[es2ls]) # Std. dev. for Seasn.
output$a_avg[J] <- mean(v[be_idx[1]:be_idx[2]], na.rm=TRUE) # Full yr mean
output$s_mag[J] <- vec_mag(mean(VX[es2ls], na.rm=TRUE), # Avg GS vector length
mean(VY[es2ls], na.rm=TRUE))
# s_mag standardized by mean NDVI during the growing season
output$s_mag_std[J] <- output$s_mag[J] /
mean(v[es2ls], na.rm=TRUE)
# Magnitude of avg vec between ES & MS thresholds
output$ems_mag[J] <- vec_mag(mean(VX[es2ms], na.rm=TRUE),
mean(VY[es2ms], na.rm=TRUE))
# Magnitude of avg vec between MS & LS thresholds
output$lms_mag[J] <- vec_mag(mean(VX[ms2ls], na.rm=TRUE),
mean(VY[ms2ls], na.rm=TRUE))
}
if (!isTRUE(return_vecs)) { # if false
return(output)
} else if (isTRUE(return_vecs)) { # if true
# If vector components requested then return by combining as list object
output <- list('metrics' <- output,
'component_vectors' <- data.frame(VX=VX, VY=VY),
'average_vectors' <- data.frame(rv_idx=rv_idx,
rv_ang=rv_ang, rv_doy=rv_doy,
rv_mag=rv_mag, av_idx=av_idx,
av_ang=av_ang, av_doy=av_doy))
return(output)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.