Skip to content

Commit

Permalink
major update: add get_pheno.rfit
Browse files Browse the repository at this point in the history
  • Loading branch information
kongdd committed Nov 21, 2021
1 parent 3e3d8a5 commit 6302763
Show file tree
Hide file tree
Showing 33 changed files with 1,094 additions and 736 deletions.
13 changes: 11 additions & 2 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -2,19 +2,27 @@

S3method(D1,fFIT)
S3method(D2,fFIT)
S3method(PhenoDeriv,default)
S3method(PhenoDeriv,fFIT)
S3method(PhenoGu,default)
S3method(PhenoGu,fFIT)
S3method(PhenoTrs,default)
S3method(PhenoTrs,fFIT)
S3method(curvature,fFIT)
S3method(get_GOF,fFITs)
S3method(get_GOF,list)
S3method(get_fitting,fFITs)
S3method(get_fitting,list)
S3method(get_param,fFITs)
S3method(get_param,list)
S3method(get_pheno,fFITs)
S3method(get_pheno,list)
S3method(get_pheno,rfit)
S3method(plot,fFITs)
S3method(print,fFIT)
export("%<-%")
export("%<>%")
export("%>%")
export(.Greenup)
export(D1)
export(D2)
export(FitDL.AG)
Expand All @@ -34,6 +42,7 @@ export(PhenoKl)
export(PhenoTrs)
export(R2_sign)
export(add_HeadTail)
export(brks2rfit)
export(check_input)
export(check_season_dt)
export(check_season_list)
Expand Down Expand Up @@ -68,7 +77,6 @@ export(get_fitting)
export(get_options)
export(get_param)
export(get_pheno)
export(get_pheno.fFITs)
export(init_AG)
export(init_AG2)
export(init_Beck)
Expand Down Expand Up @@ -182,6 +190,7 @@ importFrom(lubridate,yday)
importFrom(lubridate,year)
importFrom(lubridate,ymd)
importFrom(methods,is)
importFrom(purrr,as_mapper)
importFrom(purrr,map)
importFrom(purrr,map_df)
importFrom(purrr,map_dfc)
Expand Down
23 changes: 14 additions & 9 deletions NEWS.md
Original file line number Diff line number Diff line change
@@ -1,12 +1,4 @@
# 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 `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)
# phenofit 0.3.3 (2021-11-14)

- Fix typo error in curvefits' document.

Expand All @@ -30,6 +22,19 @@
used the strict mathematical solution to find the extreme values in the curve
of curvature's change rate.

** MAJOR updates to improve multi-GS phenology extraction **

- 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 `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.

- add `get_pheno.rfit` to extract vegetation phenology from rough fitting directly.

# phenofit 0.3.2 (2021-10-15)

- Parameters of `season_mov` and `curvefits` are wrapped into options. Scripts of phenofit v2.0 will not work anymore.
Expand Down
108 changes: 108 additions & 0 deletions R/PhenoDeriv.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,108 @@
#' Phenology extraction in Derivative method (DER)
#'
#' @inheritParams D
#' @inheritParams PhenoTrs
#' @param der1 the first order difference
#'
#' @references
#' 1. Filippa, G., Cremonese, E., Migliavacca, M., Galvagno, M., Forkel, M.,
#' Wingate, L., … Richardson, A. D. (2016). Phenopix: A R package for
#' image-based vegetation phenology. Agricultural and Forest Meteorology,
#' 220, 141–150. \doi{10.1016/j.agrformet.2016.01.006}
#'
#' @seealso [PhenoTrs()], [PhenoGu()], [PhenoKl()]
#' @export
PhenoDeriv <- function(x, t, ...) UseMethod("PhenoDeriv", x)

#' @rdname PhenoDeriv
#' @export
PhenoDeriv.fFIT <- function(x, t = NULL,
analytical = TRUE, smoothed.spline = FALSE, ...) {
if (!is.null(x$tout)) t <- x$tout
values <- last2(x$zs)
der1 <- D1.fFIT(x, t, analytical, smoothed.spline)
PhenoDeriv.default(values, t, der1, ...)
}

#' @rdname PhenoDeriv
#' @param show.legend whether show figure lelend?
#' @export
PhenoDeriv.default <- function(x, t, der1,
IsPlot = TRUE, show.legend = TRUE, ...)
{
PhenoNames <- c("SOS", "POS", "EOS")
metrics <- setNames(rep(NA, 3), c("sos", "pos", "eos")) # template
n <- length(t)

# get peak of season position
half.season <- median(which.max(x)) # deal with multiple pos x
pos <- t[half.season]

if (all(is.na(x))) return(metrics)
if (half.season < 5 || half.season > (n - 5)) return(metrics)

# get SOS and EOS according to first order derivative
# fixed 20180510, ±5 to make sure sos and eos are not in the POP.
# I_sos <- median(which.max(der1[1:(half.season - 5)]))
# I_eos <- median(which.min(der1[(half.season+5):length(der1)])) + half.season

# der.sos (eos) is impossible to occur near the pop.
nmin <- 5
I_sos <- findpeaks( der1[1:(half.season - 5)],
nups=nmin, ndowns=nmin, npeaks=1, sortstr=TRUE)$X$pos
I_eos <- findpeaks(-der1[(half.season+5):length(der1)],
nups=nmin, ndowns=nmin, npeaks=1, sortstr=TRUE)$X$pos + half.season + 4

# if half.season > length(der1), error will be occur
sos <- t[I_sos]
eos <- t[I_eos]
if (is_empty(sos)) {
I_sos <- findpeaks( der1[1:(half.season - 5)],
nups=nmin, ndowns=0, npeaks=1, sortstr=TRUE)$X$pos
sos <- t[I_sos]
if (is_empty(sos)) sos <- round((t[1] + pos)/2)
}
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
eos <- t[I_eos]
if (is_empty(eos)) eos <- round((t[n] + pos)/2)
}
# greenup <- .Greenup(x)
# if (is.na(sos)) {
# if (all(!rm_empty(greenup))) sos = pos
# }
# if (is.na(eos)) {
# if (all(rm_empty(greenup))) eos = pos
# }
metrics <- c(sos = sos, pos = pos, eos = eos)#, los = los

if (IsPlot){
main <- ifelse(all(par("mar") == 0), "", "DER")
PhenoPlot(t, x, main = main, ...)
if (show.legend) 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 = pos, col ="blue", lty = 1, lwd = linewidth)

A <- diff(range(x))
I_metrics <- match(metrics, t)
if (all(is.na(I_metrics))) {
ylons <- I_metrics
}else{
ylons <- x[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("blue", colors[4]), adj = c(0, 0))

#der1 last plot
op <- par(new = TRUE)
plot(t, der1, type= "l", lty = 2, lwd = linewidth,
col = "black", axes = FALSE)
}
return(metrics)
}
Loading

0 comments on commit 6302763

Please sign in to comment.