improve findpeaks and GS division

DESCRIPTION

@@ -39,7 +39,7 @@ Depends:
 Imports: 
     Rcpp,
     purrr,
-    dplyr,
+    dplyr, stringr,
     magrittr,
     lubridate,
     data.table,

NAMESPACE

@@ -34,9 +34,9 @@ export(PhenoKl)
 export(PhenoTrs)
 export(R2_sign)
 export(add_HeadTail)
-export(cheak_season_list)
 export(check_input)
 export(check_season_dt)
+export(check_season_list)
 export(check_ylu)
 export(curvature)
 export(curvefit)
@@ -185,6 +185,7 @@ importFrom(methods,is)
 importFrom(purrr,map)
 importFrom(purrr,map_df)
 importFrom(purrr,map_dfc)
+importFrom(stringr,str_replace_all)
 importFrom(ucminf,ucminf)
 importFrom(utils,object.size)
 importFrom(utils,str)

NEWS.md

@@ -1,8 +1,11 @@
 # phenofit 0.3.4 (2021-11-14)
 
 - Fix the bug of `findpeaks`, which lead to sharp changed growing season failed to detect.
-- Fix the bug of `PhenoKlos` where `minpeakheight` not work in previous version.
-
+- Fix the bug of `PhenoKlos`, where `minpeakheight` not work in previous version.
+- Fix the bug of `check_season_dt`, where `peak` might be able to greater than `end`.
+- Remove the parameter `check_season_dt` in `removeClosedExtreme`, which might eliminate
+  good extreme values.
+
 # phenofit 0.3.3 (2021-10-26)
 
 - Fix typo error in curvefits' document.

R/PhenoExtractMeth.R

@@ -221,7 +221,7 @@ PhenoDeriv <- function(fFIT, t = NULL,
     }
     if (is_empty(eos)) {
         I_eos <- findpeaks(-der1[(half.season+5):length(der1)],
-                           nups=nmin, ndowns=0, npeaks=1, sortstr=TRUE)$X$pos + half.season + 4
+            nups=nmin, ndowns=0, npeaks=1, sortstr=TRUE)$X$pos + half.season + 4
         eos <- t[I_eos]
         if (is_empty(eos)) eos <- round((t[n] + pos)/2)
     }

R/deprecated_season.R

@@ -212,7 +212,7 @@ season <- function(
                     (c(FALSE, diff(val) < 0) & c(FALSE, diff(type) == 2)))
             pos <- pos[I, ]
 
-            pos   <- removeClosedExtreme(pos, ypred, A, y_min = r_min*A, minpeakdistance)
+            pos   <- removeClosedExtreme(pos, ypred, A, y_min = r_min*A)
             pos$t <- t[pos$pos]
             if (nrow(pos) < 2){ # at least two points, begin and end
                 warning("Can't find a complete growing season before!"); return(res)

R/findpeaks.R

@@ -20,18 +20,24 @@
 #' @param minpeakdistance The minimum distance (in indices) peaks have to have
 #' to be counted. If the distance of two maximum extreme value less than
 #' `minpeakdistance`, only the real maximum value will be left.
-#' @param y_min Threshold is defined as the difference of peak value with
-#' trough value. There are two threshold (left and right). The minimum threshold
-#' should be greater than `y_min`.
-#' @param y_max Similar as `y_min`, The maximum threshold should
-#' be greater than `y_max`.
+#' @param h_min `h` is defined as the difference of peak value to the
+#' adjacent left and right trough value (`h_left` and `h_right` respectively).
+#' The real peaks should follow `min(h_left, h_right) >= h_min`.
+#' @param h_max Similar as `h_min`, the real peaks should follow
+#' `max(h_left, h_right) >= h_min`.
 #' @param npeaks  the number of peaks to return. If `sortstr` = true, the
 #' largest npeaks maximum values will be returned; If `sortstr` = false,
 #' just the first npeaks are returned in the order of index.
 #' @param sortstr Boolean, Should the peaks be returned sorted in decreasing oreder of
 #' their maximum value?
 #' @param IsPlot Boolean.
 #'
+#' @note
+#' In versions before v0.3.4, `findpeaks(c(1, 2, 3, 4, 4, 3, 1))` failed to detect
+#' peaks when a flat pattern exit in the middle.
+#'
+#' From version v0.3.4, the peak pattern was changed from `[+]{%d,}[-]{%d,}` to
+#' `[+]{%d,}[0]{0,}[-]{%d,}`. The latter can escape the flat part successfully.
 #' @examples
 #' x <- seq(0, 1, len = 1024)
 #' pos <- c(0.1, 0.13, 0.15, 0.23, 0.25, 0.40, 0.44, 0.65, 0.76, 0.78, 0.81)
@@ -43,14 +49,16 @@
 #' }
 #'
 #' plot(pSignal, type="l", col="navy"); grid()
-#' x <- findpeaks(pSignal, npeaks=3, y_min=4, sortstr=TRUE)
+#' x <- findpeaks(pSignal, npeaks=3, h_min=4, sortstr=TRUE)
 #' points(val~pos, x$X, pch=20, col="maroon")
 #'
 #' @export
 findpeaks <- function (x, nups = 1, ndowns = nups, zero = "0", peakpat = NULL,
-                       minpeakheight = -Inf, minpeakdistance = 1,
-                       y_min = 0, y_max = 0,
-                       npeaks = 0, sortstr = FALSE, IsPlot = F)
+    minpeakheight = -Inf, minpeakdistance = 1,
+    h_min = 0, h_max = 0,
+    npeaks = 0, sortstr = FALSE,
+    include_gregexpr = FALSE,
+    IsPlot = F)
 {
     stopifnot(is.vector(x, mode = "numeric") ||
                   is.vector(x, mode = "logical") || length(is.na(x)) == 0)
@@ -67,10 +75,13 @@ findpeaks <- function (x, nups = 1, ndowns = nups, zero = "0", peakpat = NULL,
     #     xc <- x
     # }
     xc <- sign(diff(x))
-    xc <- paste(as.character(sign(xc)), collapse = "")
-    xc <- gsub("1", "+", gsub("-1", "-", xc))
-    if (zero != "0")      xc      <- gsub("0", zero, xc)
-    if (is.null(peakpat)) peakpat <- sprintf("[+]{%d,}[-]{%d,}", nups, ndowns)
+    xc <- paste(as.character(sign(xc)), collapse = "") %>%
+        {gsub("1", "+", gsub("-1", "-", .))}
+    xc %<>% str_replace_midzero() # in v0.3.4
+
+    if (zero != "0") xc <- gsub("0", zero, xc)
+    # if (is.null(peakpat)) peakpat <- sprintf("[+]{%d,}[-]{%d,}", nups, ndowns)
+    if (is.null(peakpat)) peakpat <- sprintf("[+]{%d,}[0]{0,}[-]{%d,}", nups, ndowns)
 
     # Fatal bug found at 20211114
     # Because `diff` operation lead to the length of `x` reduced one
@@ -93,8 +104,8 @@ findpeaks <- function (x, nups = 1, ndowns = nups, zero = "0", peakpat = NULL,
         xv[i] <- x[xp[i]]
     }
     inds <- which(xv >= minpeakheight &
-                      xv - pmin(x[x1], x[x2]) >= y_max &
-                      xv - pmax(x[x1], x[x2]) >= y_min)
+                      xv - pmin(x[x1], x[x2]) >= h_max &
+                      xv - pmax(x[x1], x[x2]) >= h_min)
     X <- cbind(xv[inds], xp[inds], x1[inds], x2[inds])
 
     if (length(X) == 0) return(NULL)
@@ -128,9 +139,19 @@ findpeaks <- function (x, nups = 1, ndowns = nups, zero = "0", peakpat = NULL,
         plot(x, type = "b"); grid()
         points(val~pos, X, col = "blue")
     }
-    return(list(gregexpr = rc, X = X))
+    if (include_gregexpr) listk(X, gregexpr = rc) else listk(X)
 }
 
+#' @importFrom stringr str_replace_all
+str_replace_midzero <- function(x) {
+    replace = function(x, replacement = "+") {
+        paste(rep(replacement, nchar(x)), collapse = "")
+    }
+    str_replace_all(x, "\\++0\\++", . %>% replace("+")) %>%
+        str_replace_all("-+0-+", . %>% replace("-"))
+}
+
+
 findpeaks_season <- function(
     ypred, r_max = 0, r_min = 0,
     minpeakdistance = 0, minpeakheight = 0,
@@ -139,13 +160,13 @@ findpeaks_season <- function(
     nyear = 1)
 {
     A = max(ypred) - min(ypred)
-    y_max = r_max * A
-    y_min = r_min * A
+    h_max = r_max * A
+    h_min = r_min * A
     # local minimum values
     # peak values is small for minimum values, so can't use r_min here
     peaks <- findpeaks(-ypred,
-        zero = "-",
-        y_max = y_max, y_min = y_min * 0,
+        zero = "+",
+        h_max = h_max, h_min = h_min * 0,
         minpeakdistance = minpeakdistance,
         nups = 0
     )
@@ -161,7 +182,7 @@ findpeaks_season <- function(
     # local maximum values,
     peaks <- findpeaks(ypred,
         zero = "+",
-        y_max = y_max, y_min = y_min,
+        h_max = h_max, h_min = h_min,
         minpeakdistance = minpeakdistance,
         minpeakheight = minpeakheight,
         nups = nups, ndowns = ndowns
@@ -184,9 +205,9 @@ findpeaks_season_jl <- function(
     A = max(ypred) - min(ypred)
     ans = JuliaCall::julia_call("phenofit.findpeaks_season", ypred,
         r_max = r_max, r_min = r_min,
-        # r_max = y_max/A, r_min = y_min/A,
+        # r_max = h_max/A, r_min = h_min/A,
         minpeakdistance = as.integer(minpeakdistance), minpeakheight = minpeakheight,
         nups = nups, ndowns = ndowns)
-    ans$threshold <- data.table(y_max = r_max*A, y_min = r_min*A, r_max, r_min)
+    ans$threshold <- data.table(h_max = r_max * A, h_min = r_min * A, r_max, r_min)
     ans
 }

R/global_options.R

@@ -82,19 +82,19 @@ options_fitting <- list(
 #' @param ... list of phenofit options
 #' FUN_season: character, `season_mov` or `season`
 #' rFUN: character, rough fitting function. `smooth_wWHIT`, `smooth_wSG` or `smooth_wHANTs`.
-#'
+#' @param options If not NULL, `options` will be used and `...` will be ignored.
 #' @examples
 #' set_options(verbose = FALSE)
 #' get_options("season") %>% str()
 #' @export
-set_options <- function(...) {
-    opts = list(...)
+set_options <- function(..., options = NULL) {
+    if (is.null(options)) options = list(...)
     # rm unrelated parameters
-    ind = match(names(opts), names(.options)) %>% which.notna()
+    ind = match(names(options), names(.options)) %>% which.notna()
     # This step might lead to error, but will improve performance
     if (.options$initialized && length(ind) == 0) return()
-
-    .options %<>% modifyList(opts[ind])
+    
+    .options %<>% modifyList(options[ind])
     # `season` and `fitting` share the same parameter
     pars_comm = c("wFUN", "wmin", "verbose") 
     for (par in pars_comm) {

R/qcFUN.R

@@ -245,8 +245,8 @@ qc_SPOT <- function (QA, wmin = 0.2, wmid = 0.5, wmax = 1) {
 #' @examples
 #' qc_sentinel2(1:11)
 #' @export
-qc_sentinel2 <- function(SLC, wmin = 0.2, wmid = 0.5, wmax = 1) {
-    QA <- SLC
+qc_sentinel2 <- function(SCL, wmin = 0.2, wmid = 0.5, wmax = 1) {
+    QA <- SCL
     qc_good = c(4, 5, 6, 7, 10)
     qc_bad = c(1, 2, 3, 9)
     qc_mid = c(8)

R/season_input.R

@@ -62,7 +62,7 @@ find_season.peaks <- function(
                     (c(FALSE, diff(val) < 0) & c(FALSE, diff(type) == 2)))
             pos <- pos[I, ]
 
-            pos   <- removeClosedExtreme(pos, ypred, A, y_min = r_min*A, minpeakdistance)
+            pos   <- removeClosedExtreme(pos, ypred, A, y_min = r_min*A)
             pos$t <- t[pos$pos]
             if (nrow(pos) < 2){ # at least two points, begin and end
                 warning("Can't find a complete growing season before!"); return(NULL)
@@ -146,11 +146,9 @@ season_input <- function(INPUT, options = NULL, verbose = FALSE, ...)
     if (verbose) print(str(.options$season))
 
     opt = get_options()
-#     browser()
-# assign("opt", opt, envir = globalenv())
-# assign("INPUT", INPUT, envir = globalenv())
-    c(d_fit, info_peak) %<-% roughFit(INPUT)
-    d_season = find_season.peaks(d_fit, info_peak)
-    listk(fit = d_fit, dt = d_season)
+    # export2glob(opt, INPUT)
+    c(rfit, info_peak) %<-% roughFit(INPUT)
+    d_season = find_season.peaks(rfit, info_peak)
+    listk(fit = rfit, dt = d_season)
 }
 

R/season_mov.R

@@ -85,13 +85,12 @@ season_mov <- function(INPUT,
 
     brks  <- list()
     vcs   <- vector("list", length(years.run)) %>% set_names(years.run)
-
     for (year_i in years.run) {
-        i = which(year_i == years)
+        i = which(year_i == years.run)
         if (.options$season$verbose) fprintf("  [season_mov] running %d ... \n", i)
-        
+
         # I <- which(date_year %in% years[(i-ny_extend):(i+ny_extend)]) # 3y index
-        I   <- which(date_year %in% years[i])
+        I   <- which(date_year %in% years.run[i])
         ylu <- get_ylu(INPUT$y, date_year, INPUT$w, width = width_ylu, I, Imedian = TRUE,
             opt$wmin)
         ylu %<>% merge_ylu(INPUT$ylu, .) # curvefits.R
@@ -118,7 +117,15 @@ season_mov <- function(INPUT,
         if (!is_empty(brk$dt))
             brk$dt %<>% subset(year == year_i) %>% mutate(lambda = lambda)
 
-        ans = list(fit = brk$fit[date_year[I] == year_i, ], dt = brk$dt)
+        # improve for head and tail, v0.3.4
+        if (i == 1) {
+            rfit = brk$fit[date_year[I] <= year_i, ]
+        } else if (i == length(years.run)) {
+            rfit = brk$fit[date_year[I] >= year_i, ]
+        } else {
+            rfit = brk$fit[date_year[I] == year_i, ]
+        }
+        ans = list(fit = rfit, dt = brk$dt)
         if (.options$debug) ans %<>% c(list(fit.raw = brk$fit), .)
         brks[[i]] <- ans
     }
@@ -133,11 +140,12 @@ season_mov <- function(INPUT,
         if (is_empty(dt)) {
             warning( 'No growing season found!'); return(NULL)
         }
-        if (opt$.check_season) brks$dt %<>% cheak_season_list()
+        if (opt$.check_season) brks$dt %<>% check_season_list()
     }
     brks$GOF <- stat_season(INPUT, brks$fit)
-
-    if (opt$rFUN == "smooth_wWHIT" && !has_lambda && opt$.lambda_vcurve) 
+    # compatible for previous versions
+    if (is.character(opt$rFUN) && opt$rFUN == "smooth_wWHIT" &&
+        !has_lambda && opt$.lambda_vcurve)
         brks$optim <- vcs
     return(brks)
 }
@@ -174,7 +182,8 @@ check_season_dt <- function(dt) {
     rcpp_season_filter(dt, .options$season$rm.closed,
         rtrough_max = .options$season$rtrough_max,
         r_min = .options$season$r_min)
-    dt[y_peak >= .options$season$ypeak_min & (len > len_min & len < len_max), ]
+    dt[y_peak >= .options$season$ypeak_min & (len > len_min & len < len_max), ] %>%
+        .[peak <= end & beg < end] # a temporary resolution
 }
 
 #' @param dt data.table of growing season dividing info
@@ -183,7 +192,7 @@ check_season_dt <- function(dt) {
 #' @inheritParams season_mov
 #' @rdname rcpp_season_filter
 #' @export
-cheak_season_list <- function(lst_dt) {
+check_season_list <- function(lst_dt) {
     if (is.data.frame(lst_dt)) lst_dt = list(lst_dt)
     lst_dt %<>% rm_empty()
     dt2 = map(lst_dt, ~ check_season_dt(.x)) %>%