Nothing
# colors <- c("blue", "green3", "orange", "red")
colors <- c("green3", "darkgreen", "darkorange3", "red")
linewidth <- 1.2
# identify greenup or dormancy(brown) period
.Greenup <- function(x, ...) {
ratio.deriv <- c(NA, diff(x))
greenup <- rep(NA, length(x))
greenup[ratio.deriv > 0] <- TRUE
greenup[ratio.deriv < 0] <- FALSE
return(greenup)
}
#' Phenology extraction in Threshold method (TRS)
#'
#' @param x numeric vector, or `fFIT` object returned by [curvefit()].
#' @param t `doy` vector, corresponding doy of vegetation index.
#' @param approach to be used to calculate phenology metrics.
#' 'White' (White et al. 1997) or 'Trs' for simple threshold.
#' @param trs threshold to be used for approach "Trs", in (0, 1).
#' @param IsPlot whether to plot?
#' @param ... other parameters to PhenoPlot
#' @param asymmetric If true, background value in spring season and autumn season
#' is regarded as different.
#'
#' @example R/examples/ex-PhenoTrs.R
#' @seealso [PhenoDeriv()], [PhenoGu()], [PhenoKl()]
#' @export
PhenoTrs <- function(x, t = NULL,
approach = c("White", "Trs"), trs = 0.5, #, min.mean = 0.1
asymmetric = TRUE, IsPlot = TRUE, ...)
{
UseMethod("PhenoTrs", x)
}
#' @rdname PhenoTrs
#' @export
PhenoTrs.fFIT <- function(x, t = NULL, ...) {
if (!is.null(x$tout)) t <- x$tout
values <- last2(x$zs)
PhenoTrs.default(values, t, ...)
}
#' @export
#' @rdname PhenoTrs
PhenoTrs.default <- function(x, t = NULL,
approach = c("White", "Trs"), trs = 0.5, # , min.mean = 0.1
asymmetric = TRUE,
IsPlot = TRUE, ...) {
metrics <- c(sos = NA, eos = NA)
n <- length(x)
# get peak of season position
half.season <- median(which.max(x)) %>% round() # + 20, half season + 20 was unreasonable
pos <- t[half.season]
if (all(is.na(x))) {
return(metrics)
}
if (half.season < 5 || half.season > (n - 5)) {
return(metrics)
}
# get statistical x
# avg <- mean(x, na.rm = TRUE)
# avg2 <- mean(x2[x2 > min.mean], na.rm = TRUE)
peak <- max(x, na.rm = TRUE)
if (asymmetric) {
mn_a <- min(x[1:half.season], na.rm = T)
mn_b <- min(x[-(1:half.season)], na.rm = T)
mn <- c(
rep(mn_a, half.season),
rep(mn_b, n - half.season)
)
} else {
mn <- min(x, na.rm = TRUE)
}
ampl <- peak - mn
# select (or scale) x and thresholds for different methods
approach <- approach[1]
if (approach == "White") {
# scale annual time series to 0-1
ratio <- (x - mn) / ampl
# trs <- 0.5
trs.low <- trs - 0.05
trs.up <- trs + 0.05
}
if (approach == "Trs") {
ratio <- x
a <- diff(range(ratio, na.rm = TRUE)) * 0.05
trs.low <- trs - a
trs.up <- trs + a
}
greenup <- .Greenup(ratio)
## distinguish the first half year and second year
# select time where SOS and EOS are located (around trs value)
bool <- ratio >= trs.low & ratio <= trs.up
# get SOS, EOS, LOS
# fixed 2017-01-04, according to TP phenology property
# sos <- round(median(sose os[greenup & bool], na.rm = TRUE))
# eos <- round(median(soseos[!greenup & bool], na.rm = TRUE))
sos <- round(median(t[greenup & bool & t < pos], na.rm = TRUE))
eos <- round(median(t[!greenup & bool & t > pos], na.rm = TRUE))
# greenup <- .Greenup(ratio)
if (is.na(sos + eos)) {
bool2 <- ratio >= (trs - 0.1) & ratio <= (trs + 0.1)
if (is.na(sos)) sos <- round(median(t[greenup & bool2 & t < pos], na.rm = TRUE))
if (is.na(eos)) eos <- round(median(t[!greenup & bool2 & t > pos], na.rm = TRUE))
# if still is NA
if (is.na(sos)) sos <- round((t[1] + pos) / 2)
if (is.na(eos)) eos <- round((t[n] + pos) / 2)
}
# monotonous = length(unique(greenup[2:n]) %>% rm_empty()) == 1
# los <- eos - sos#los < 0 indicate that error
# los[los < 0] <- n + (eos[los < 0] - sos[los < 0])
# get MGS, MSP, MAU
# mgs <- mean(x[ratio > trs], na.rm = TRUE)
# msp <- mau <- NA
# if (!is.na(sos)) {
# id <- (sos - 10):(sos + 10)
# id <- id[(id > 0) & (id < n)]
# msp <- mean(x[which(index(x) %in% id == TRUE)], na.rm = TRUE)
# }
# if (!is.na(eos)) {
# id <- (eos - 10):(eos + 10)
# id <- id[(id > 0) & (id < n)]
# mau <- mean(x[which(index(x) %in% id == TRUE)], na.rm = TRUE)
# }
# metrics <- c(sos = sos, eos = eos, los = los, pos = pos, mgs = mgs,
# rsp = NA, rau = NA, peak = peak, msp = msp, mau = mau)
metrics <- c(sos = sos, eos = eos) # , los = los
if (IsPlot) {
main <- ifelse(all(par("mar") == 0), "", sprintf("TRS%d", trs * 10))
PhenoPlot(t, x, main = main, ...)
lines(t, trs * ampl + mn, lwd = linewidth)
# lines(t, trs.low*ampl + mn, lty = 2, lwd = linewidth)
# lines(t, trs.up*ampl + mn, lty = 2, lwd = linewidth)
abline(v = metrics, col = colors[c(1, 4)], lwd = linewidth)
text(metrics[1] - 5, min(trs + 0.15, 1) * ampl[1] + mn[1], "SOS", col = colors[1], adj = c(1, 0))
text(metrics[2] + 5, min(trs + 0.15, 1) * last(ampl) + last(mn), "EOS", col = colors[4], adj = c(0, 0))
}
return(metrics)
### The function returns a vector with SOS, EOS, LOS, POS, MGS, rsp, rau, PEAK, MSP and MAU. }
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.