R/extract_pheno.R

In ranghetti/sen2rts: Build and Analyse Sentinel-2 Time Series

#' @title Extract phenological metrics
#' @description Extract phenological metrics from a fitted time series.
#' @param data List of fitted time series as generated by function `fit_curve()`.
#' @param method Thresholding method among `"trs"`, `"derivatives"`, 
#'  `"klosterman"` and `"gu"` (see `phenopix::PhenoExtract()`).
#' @param trs Argument passed to `phenopix::PhenoExtract()` if `method = "trs"`.
#'  It is valid also for `method = "derivatives"` thanks to a `{sen2rts}`
#'  implementation (it will be documented soon).
#'  Other methods do not make use of this argument.
#' @param ... Additional arguments passed to `phenopix::PhenoExtract()`.
#' @return A data table with the following fields:
#'  - `id`: the time series ID (see `s2ts`);
#'  - `year`: the year (integer);
#'  - `cycle`: the cycle ID (integer);
#'  - `begin`, `end`, `maxval`: the dates of begin, end and maximum value in
#'      the cycle as computed by `cut_cycles()`;
#'  - other phenology metrics, depending on `method`:
#'      - if `method = "trs"` or `"derivatives"`: `sos`, `eos`, `los`, `pop`,
#'          `mgs`, `rsp`, `rau`, `peak`, `msp`, `mau` (see
#'          `phenopix::PhenoTrs()` or `phenopix::PhenoDeriv()`);
#'      - if `method = "gu"`: `UD`, `SD`, `DD`, `RD`, `maxline`, `baseline`,
#'          `prr`, `psr`, `plateau.slope` (see `phenopix::PhenoGu()`);
#'      - if `method = "klosterman"`: `Greenup`, `Maturity`, `Senescence`,
#'          `Dormancy` (see `phenopix::PhenoKl()`).
#' @author Luigi Ranghetti, PhD (2020) \email{luigi@@ranghetti.info}
#' @import data.table
#' @export
#' @examples
#' # Load input data
#' data("cf")
#' data("ts_filled") # used for plots
#' 
#' # Default extraction ("trs" method with 50% threshold)
#' dt_pheno <- extract_pheno(cf)
#' plot(ts_filled, fitted = cf, pheno = dt_pheno, plot_dates = TRUE)
#' 
#' \donttest{
#' # Customize parameters (e.g. "derivatives" method with 30% threshold)
#' dt_pheno_2 <- extract_pheno(cf, method = "derivatives", trs = 0.3)
#' plot(ts_filled, fitted = cf, pheno = dt_pheno_2, plot_dates = TRUE)
#' 
#' # Other methods: Klosterman
#' dt_pheno_kl <- extract_pheno(cf, method = "klosterman")
#' plot(ts_filled, fitted = cf, pheno = dt_pheno_kl, plot_dates = TRUE)
#' 
#' # Other methods: Gu
#' dt_pheno_gu <- extract_pheno(cf, method = "gu")
#' plot(ts_filled, fitted = cf, pheno = dt_pheno_gu, plot_dates = TRUE)
#' }


extract_pheno <- function(
  data, 
  method = "trs",
  trs = 0.5, ...
) {
  
  # Avoid check notes for data.table related variables
  sos <- begin <- eos <- pop <- mgs <- peak <- msp <- mau <- los <- 
    Greenup <- Maturity <- Senescence <- Dormancy <- 
    UD <- SD <- DD <- maxline <- baseline <- NULL
  
  ## Check arguments
  # TODO
  
  # Check threshold
  if (!method %in% c("trs", "derivatives", "klosterman", "gu")) {
    print_message(
      type = "error",
      "Argument 'method' only accepts values \"trs\", \"derivatives\", ",
      "\"klosterman\" and \"gu\"."
    )
  }
  
  # # Define which metrics are quantitative (y) and which refer to dates (x)
  # metrics_y <- c(
  #   "mgs", "rsp", "rau", "peak", "msp", "mau", # trs / derivatives
  #   "maxline", "baseline", "prr", "psr", "plateau.slope" # gu
  # )
  # metrics_x <- c(
  #   "sos", "eos", "los", "pop", # trs / derivatives
  #   "UD", "SD", "DD", "RD", # gu
  #   "Greenup", "Maturity", "Senescence", "Dormancy" # klosterman
  # )
  
  # Fit for each ID/cycle
  pheno_list <- list()
  for (sel_id in names(data)) {
    pheno_list[[sel_id]] <- list()
    for (sel_year in names(data[[sel_id]])) {
    for (sel_cycle in names(data[[sel_id]][[sel_year]])) {
      sel_metrics <- PhenoExtract(
        data[[sel_id]][[sel_year]][[sel_cycle]], 
        method = method,
        uncert = FALSE, plot = FALSE,
        sf = c(0,1),
        trs = trs, ...
      )$metrics
      pheno_list[[sel_id]][[sel_year]][[sel_cycle]] <- data.table(
        "id" = sel_id,
        "year" = sel_year,
        "cycle" = sel_cycle,
        "begin" = min(data[[sel_id]][[sel_year]][[sel_cycle]]$ts$date),
        "end" = max(data[[sel_id]][[sel_year]][[sel_cycle]]$ts$date),
        "maxval" = data[[sel_id]][[sel_year]][[sel_cycle]]$maxval,
        "weight" = data[[sel_id]][[sel_year]][[sel_cycle]]$weight,
        as.data.frame(as.list(sel_metrics))
      )
      # sel_metrics_x <- sel_metrics[names(sel_metrics) %in% metrics_x]
      # sel_metrics_y <- sel_metrics[names(sel_metrics) %in% metrics_y]
      # sel_metrics[names(sel_metrics) %in% metrics_x]
      # if (!all(is.na(sel_metrics_x))) {
      #   # dates of the metrics: if metrics refer to dates in the cycle range,
      #   # use them (rounded) as indices of the vector of dates (safer method);
      #   # otherwise, add them to the first date (allowing catching the "outliers" dates)
      #   # (the two methods should be equal, since input "ts" must be a daily-filled s2ts).
      #   pheno_x_dates <- if (
      #     any(!round(sel_metrics_x) %in% seq_along(data[[sel_id]][[sel_year]][[sel_cycle]]$ts$date))
      #   ) {
      #     round(sel_metrics_x)-1 + data[[sel_id]][[sel_year]][[sel_cycle]]$ts$date[1]
      #   } else {
      #     data[[sel_id]][[sel_year]][[sel_cycle]]$ts$date[round(sel_metrics_x)]
      #   }
      #   pheno_x_list[[sel_id]][[sel_cycle]] <- data.table(
      #     "id" = sel_id,
      #     "cycle" = sel_cycle,
      #     "date" = c(
      #       range(data[[sel_id]][[sel_year]][[sel_cycle]]$ts$date), # begin - end of the cycle
      #       data[[sel_id]][[sel_year]][[sel_cycle]]$maxval, # date of the peak
      #       pheno_x_dates # dates of the metrics
      #     ),
      #     "pheno" = c(
      #       "begin", "end", 
      #       if (!is.null(data[[sel_id]][[sel_year]][[sel_cycle]]$maxval)) "maxval",
      #       names(sel_metrics_x)
      #     )
      #   )
      # }
      # if (!all(is.na(sel_metrics_y))) {
      #   pheno_y_list[[sel_id]][[sel_cycle]] <- data.table(
      #     "id" = sel_id,
      #     "cycle" = sel_cycle,
      #     "value" = attr(data, "info")$rescale[1] + 
      #       sel_metrics_y * attr(data, "info")$rescale[2], # rescaled metrics
      #     "pheno" = names(sel_metrics_y)
      #   )
      # }
    } # end of sel_cycle FOR cycle
    } # end of sel_year FOR cycle
  } # end of id FOR cycle

  pheno_dt <- rbindlist(lapply(lapply(pheno_list, lapply, rbindlist), rbindlist))
  
  # Format metrics basing on method
  if (method %in% c("trs", "derivatives") ) {
    pheno_dt[,c("sos", "eos", "pop", "mgs", "peak", "msp", "mau") := list(
      as.Date(sos, origin = begin)-1,
      as.Date(eos, origin = begin)-1,
      as.Date(pop, origin = begin)-1,
      attr(data, "info")$rescale[1] + mgs * attr(data, "info")$rescale[2],
      attr(data, "info")$rescale[1] + peak * attr(data, "info")$rescale[2],
      attr(data, "info")$rescale[1] + msp * attr(data, "info")$rescale[2],
      attr(data, "info")$rescale[1] + mau * attr(data, "info")$rescale[2]
    )]
    pheno_dt[,los := as.integer(los)]
  } else if (method %in% c("klosterman") ) {
    pheno_dt[,c("Greenup", "Maturity", "Senescence", "Dormancy") := list(
      as.Date(Greenup, origin = begin)-1,
      as.Date(Maturity, origin = begin)-1,
      as.Date(Senescence, origin = begin)-1,
      as.Date(Dormancy, origin = begin)-1
    )]
  } else if (method %in% c("gu") ) {
    pheno_dt[,c("UD", "SD", "DD", "RD", "maxline", "baseline") := list(
      as.Date(UD, origin = begin)-1,
      as.Date(SD, origin = begin)-1,
      as.Date(DD, origin = begin)-1,
      as.Date(SD, origin = begin)-1,
      attr(data, "info")$rescale[1] + maxline * attr(data, "info")$rescale[2],
      attr(data, "info")$rescale[1] + baseline * attr(data, "info")$rescale[2]
    )]
  }

  # pheno_dt <- list(
  #   "dates" = rbindlist(lapply(pheno_x_list, rbindlist))[!is.na(date),]
  # )
  # if (sum(sapply(pheno_y_list, length)) > 0) {
  #   pheno_dt[["values"]] <- rbindlist(lapply(pheno_y_list, rbindlist))[!is.na(value),]
  # }
  attr(pheno_dt, "gen_by") <- "extract_pheno"
  attr(pheno_dt, "info") <- list(
    "method" = method
  )
  if (method == "trs") {attr(pheno_dt, "info")$trs <- trs}
  pheno_dt
  
}