R/fit_curve.R

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

#' @title Covert a `zoo` into a `{phenopix}` format
#' @description Internal function to convert a `zoo` time series in the format
#'  returned by functions like `phenopix::ElmoreFit()` or `phenopix::GuFit()`.
#' @param ts A `ts` or `zoo` object with `gcc` data. 
#'  index (`ts`) must be numeric days of year (`DOY`) or a `POSIXct` vector
#' @param uncert Unused (left for compatibility).
#' @param nrep Unused (left for compatibility).
#' @return A list containing the following items.
#'  - `fit`: A list as returned by the function `FitDoubleLogElmore`
#'  - `uncertainty`: A list containing a `zoo` `data.frame` with the
#'      uncertainty predicted values, and a dataframe containing the 
#'      respective uncertainty curve parameters
#' @author Luigi Ranghetti, PhD (2020) \email{luigi@@ranghetti.info}
fakeFit <- function (ts, uncert=NA, nrep=NA) {
  list(
    "fit" = list(
      "predicted" = ts,
      "params" = NULL,
      "formula" = NULL,
      "sf" = NULL
    ),
    "uncertainty" = NULL
  )
}


#' @title Fit curve over cycles
#' @description Fit a curve using a parametric function from 
#'  `phenopix::greenProcess()`.
#' @param ts Time series in `s2ts` format (generated using [`fill_s2ts()`]).
#' @param cycles Cycles allocation (data table generated by [`cut_cycles()`],
#'  eventually filtered using [`assign_season()`]).
#' @param fit Fitting function among `"klosterman"`, `"beck"`,
#'  `"elmore"`, `"gu"` and `"no"` (see `phenopix::greenProcess()`).
#'  Value `"no"` can be used to convert a `s2ts` time series in the format
#'  returned by this function without making any interpolation.
#'  A vector of length > 1 can be specified, in which case the second element
#'  is used if the first method fails.
#' @return A named list of `n` elements (being `n` the number of IDs included
#'  in `ts`). Each element is a named list of `m` elements (being `m` the
#'  number of detected cycles).
#'  Each element is a list including the following elements:
#'  - `fit`: the fitted curve as returned by `phenopix` fitting functions;
#'  - `ts`: a data table of fitted values, containing the columns
#'      `date` (date of the record) and
#'      `value` (fitted value);
#'  - `maxval`: the date of the maximum value in the cycles
#'      (obtained from input `cycles` and here reported).
#' @author Luigi Ranghetti, PhD (2020) \email{luigi@@ranghetti.info}
#' @import data.table
#' @importFrom zoo zoo
#' @importFrom phenopix KlostermanFit BeckFit ElmoreFit GuFit
#' @importFrom utils txtProgressBar setTxtProgressBar
#' @export
#' @examples
#' # Load input data
#' data("ts_filled")
#' data("dt_cycles")
#' 
#' \donttest{
#' # Standard interpolation
#' cf <- fit_curve(ts_filled, dt_cycles)
#' plot(ts_filled, fitted = cf)
#' }
#' 
#' # Interpolation using a different method
#' cf_2 <- fit_curve(ts_filled, dt_cycles, fit = "no")


fit_curve <- function(
  ts,
  cycles,
  fit = c("gu", "klosterman")
) {
  
  # Avoid check notes for data.table related variables
  id <- value <- relval <- cycle <- begin <- maxval <- end <- NULL
  
  ## Check arguments
  # TODO
  
  # Check fit
  if (any(!fit %in% c("beck", "elmore", "klosterman", "gu", "no"))) {
    print_message(
      type = "error",
      "Argument 'fun' only accepts values \"klosterman\", \"beck\", ",
      "\"elmore\", \"gu\" and \"no\"."
    )
  }
  
  ## Check s2ts format
  # (must contain date, id, orbit, sensor, value, opt. quality)
  if (!inherits(ts, "s2ts")) {
    print_message(
      type = "error",
      "Argument 'ts' is not in the right format."
    )
  }
  # TODO
  ts_dt <- as.data.table(ts)
  
  # Check the frequency to be daily
  if (any(unique(ts_dt[,list(freq = diff(date)), by=id]$freq) != 1)) {
    print_message(
      type = "error",
      "Argument 'ts' must be a daily s2ts object, ",
      "generated with fill_s2ts(..., frequency = \"daily\")."
    )
  }
  
  # Compute relative values
  rescale <- ts_dt[,c(min(value, na.rm=TRUE), diff(range(value, na.rm=TRUE)))]
  ts_dt[,relval := (value - rescale[1]) / rescale[2]]
  
  # Assign fit.funs
  fit.funs <- lapply(fit, function(f) {
    if (f=="beck") phenopix::BeckFit else
      if (f=="elmore") phenopix::ElmoreFit else
        if (f=="klosterman") phenopix::KlostermanFit else
          if (f=="gu") phenopix::GuFit else
            if (f=="no") fakeFit
  })
  fit.funs <- list(
    "beck" = phenopix::BeckFit,
    "elmore" = phenopix::ElmoreFit,
    "klosterman" = phenopix::KlostermanFit,
    "gu" = phenopix::GuFit,
    "no" = fakeFit
  )
  
  # Set progress bar (time consuming function)
  use_pb <- inherits(stdout(), "terminal") && cycles[,length(unique(paste(id,year,cycle)))] > 1
  if (use_pb) {
    pb <- txtProgressBar(0, cycles[,length(unique(paste(id,year,cycle)))], style = 3)
  }
  
  # Fit for each ID/year/cycle
  fit_out <- list()
  for (sel_id in unique(ts_dt$id)) {
    fit_out[[sel_id]] <- list()
    for (sel_year in cycles[id == sel_id, unique(year)]) {
    fit_out[[sel_id]][[as.character(sel_year)]] <- list()
    for (sel_cycle in cycles[id == sel_id & year == sel_year, unique(cycle)]) {
      sel_cut_date <- as.Date(as.matrix(cycles[id == sel_id & year == sel_year & cycle == sel_cycle,][,list(begin, maxval, end)])[1,])
      sel_ts_zoo <- zoo::zoo(
        ts_dt[id == sel_id & date >= sel_cut_date[1] & date < sel_cut_date[3], relval]
      )
      sel_fit_success <- FALSE
      for (sel_fit in fit) {
        if (!sel_fit_success) {
          sel_fit_out <- try(fit.funs[[sel_fit]](ts = sel_ts_zoo, uncert = FALSE), silent = TRUE)
          attr(sel_fit_out, "info") <- list("method" = sel_fit)
          sel_fit_success <- !inherits(sel_fit_out, "try-error") &&
            !all(is.na(sel_fit_out$fit$predicted))
        }
      }
      if (sel_fit_success) {
        fit_out[[sel_id]][[as.character(sel_year)]][[as.character(sel_cycle)]] <- list(
          "fit" = sel_fit_out$fit,
          "ts" = data.table(
            "date" = ts_dt[
              id == sel_id & date >= sel_cut_date[1] & date < sel_cut_date[3], 
              date
              ],
            "value" = rescale[1] + as.numeric(sel_fit_out$fit$predicted) * rescale[2]
          ),
          "maxval" = sel_cut_date[2],
          "weight" = cycles[id == sel_id & year == sel_year & cycle == sel_cycle,]$weight,
          "info" = attr(sel_fit_out, "info")
        )
      }
      if (use_pb) {setTxtProgressBar(pb, pb$getVal()+1)} # update progress bar
    } # end of sel_cycle FOR cycle
    } # end of sel_year FOR cycle
  } # end of id FOR cycle
  
  if (use_pb) {message("")} # separate progress bar
  
  attr(fit_out, "gen_by") <- "fit_curve"
  attr(fit_out, "info") <- list("rescale" = rescale)
  fit_out
  
}