revise according to Tim's advice

.gitmodules

@@ -1,3 +1,3 @@
 [submodule "scripts"]
 	path = scripts
-	url = https://github.com/kongdd/phenofit-scripts
+	url = https://github.com/eco-hydro/phenofit-scripts

NAMESPACE

@@ -74,6 +74,7 @@ export(lambda_cv_jl)
 export(lambda_vcurve)
 export(lambda_vcurve_jl)
 export(logistic)
+export(make_date)
 export(melt_list)
 export(movmean)
 export(opt_nlm)
@@ -164,6 +165,7 @@ importFrom(lubridate,ddays)
 importFrom(lubridate,dyears)
 importFrom(lubridate,is.Date)
 importFrom(lubridate,leap_year)
+importFrom(lubridate,make_date)
 importFrom(lubridate,month)
 importFrom(lubridate,yday)
 importFrom(lubridate,year)

R/PhenoExtractMeth.R

@@ -1,3 +1,7 @@
+# colors    <- c("blue", "green3", "orange", "red")
+colors <- c("green3", "darkgreen", "darkorange3", "red")
+linewidth <- 1.2
+
 #' @name PhenoExtractMeth
 #' @title Phenology Extraction methods
 #'
@@ -16,7 +20,7 @@
 #' # simulate vegetation time-series
 #' fFUN = doubleLog.Beck
 #' par  = c( mn  = 0.1 , mx  = 0.7 , sos = 50 , rsp = 0.1 , eos = 250, rau = 0.1)
-#' 
+#'
 #' t    <- seq(1, 365, 8)
 #' tout <- seq(1, 365, 1)
 #' y <- fFUN(par, t)
@@ -57,7 +61,7 @@ PhenoTrs <- function(fFIT, approach = c("White", "Trs"), trs = 0.5, #, min.mean
 
     # get peak of season position
     half.season <- median(which.max(values)) %>% round() # + 20, half season + 20 was unreasonable
-    pop <- t[half.season]
+    pos <- t[half.season]
 
     if (all(is.na(values))) return(metrics)
     if (half.season < 5 || half.season > (n - 5)) return(metrics)
@@ -107,8 +111,8 @@ PhenoTrs <- function(fFIT, approach = c("White", "Trs"), trs = 0.5, #, min.mean
     # 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 < pop], na.rm = TRUE))
-    eos <- round(median(t[!greenup & bool & t > pop], 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))
     # los <- eos - sos#los < 0 indicate that error
     # los[los < 0] <- n + (eos[los < 0] - sos[los < 0])
 
@@ -125,7 +129,7 @@ PhenoTrs <- function(fFIT, approach = c("White", "Trs"), trs = 0.5, #, min.mean
     #     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, pop = pop, mgs = mgs,
+    # 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
 
@@ -134,15 +138,14 @@ PhenoTrs <- function(fFIT, approach = c("White", "Trs"), trs = 0.5, #, min.mean
         PhenoPlot(t, values, 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)
-
+        # 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, POP, MGS, rsp, rau, PEAK, MSP and MAU. }
+    ### The function returns a vector with SOS, EOS, LOS, POS, MGS, rsp, rau, PEAK, MSP and MAU. }
 }
 
 # identify greenup or dormancy(brown) period
@@ -164,16 +167,16 @@ PhenoDeriv <- function(fFIT,
     analytical = TRUE, smoothed.spline = FALSE,
     IsPlot = TRUE, show.lgd = TRUE, ...)
 {
-    PhenoNames <- c("SOS", "POP", "EOS")
-    metrics <- setNames(rep(NA, 3), c("sos", "pop", "eos")) # template
+    PhenoNames <- c("SOS", "POS", "EOS")
+    metrics <- setNames(rep(NA, 3), c("sos", "pos", "eos")) # template
 
     t      <- fFIT$tout
     values <- last(fFIT$zs)
     n      <- length(t)
 
     # get peak of season position
-    half.season <- median(which.max(values)) # deal with multiple pop values
-    pop <- t[half.season]
+    half.season <- median(which.max(values)) # deal with multiple pos values
+    pos <- t[half.season]
 
     if (all(is.na(values))) return(metrics)
     if (half.season < 5 || half.season > (n - 5)) return(metrics)
@@ -197,29 +200,29 @@ PhenoDeriv <- function(fFIT,
     if (is_empty(sos)) sos <- NA
     if (is_empty(eos)) eos <- NA
 
-    metrics <- c(sos = sos, pop = pop, eos = eos)#, los = los
+    metrics <- c(sos = sos, pos = pos, eos = eos)#, los = los
 
     if (IsPlot){
         main  <- ifelse(all(par("mar") == 0), "", "DER")
         PhenoPlot(t, values, main = main, ...)
         if (show.lgd) legend('topright', c("f(t)'"), lty = 2, col = "black", bty='n')
 
         abline(v = c(sos, eos), col = colors[c(1, 4)], lwd = linewidth)
-        abline(v = pop, col ="darkgreen", lty = 1, lwd = linewidth)
+        abline(v = pos, col ="blue", lty = 1, lwd = linewidth)
 
         A <- diff(range(values))
         I_metrics <- match(metrics, t)
         if (all(is.na(I_metrics))) {
             ylons <- I_metrics
         }else{
-            ylons <- values[I_metrics] + c(1, -1, 1)*0.1*A
+            ylons <- values[I_metrics] + c(1, -2, 1)*0.1*A
         }
         xlons <- metrics + c(-1, 1, 1)*5
         xlons[xlons < min(t)] <- min(t)
         xlons[xlons > max(t)] <- max(t)
 
         I <- c(1); text(xlons[I], ylons[I], PhenoNames[I], col = colors[I], adj = c(1, 0))
-        I <- 2:3 ; text(xlons[I], ylons[I], PhenoNames[I], col = c("darkgreen", colors[3]), adj = c(0, 0))
+        I <- 2:3 ; text(xlons[I], ylons[I], PhenoNames[I], col = c("blue", colors[4]), adj = c(0, 0))
 
         #der1 last plot
         op <- par(new = TRUE)
@@ -245,8 +248,8 @@ PhenoGu <- function(fFIT,
     n      <- length(t)
 
     # get peak of season position
-    half.season <- median(which.max(values)) # deal with multiple pop values
-    pop <- t[half.season]
+    half.season <- median(which.max(values)) # deal with multiple pos values
+    pos <- t[half.season]
 
     if (all(is.na(values))) return(metrics)
     if (half.season < 5 || half.season > (n - 5)) return(metrics)
@@ -329,6 +332,7 @@ PhenoGu <- function(fFIT,
         xlons <- metrics[1:4] + c(-1, -1, 1, 1)*5
         xlons[xlons < min(t)] <- min(t)
         xlons[xlons > max(t)] <- max(t)
+
         I <- c(1, 2); text(xlons[I], ylons[I], PhenoNames[I], col = colors[I], adj = c(1, 0))
         I <- c(3, 4); text(xlons[I], ylons[I], PhenoNames[I], col = colors[I], adj = c(0, 0))
     }
@@ -352,7 +356,7 @@ PhenoKl <- function(fFIT,
 
     # get peak of season position
     half.season <- median(which.max(values)) # + 20, half season + 20 was unreasonable
-    pop <- t[half.season]
+    pos <- t[half.season]
 
     if (all(is.na(values))) return(metrics)
     if (half.season < 5 || half.season > (n - 5)) return(metrics)
@@ -429,12 +433,15 @@ PhenoKl <- function(fFIT,
             legend('topright', c("K'"), lty = c(3), col = c("black"), bty='n') ##pch =c(20, 1),
         }
 
-        pop     <- t[half.season]
+        pos     <- t[half.season]
         abline(v = metrics, col = colors, lwd = linewidth)
-        # abline(v = pop, col ="darkgreen", lty = 1, lwd = linewidth)
-        # abline(v = pop + 20, col ="darkgreen", lty = 2, lwd = linewidth)
-        I <- c(1, 3); text(xlons[I], ylons[I], PhenoNames[I], col = colors[I], adj = c(0, 0))
-        I <- c(2, 4); text(xlons[I], ylons[I], PhenoNames[I], col = colors[I], adj = c(1, 0))
+        # abline(v = pos, col ="darkgreen", lty = 1, lwd = linewidth)
+        # abline(v = pos + 20, col ="darkgreen", lty = 2, lwd = linewidth)
+        PhenoNames2 <- c("Greenup", "Maturity", "Senescence", "Dormancy")
+        # PhenoNames2 <- c("G", "M", "S", "D")
+
+        I <- c(1, 3); text(xlons[I], ylons[I], PhenoNames2[I], col = colors[I], adj = c(0, 0))
+        I <- c(2, 4); text(xlons[I], ylons[I], PhenoNames2[I], col = colors[I], adj = c(1, 0))
         # der.k last plot
         op <- par(new = TRUE)
         plot(t, der.k, xlim = xlim, type= "l",

R/S3_fFITs.R

@@ -72,7 +72,7 @@ plot.fFITs <- function(x, method, ...){
 
         pred   <- last(fFIT$zs)
 
-        pop    <- t[which.max(pred)]
+        pos    <- t[which.max(pred)]
         derivs <- curvature(fFIT)
 
         # plot for every method
@@ -84,22 +84,22 @@ plot.fFITs <- function(x, method, ...){
             type= "b", pch = 20, cex = 1.3, col = "grey",
             main = "curve fitting VI", xlab = "Index", ylab = "VI")
         lines(t, pred); grid()
-        abline(v = pop, col ="green")
+        abline(v = pos, col ="green")
 
         maxd_der1 <- t[which.max(derivs$der1)]
         mind_der1 <- t[which.min(derivs$der1)]
 
         plot(t, derivs$der1, main = "D1"); grid()
         abline(v = maxd_der1, col ="blue")
         abline(v = mind_der1, col ="red")
-        abline(v = pop, col ="green")
+        abline(v = pos, col ="green")
 
         plot(t, derivs$der2, main = "D2"); grid()
         plot(t, derivs$k, main = "k")    ; grid()
         abline(v = maxd_der1, col ="blue")
         abline(v = mind_der1, col ="red")
-        abline(v = pop, col ="green")
-        abline(v = pop + 20, col ="green", lty = 2)
+        abline(v = pos, col ="green")
+        abline(v = pos + 20, col ="green", lty = 2)
         # plot(diff(der1_diff), main = "diff2")
 
         # k <- derivs$k

R/curvefit.R

@@ -1,5 +1,5 @@
 phenonames <- c('TRS2.SOS', 'TRS2.EOS', 'TRS5.SOS', 'TRS5.EOS', 'TRS6.SOS', 'TRS6.EOS',
-    'DER.SOS', 'DER.POP', 'DER.EOS',
+    'DER.SOS', 'DER.POS', 'DER.EOS',
     'UD', 'SD', 'DD','RD',
     'GreenUp', 'Maturity', 'Senescence', 'Dormancy')
 

R/findpeaks.R

@@ -176,12 +176,11 @@ findpeaks_season_jl <- function(
     nyear = 1)
 {
     A = max(ypred) - min(ypred)
-    ans = JuliaCall::julia_call("findpeaks_season", 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,
         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
 }
-

R/get_pheno.R

@@ -1,6 +1,3 @@
-colors    <- c("blue", "green3", "orange", "red")
-linewidth <- 1.2
-
 # ' PhenoPlot
 # '
 # ' @inheritParams check_input
@@ -11,8 +8,8 @@ linewidth <- 1.2
 PhenoPlot <- function(t, y, main = "", ...){
     plot(t, y, main = main, ...,
              type= "l", cex = 2, col = "black", lwd = linewidth) #pch = 20,
-    # grid(nx = NA)
-    grid(ny = 4, nx = NA)
+    grid(nx = NA)
+    # grid(ny = 4, nx = NA)
 }
 
 #' get_pheno
@@ -52,9 +49,13 @@ get_pheno <- function(fits, method,
     res <- lapply(set_names(seq_along(methods), methods), function(k){
         method <- methods[k]
         if (IsPlot){
-            oma <- if (show_title) c(1, 2, 4, 1) else c(1, 2, 2, 1)
-            op <- par(mfrow = c(length(fits), 5), oma = oma,
-                mar = rep(0, 4), yaxt = "n", xaxt = "n")
+            op <- par(mfrow = c(length(fits), 5),
+                mgp = c(3, 0.6, 0), mar = rep(0, 4), yaxt = "n", xaxt = "n")
+            if (isTRUE(all.equal(par("oma"), c(0, 0, 0, 0)))) {
+                margin_l = 5.5
+                oma <- if (show_title) c(1, margin_l, 4, 1) else c(1, margin_l, 2, 1)
+                par(oma = oma)
+            }
         }
 
         # fFITs
@@ -121,7 +122,7 @@ get_pheno.fFITs <- function(fFITs, method,
         ylim     <- ylim0 + c(-1, 0.2) * 0.05 *A
         ylim_trs <- (ylim - ylim0) / A # TRS:0-1
 
-        PhenoPlot(fFITs$tout, ypred, ylim = ylim)
+        PhenoPlot(fFITs$tout, ypred, ylim = ylim, yaxt = "s")
         lines(ti, yi, lwd = 1, col = "grey60")
 
         QC_flag <- fFITs$data$QC_flag
@@ -153,9 +154,10 @@ get_pheno.fFITs <- function(fFITs, method,
         )
 
         legend('topleft', do.call(expression, exprs), adj = c(0.2, 0.2), bty='n', text.col = "red")
-        mtext(title_left, side = 2, line = 0.2)
+        # mtext(title_left, side = 2, line = 0.2)
+        mtext(title_left, side = 2, line = 1.8)
     }
-    if (showName_pheno && IsPlot) mtext("Fitting", line = 0.2)
+    if (showName_pheno && IsPlot) mtext("Fine fitting", line = 0.2)
 
     p_TRS <- lapply(TRS, function(trs) {
         PhenoTrs(fFIT, approach = "White", trs = trs, IsPlot = FALSE)
@@ -171,9 +173,14 @@ get_pheno.fFITs <- function(fFITs, method,
         IsPlot, ylim = ylim)
     param_common2 <- c(param_common, list(show.lgd = show.lgd))
 
-    der   <- do.call(PhenoDeriv, param_common2);  if (showName_pheno && IsPlot) mtext("DER", line = 0.2)
-    gu    <- do.call(PhenoGu, param_common)[1:4]; if (showName_pheno && IsPlot) mtext("GU", line = 0.2)
-    zhang <- do.call(PhenoKl, param_common2);     if (showName_pheno && IsPlot) mtext("ZHANG", line = 0.2)
+    der <- do.call(PhenoDeriv, param_common2)
+    if (showName_pheno && IsPlot) mtext("DER", line = 0.2)
+
+    zhang <- do.call(PhenoKl, param_common2)
+    if (showName_pheno && IsPlot) mtext("Inflexion ", line = 0.2)
+
+    gu <- do.call(PhenoGu, param_common)[1:4]
+    if (showName_pheno && IsPlot) mtext("Gu", line = 0.2)
 
     c(p_TRS, list(der, gu, zhang)) %>% set_names(methods)
 }

R/plot_curvefits.R

@@ -1,3 +1,5 @@
+# ' @param theme ggplot theme to be applied
+
 #' plot_curvefits
 #'
 #' @param d_fit data.frame of curve fittings returned by [get_fitting()].
@@ -10,11 +12,12 @@
 #' @param yticks ticks of y axis
 #' @param font.size Font size of axis.text
 #' @param show.legend Boolean
-#' @param theme ggplot theme to be applied
 #' @param shape the shape of input VI observation? `line` or `point`
 #' @param cex point size for VI observation.
 #' @param angle `text.x` angle
-#'
+#' @param layer_extra (not used) extra ggplot layers
+#' @param ... ignored
+#' 
 #' @example inst/examples/ex-curvefits.R
 #' 
 #' @export
@@ -29,7 +32,9 @@ plot_curvefits <- function(
     theme = NULL,
     cex = 2,
     shape = "point", angle = 30,
-    show.legend = TRUE)
+    show.legend = TRUE, 
+    layer_extra = NULL, 
+    ...)
 {
     methods <- d_fit$meth %>% unique() %>% rm_empty() # in case of NA
     nmethod <- length(methods) # how many curve fitting methods?
@@ -69,8 +74,7 @@ plot_curvefits <- function(
         p <- p + geom_point(
             data = d_obs,
             aes_string("t", "y", shape = "QC_flag", color = "QC_flag", fill = "QC_flag"),
-            size = cex, alpha = 0.7
-        )
+            size = cex, alpha = 0.7)
     } else {
         p <- if (shape == "point") {
             p + geom_point(data = d_obs, aes_string("t", "y"), size = cex, alpha = 0.6, color = "grey60")
@@ -92,13 +96,13 @@ plot_curvefits <- function(
         scale_shape_manual(values = qc_shapes, drop = F) +
         coord_cartesian(xlim = xlim)
 
-    if (!is.null(theme)) p <- p + theme
     if (!is.null(yticks)) p <- p + scale_y_continuous(breaks = yticks)
     if  (is.null(title)) p <- p + theme(plot.title = element_blank())
-
+    if (!is.null(layer_extra)) p <- p + layer_extra
+
     if (show.legend) {
-        iters_name_fine = c("Rough fitting", "Fine fitting")
-        lines_colors = c("black", "red")
+        iters_name_fine = c("Rough fitting", "iter1", "iter2")
+        lines_colors = c("black", "blue", "red")
         lgd <- make_legend_nmax(iters_name_fine, lines_colors, d_obs$QC_flag)
         p <- p + theme(legend.position = "none")
         p <- arrangeGrob(p, lgd,