R/detect_forages.R

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

#' @title Detect forages
#' @description Experimental function to detect forages (time series characterised
#'  by a high "instability").
#'  Already detected cycles are labelled as "forage" (adjacent ones are dissolved)
#'  basing on arguments `diff_thr` and `ncuts_thr` (see).
#' @param ts Time series in `s2ts` format (generated using `fill_s2ts()`).
#' @param pheno Cycle allocation (data table generated by `extract_pheno()`
#'  or `cut_cycles()`).
#' @param win List with the allowed ranges for the dates of
#'  start of cycle (see `pop_seasons` for details about the format).
#' @param diff_thr First threshold to detect forages (TODO improve description)
#' @param ncuts_thr Second threshold to detect forages (TODO improve description)
#' @param sliding Length (in days) of the sliding window in which `diff_thr` 
#'  and `ncuts_thr` must be evaluated.
#' @param reldiff Logical: should argument `diff_thr` refer to relative (default)
#'  or absolute values?
#' @param byyear Logical: compute aggregation by year (default) or globally
#'  (this should make sense in case of permanent crops).
#' @return A data table in the format of the input `pheno`, 
#' in which adjacent forage cycles are merged.
#' An additional logical field `forage` discriminates forages from non-forages.
#' @author Luigi Ranghetti, PhD (2021) \email{luigi@@ranghetti.info}
#' @import data.table


detect_forages <- function(
  ts,
  pheno,
  win = NULL, # window,
  diff_thr = 10,
  ncuts_thr = 10,
  sliding = 90,
  reldiff = TRUE,
  byyear = TRUE
) {
  
  # Avoid check notes for data.table related variables
  # ...
  
  ## Check arguments
  # TODO
  if (is.null(win)) {win <- c("01-01","12-31")}
  
  # Check that ts is a filled ts
  if (is.null(attr(ts, "gen_by")) || attr(ts, "gen_by") != "fill_s2ts") {
    print_message(type="error", "'ts' must be generated by function fill_s2ts().")
  }
  # Check that pheno is generated by cut_cycles
  if (is.null(attr(pheno, "gen_by")) || 
      attr(pheno, "gen_by") != "cut_cycles" ||
      is.null(attr(pheno, "weight"))) {
    print_message(type="error", "'pheno' must be generated by function cut_cycles().")
  }
  
  
  ## Filter relevant cycles
  pheno_dt <- copy(pheno)
  `%op%` <- if (package_version(win)[1] < package_version(win)[2]) {`&`} else {`|`}
  pheno_sub <- pheno_dt[
    (package_version(strftime(maxval, "%m.%d")) >= package_version(win)[1]) %op%
      (package_version(strftime(maxval, "%m.%d")) <= package_version(win)[2]),
  ]
  
  ## Build relative TS
  tsd <- as.data.table(copy(ts))
  setnames(tsd, "value", "indexvalue", skip_absent = TRUE)
  tsd[,relindexval := (indexvalue - min(indexvalue, na.rm=TRUE)) / 
        diff(range(indexvalue, na.rm=TRUE)), by = id]
  tsd[,value := c(NA, diff(indexvalue)), by = id]
  tsd[,relval := c(NA, diff(relindexval)), by = id]
  tsd[, bg0 := relval <= 1e-2 & c(relval[-1], NA) >= 1e-2, by = id] # begin of growing
  tsd[, es0 := relval <= -1e-2 & c(relval[-1], NA) >= -1e-2, by = id] # end of senescence
  tsd[bg0|es0, bg1 := bg0 & !c(NA,bg0[-length(bg0)]), by = id]
  tsd[bg0|es0, es1 := es0 & !c(es0[-1],NA), by = id]
  tsd[, c("cut","bg0","es0") := list(bg1|es1,NULL,NULL)]

  ## Sliding windows
  ref_variable <- if (reldiff==TRUE) {"relval"} else {"value"}
  tsd[,rollval:=frollmean(abs(get(ref_variable)), sliding, align="center"), by = list(id)]

  ## Subsample ts into relevant cycles
  tsd_sub_l <- lapply(seq_len(nrow(pheno_sub)), function(j) {
    tsd[id == pheno_sub[j,id] & date>=pheno_sub[j,begin] & date<pheno_sub[j,end],][
      ,c("year","cycle"):=list(pheno_sub[j,year],pheno_sub[j,cycle])
    ]
  })
  tsd_sub <- rbindlist(tsd_sub_l)
  
  
  # TODO assunti (da passare a argomenti): relval/value, byyear
  diffsum <- tsd_sub[,list(
    diff = quantile(rollval,0.75,na.rm=TRUE)*365,
    ncuts = sum(bg1|es1,na.rm=TRUE)/.N*365
  ),by = list(id,year)]
  # diffsum[,aa := agreah_sel[match(id, agreah_sel$uid),]$aa]
  # diffsum[,colt := agreah_sel[match(id, agreah_sel$uid),]$DESC_SUOLO]
  # diffsum[,classerappr := (.N>10), by=colt]

  diffsum[,forage1:=(diff>diff_thr & ncuts>ncuts_thr)] # selezionato
  # diffsum[,table(aa,forage1)]
  # diffsum[,(sum(aa&forage1)+sum(!aa&!forage1))/.N]
  
  # Associazione diffsum-pheno
  pheno_sub1 <- merge(pheno_sub, diffsum[,list(id, year, forage1)], by=c("id","year"), all.x=TRUE)
  
  pheno_diffsum <- merge(
    pheno_sub1, 
    tsd_sub[,list(
      diff = quantile(rollval,0.75,na.rm=TRUE)*365
    ),by = list(id,year,cycle)], 
    by=c("id","year","cycle"), all.x=TRUE
  )
  
  weight_aggr <- if (attr(pheno_dt, "weight") %in% c("maxval")) {
    max
  } else {
    sum
  }
  pheno_diffsum_minmax <- pheno_diffsum[
    diff>diff_thr, 
    list(
      cycle = min(cycle), begin=min(begin), end=max(end), maxval = as.Date(NA), 
      weight = weight_aggr(weight), forage = TRUE, 
      min_maxval=min(maxval), max_maxval=max(maxval)
    ), 
    by=list(id, year)
  ]
  for (j in seq_len(nrow(pheno_diffsum_minmax))) {
    pheno_diffsum_minmax[j, maxval := tsd[
      id==id & date>=begin & date<=end,][
        relindexval==max(relindexval), median(date)
      ]] 
  }
  pheno_out <- rbind(
    pheno_diffsum_minmax[,list(id, year, cycle, begin, end, maxval, weight, forage)],
    merge(
      pheno_dt, 
      pheno_diffsum_minmax[,list(id, year, min_maxval, max_maxval)], 
      by = c("id", "year"), 
      all.x=TRUE
    )[
      is.na(min_maxval) | maxval<min_maxval | maxval>max_maxval, 
      list(id, year, cycle, begin, end, maxval, weight, forage=FALSE)
    ]
  )[order(id, year, begin),]
  
  # Return
  attr(pheno_out, "gen_by") <- "detect_forages"
  pheno_out  
  
}