phenofit

Description

Vegetation phenology package

Author(s)

Maintainer: Dongdong Kong kongdd.sysu@gmail.com

Authors:

• Mingzhong Xiao xmingzh@mail2.sysu.edu.cn

• Yongqiang Zhang yongqiang.zhang2014@gmail.com

• Xihui Gu guxh@cug.edu.cn

• Jianjian Cui cuijj6@mail2.sysu.edu.cn

See Also

Useful links:

• https://github.com/eco-hydro/phenofit

• Report bugs at https://github.com/eco-hydro/phenofit/issues

MOD13A1 EVI observations at flux site CA-NS6

Description

Variables in CA-NS6:

• site: site name

• y: EVI

• date: date of image

• t: date of compositing image

• w: weights of data point

• QC_flag: QC flag of y, in the range of c("snow", "cloud", "shadow", "aerosol", "marginal", "good")

Usage

data('CA_NS6')

Format

An object of class data.table (inherits from data.frame) with 161 rows and 6 columns.

D

Description

Get derivative of phenofit object. D1 first order derivative, D2 second order derivative, n curvature curvature.

Usage

D1(fit, t = NULL, analytical = FALSE, smoothed.spline = FALSE, ...)

D2(fit, t = NULL, analytical = FALSE, smoothed.spline = FALSE, ...)

## S3 method for class 'fFIT'
D1(fit, t = NULL, analytical = FALSE, smoothed.spline = FALSE, ...)

## S3 method for class 'fFIT'
D2(fit, t = NULL, analytical = FALSE, smoothed.spline = FALSE, ...)

curvature(fit, t = NULL, analytical = FALSE, smoothed.spline = FALSE, ...)

## S3 method for class 'fFIT'
curvature(fit, t = NULL, analytical = FALSE, smoothed.spline = FALSE, ...)

Arguments

Details

If fit$fun has no gradient function or smoothed.spline = TRUE, time-series smoothed by spline first, and get derivatives at last. If fit$fun exists and analytical = TRUE, smoothed.spline will be ignored.

Value

• der1 First order derivative

• der2 Second order derivative

• k Curvature

Examples

# doubleLog.Beck
t    <- seq(1, 365, 8)
tout <- seq(1, 365, 1)
par  = c(mn  = 0.1, mx  = 0.7, sos = 50, rsp = 0.1, eos = 250, rau = 0.1)
y <- doubleLog.Beck(par, t)

methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
fit <- curvefit(y, t, tout, methods)
x  <- fit$model$AG
d1 <- D1(x)
d2 <- D2(x)
d_k <- curvature(x)

Fine fitting

Description

Fine curve fitting function is used to fit vegetation time-series in every growing season.

Usage

FitDL.Zhang(y, t = index(y), tout = t, method = "nlm", w, type = 1L, ...)

FitDL.AG(y, t = index(y), tout = t, method = "nlminb", w, type = 1L, ...)

FitDL.AG2(y, t = index(y), tout = t, method = "nlminb", w, type = 1L, ...)

FitDL.Beck(y, t = index(y), tout = t, method = "nlminb", w, type = 1L, ...)

FitDL.Elmore(y, t = index(y), tout = t, method = "nlminb", w, type = 1L, ...)

FitDL.Gu(y, t = index(y), tout = t, method = "nlminb", w, type = 1L, ...)

FitDL.Klos(y, t = index(y), tout = t, method = "BFGS", w, type = 1L, ...)

Arguments

Value

• tout: The time of output curve fitting time-series.

• zs : Smoothed vegetation time-series of every iteration.

• ws : Weights of every iteration.

• par : Final optimized parameter of fine fitting.

• fun : The name of fine fitting.

References

1. Beck, P.S.A., Atzberger, C., Hogda, K.A., Johansen, B., Skidmore, A.K., 2006. Improved monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI. Remote Sens. Environ. https://doi.org/10.1016/j.rse.2005.10.021.

2. Elmore, A.J., Guinn, S.M., Minsley, B.J., Richardson, A.D., 2012. Landscape controls on the timing of spring, autumn, and growing season length in mid-Atlantic forests. Glob. Chang. Biol. 18, 656-674. https://doi.org/10.1111/j.1365-2486.2011.02521.x.
   

3. Gu, L., Post, W.M., Baldocchi, D.D., Black, TRUE.A., Suyker, A.E., Verma, S.B., Vesala, TRUE., Wofsy, S.C., 2009. Characterizing the Seasonal Dynamics of Plant Community Photosynthesis Across a Range of Vegetation Types, in: Noormets, A. (Ed.), Phenology of Ecosystem Processes: Applications in Global Change Research. Springer New York, New York, NY, pp. 35-58. https://doi.org/10.1007/978-1-4419-0026-5_2.
   

4. https://github.com/cran/phenopix/blob/master/R/FitDoubleLogGu.R

Examples

# simulate vegetation time-series
t    <- seq(1, 365, 8)
par  <- c(mn = 0.1, mx = 0.7, sos = 50, rsp = 0.1, eos = 250, rau = 0.1)
y    <- doubleLog.Beck(par, t)
data <- data.frame(t, y)
# methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang")
tout <- seq(1, 365, 1)
r <- FitDL.Elmore(y, t, tout)

plot(r, data)
get_GOF(r, data)
get_param(r)

GOF

Description

Good of fitting

Usage

GOF(Y_obs, Y_sim, w, include.r = TRUE, include.cv = FALSE)

Arguments

Value

• RMSE root mean square error

• NSE NASH coefficient

• MAE mean absolute error

• AI Agreement index (only good points (w == 1)) participate to calculate. See details in Zhang et al., (2015).

• Bias bias

• Bias_perc bias percentage

• n_sim number of valid obs

• cv Coefficient of variation

• R2 correlation of determination

• R pearson correlation

• pvalue pvalue of R

References

Zhang Xiaoyang (2015), http://dx.doi.org/10.1016/j.rse.2014.10.012

Examples

Y_obs = rnorm(100)
Y_sim = Y_obs + rnorm(100)/4
GOF(Y_obs, Y_sim)

Interface of unified optimization functions.

Description

optimx speed is not satisfied. So I_optim is present.

• I_optim: Interface of unified optimization functions.

• I_optimx: deprecated, which is about 10 times slower than I_optim.

Usage

I_optim(prior, FUN, y, t, method = "BFGS", ..., use.cpp = FALSE)

I_optimx(prior, FUN, y, t, method, verbose = FALSE, ..., use.cpp = FALSE)

Arguments

Value

• convcode: An integer code. 0 indicates successful convergence. Various methods may or may not return sufficient information to allow all the codes to be specified. An incomplete list of codes includes

  • 1: indicates that the iteration limit maxit had been reached.

  • 20: indicates that the initial set of parameters is inadmissible, that is, that the function cannot be computed or returns an infinite, NULL, or NA value.

  • 21: indicates that an intermediate set of parameters is inadmissible.

  • 10: indicates degeneracy of the Nelder–Mead simplex.

  • 51: indicates a warning from the "L-BFGS-B" method; see component message for further details.

  • 52: indicates an error from the "L-BFGS-B" method; see component message for further details.

  • 9999: error

• value: The value of fn corresponding to par

• par: The best parameter found

• nitns: the number of iterations

• fevals: The number of calls to objective.

See Also

stats::optim(), stats::nlminb(), stats::nlm(), optimx::optimx(), ucminf::ucminf()

Examples

# simulate vegetation time-series
FUN  = doubleLog_Beck
par  = c( mn  = 0.1 , mx  = 0.7 , sos = 50 , rsp = 0.1 , eos = 250, rau = 0.1)
par0 = c( mn  = 0.15, mx  = 0.65, sos = 100, rsp = 0.12, eos = 200, rau = 0.12)

t <- seq(1, 365, 8)
y <- FUN(par, t)

methods = c("BFGS", "ucminf", "nlm", "nlminb")
opt1 <- I_optim (par0, FUN, y, t, methods)
opt2 <- I_optimx(par0, FUN, y, t, methods)

# \dontrun{
# microbenchmark::microbenchmark(
#     opt1 = I_optim (par0, FUN, y, t, methods),
#     opt2 = I_optimx(par0, FUN, y, t, methods),
#     times = 2
# )
# }

Fine fitting functions

Description

double logistics, piecewise logistics and many other functions to curve fit VI time-series.

Usage

Logistic(par, t)

doubleLog.Zhang(par, t)

doubleLog.AG(par, t)

doubleLog.AG2(par, t)

doubleLog.Beck(par, t)

doubleLog.Elmore(par, t)

doubleLog.Gu(par, t)

doubleLog.Klos(par, t)

Arguments

Details

• Logistic The traditional simplest logistic function. It can be only used in half growing season, i.e. vegetation green-up or senescence period.

• doubleLog.Zhang Piecewise logistics, (Zhang Xiaoyang, RSE, 2003).

• doubleAG Asymmetric Gaussian.

• doubleLog.Beck Beck logistics.

• doubleLog.Gu Gu logistics.

• doubleLog.Elmore Elmore logistics.

• doubleLog.Klos Klos logistics.

All of those function have par and formula attributes for the convenience for analytical D1 and D2

References

1. Beck, P.S.A., Atzberger, C., Hogda, K.A., Johansen, B., Skidmore, A.K., 2006. Improved monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI. Remote Sens. Environ. https://doi.org/10.1016/j.rse.2005.10.021.

2. Elmore, A.J., Guinn, S.M., Minsley, B.J., Richardson, A.D., 2012. Landscape controls on the timing of spring, autumn, and growing season length in mid-Atlantic forests. Glob. Chang. Biol. 18, 656-674. https://doi.org/10.1111/j.1365-2486.2011.02521.x.
   

3. Gu, L., Post, W.M., Baldocchi, D.D., Black, TRUE.A., Suyker, A.E., Verma, S.B., Vesala, TRUE., Wofsy, S.C., 2009. Characterizing the Seasonal Dynamics of Plant Community Photosynthesis Across a Range of Vegetation Types, in: Noormets, A. (Ed.), Phenology of Ecosystem Processes: Applications in Global Change Research. Springer New York, New York, NY, pp. 35-58. https://doi.org/10.1007/978-1-4419-0026-5_2.
   

4. Peter M. Atkinson, et al., 2012, RSE, 123:400-417

5. https://github.com/cran/phenopix/blob/master/R/FitDoubleLogGu.R

Examples

# simulate vegetation time-series
t    <- seq(1, 365, 8)
par  <- c(mn = 0.1, mx = 0.7, sos = 50, rsp = 0.1, eos = 250, rau = 0.1)
y    <- doubleLog.Beck(par, t)
data <- data.frame(t, y)
# methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang")
tout <- seq(1, 365, 1)
r <- FitDL.Elmore(y, t, tout)

plot(r, data)
get_GOF(r, data)
get_param(r)

MOD13A1

Description

A data.table dataset, raw data of MOD13A1 data, clipped in 10 representative points ('DE-Obe', 'IT-Col', 'CN-Cha', 'AT-Neu', 'ZA-Kru', 'AU-How', 'CA-NS6', 'US-KS2', 'CH-Oe2', 'CZ-wet').

Usage

data('MOD13A1')

Format

An object of class list of length 2.

Details

Variables in MOD13A1:

• dt: vegetation index data

  • system:index: image index

  • DayOfYear: Numeric, Julian day of year

  • DayOfYear: corresponding doy of compositing NDVI and EVI

  • DetailedQA: VI quality indicators

  • SummaryQA: Quality reliability of VI pixel

  • EVI: Enhanced Vegetation Index

  • NDVI: Normalized Difference Vegetation Index

  • date: Date, corresponding date

  • site: String, site name

  • sur_refl_b01: Red surface reflectance

  • sur_refl_b02: NIR surface reflectance

  • sur_refl_b03: Blue surface reflectance

  • sur_refl_b07: MIR surface reflectance

  • .geo: geometry

• st: station info

  • ID: site ID

  • site: site name

  • lat: latitude

  • lon: longitude

  • IGBPname: IGBP land cover type

References

1. https://code.earthengine.google.com/dataset/MODIS/006/MOD13A1

Phenology extraction in Derivative method (DER)

Description

Phenology extraction in Derivative method (DER)

Usage

PhenoDeriv(x, t, ...)

## S3 method for class 'fFIT'
PhenoDeriv(x, t = NULL, analytical = FALSE, smoothed.spline = FALSE, ...)

## Default S3 method:
PhenoDeriv(x, t, der1, IsPlot = TRUE, show.legend = TRUE, ...)

Arguments

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

See Also

PhenoTrs(), PhenoGu(), PhenoKl()

Phenology extraction in GU method (GU)

Description

Phenology extraction in GU method (GU)

Usage

PhenoGu(x, t, ...)

## S3 method for class 'fFIT'
PhenoGu(x, t = NULL, analytical = FALSE, smoothed.spline = FALSE, ...)

## Default S3 method:
PhenoGu(x, t, der1, IsPlot = TRUE, ...)

Arguments

Value

A numeric vector, with the elements of:

• UD: upturn date

• SD: stabilisation date

• DD: downturn date

• RD: recession date

References

1. Gu, L., Post, W. M., Baldocchi, D. D., Black, T. A., Suyker, A. E., Verma, S. B., … Wofsy, S. C. (2009). Characterizing the Seasonal Dynamics of Plant Community Photosynthesis Across a Range of Vegetation Types. In A. Noormets (Ed.), Phenology of Ecosystem Processes: Applications in Global Change Research (pp. 35–58). New York, NY: Springer New York. doi:10.1007/978-1-4419-0026-5_2

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

Examples

# `doubleLog.Beck` simulate vegetation time-series
t    <- seq(1, 365, 8)
tout <- seq(1, 365, 1)
par  = c( mn  = 0.1 , mx  = 0.7 , sos = 50 , rsp = 0.1 , eos = 250, rau = 0.1)
y <- doubleLog.Beck(par, t)

methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
fit <- curvefit(y, t, tout, methods)
x <- fit$model$AG # one model

par(mfrow = c(2, 2))
PhenoTrs(x)
PhenoDeriv(x)
PhenoGu(x)
PhenoKl(x)

Phenology extraction in Inflection method (Zhang)

Description

Phenology extraction in Inflection method (Zhang)

Usage

PhenoKl(
  fFIT,
  t = NULL,
  analytical = FALSE,
  smoothed.spline = FALSE,
  IsPlot = TRUE,
  show.legend = TRUE,
  ...
)

Arguments

Value

A numeric vector, with the elements of: Greenup, Maturity, Senescence, Dormancy.

References

1. Zhang, X., Friedl, M. A., Schaaf, C. B., Strahler, A. H., Hodges, J. C. F. F., Gao, F., … Huete, A. (2003). Monitoring vegetation phenology using MODIS. Remote Sensing of Environment, 84(3), 471–475. doi:10.1016/S0034-4257(02)00135-9

Examples

# `doubleLog.Beck` simulate vegetation time-series
t    <- seq(1, 365, 8)
tout <- seq(1, 365, 1)
par  = c( mn  = 0.1 , mx  = 0.7 , sos = 50 , rsp = 0.1 , eos = 250, rau = 0.1)
y <- doubleLog.Beck(par, t)

methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
fit <- curvefit(y, t, tout, methods)
x <- fit$model$AG # one model

par(mfrow = c(2, 2))
PhenoTrs(x)
PhenoDeriv(x)
PhenoGu(x)
PhenoKl(x)

Phenology extraction in Threshold method (TRS)

Description

Phenology extraction in Threshold method (TRS)

Usage

PhenoTrs(
  x,
  t = NULL,
  approach = c("White", "Trs"),
  trs = 0.5,
  asymmetric = TRUE,
  IsPlot = TRUE,
  ...
)

## S3 method for class 'fFIT'
PhenoTrs(x, t = NULL, ...)

## Default S3 method:
PhenoTrs(
  x,
  t = NULL,
  approach = c("White", "Trs"),
  trs = 0.5,
  asymmetric = TRUE,
  IsPlot = TRUE,
  ...
)

Arguments

See Also

PhenoDeriv(), PhenoGu(), PhenoKl()

Examples

# `doubleLog.Beck` simulate vegetation time-series
t    <- seq(1, 365, 8)
tout <- seq(1, 365, 1)
par  = c( mn  = 0.1 , mx  = 0.7 , sos = 50 , rsp = 0.1 , eos = 250, rau = 0.1)
y <- doubleLog.Beck(par, t)

methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
fit <- curvefit(y, t, tout, methods)
x <- fit$model$AG # one model

par(mfrow = c(2, 2))
PhenoTrs(x)
PhenoDeriv(x)
PhenoGu(x)
PhenoKl(x)

Critical value of determined correlation

Description

Critical value of determined correlation

Usage

R2_sign(n, NumberOfPredictor = 2, alpha = 0.05)

Arguments

Value

F statistic and R2 at significant level.

References

Chen Yanguang (2012), Geographical Data analysis with MATLAB.

Examples

R2_critical <- R2_sign(30, NumberOfPredictor = 2, alpha = 0.05)

Add one year data in the head and tail

Description

Add the data of the year of year_start - 1 to the head, add the data of the year of year_end - 1 to the tail.

Usage

add_HeadTail(d, south = FALSE, nptperyear, trs = 0.45)

Arguments

Value

data.table

Note

date is image date; t is compositing date.

Examples

library(phenofit)
data("CA_NS6"); d = CA_NS6

nptperyear = 23
dnew     <- add_HeadTail(d, nptperyear = nptperyear) # add one year in head and tail

get rough fitting

Description

get rough fitting

Usage

brks2rfit(brks)

Arguments

Value

• data:

  • t

  • y

  • QC_flag

• tout:

• zs: list of iter1, ..., itern

• ws: list of iter1, ..., itern

Check growing season head and tail minimum values

Description

Check growing season head and tail minimum values

Usage

check_GS_HeadTail(pos, ypred, minlen, A = NULL, ...)

Arguments

Value

di

check_input

Description

Check input data, interpolate NA values in y, remove spike values, and set weights for NA in y and w.

Usage

check_input(
  t,
  y,
  w,
  QC_flag,
  nptperyear,
  south = FALSE,
  wmin = 0.2,
  wsnow = 0.8,
  ymin,
  missval,
  maxgap,
  alpha = 0.02,
  alpha_high = NULL,
  date_start = NULL,
  date_end = NULL,
  mask_spike = TRUE,
  na.rm = FALSE,
  ...
)

Arguments

Value

A list object returned:

• t : Numeric vector

• y0: Numeric vector, original vegetation time-series.

• y : Numeric vector, checked vegetation time-series, NA values are interpolated.

• w : Numeric vector

• Tn: Numeric vector

• ylu: = ⁠[ymin, ymax]⁠. w_critical is used to filter not too bad values.

  If the percentage good values (w=1) is greater than 30\

  The else, if the percentage of w >= 0.5 points is greater than 10\ w_critical=0.5. In boreal regions, even if the percentage of w >= 0.5 points is only 10\

  We can't rely on points with the wmin weights. Then,
  y_good = y[w >= w_critical],
  ymin = pmax( quantile(y_good, alpha/2), 0)
ymax = max(y_good).

Examples

data("CA_NS6")
d = CA_NS6
# head(d)

nptperyear = 23
INPUT <- check_input(d$t, d$y, d$w, QC_flag = d$QC_flag,
     nptperyear = nptperyear, south = FALSE, 
     maxgap = nptperyear/4, alpha = 0.02, wmin = 0.2)
plot_input(INPUT)

check_ylu

Description

Curve fitting values are constrained in the range of ylu. Only constrain trough value for a stable background value. But not for peak value.

Usage

check_ylu(yfit, ylu)

Arguments

Value

yfit, the numeric vector in the range of ylu.

Examples

check_ylu(1:10, c(2, 8))

Fine curve fitting

Description

Curve fit vegetation index (VI) time-series of every growing season using fine curve fitting methods.

Usage

curvefit(
  y,
  t = index(y),
  tout = t,
  methods = c("AG", "Beck", "Elmore", "Gu", "Klos", "Zhang"),
  w = NULL,
  ...,
  type = 1L,
  use.cpp = FALSE
)

Arguments

Value

fFITs S3 object, see fFITs() for details.

Note

'Klos' have too many parameters. It will be slow and not stable.

See Also

fFITs()

Examples

library(phenofit)
# simulate vegetation time-series
FUN = 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 <- FUN(par, t)

methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
fit <- curvefit(y, t, tout = tout, methods)

curvefit0

Description

curvefit0

Usage

curvefit0(
  y,
  t = index(y),
  tout = t,
  methods = c("AG", "Beck", "Elmore", "Gu", "Klos", "Zhang"),
  w = NULL,
  ...
)

Fine Curve fitting

Description

Fine Curve fitting for INPUT time-series.

Usage

curvefits(INPUT, brks, options = list(), ...)

Arguments

Value

List of phenofit fitting object.

options for fitting

• methods (default ⁠c('AG', 'Beck', 'Elmore', 'Zhang')``): Fine curve fitting methods, can be one or more of ⁠c('AG', 'Beck', 'Elmore', 'Zhang', 'Gu', 'Klos')‘. Note that ’Gu' and 'Klos' are very slow.

• iters (default 2): max iterations of fine fitting.

• wFUN (default wTSM): Character or function, weights updating function of fine fitting function.

• wmin (default 0.1): min weights in the weights updating procedure.

• use.rough (default FALSE): Whether to use rough fitting smoothed time-series as input? If false, smoothed VI by rough fitting will be used for Phenological metrics extraction; If true, original input y will be used (rough fitting is used to divide growing seasons and update weights.

• use.y0 (default TRUE): boolean. whether to use original y0 as the input of plot_input, note that not for curve fitting. y0 is the original value before the process of check_input.

• nextend (default 2): Extend curve fitting window, until nextend good or marginal points are found in the previous and subsequent growing season.

• maxExtendMonth (default 1): Search good or marginal good values in previous and subsequent maxExtendMonth period.

• minExtendMonth (default 0.5): Extend period defined by nextend and maxExtendMonth, should be no shorter than minExtendMonth. When all points of the input time-series are good value, then the extending period will be too short. In that situation, we can't make sure the connection between different growing seasons is smoothing.

• minPercValid: (default 0, not use). If the percentage of good- and marginal- quality points is less than minPercValid, curve fiting result is set to NA.

• minT: (not use). If Tn not provided in INPUT, minT will not be used. minT use night temperature Tn to define backgroud value (days with Tn < minT treated as ungrowing season).

See Also

FitDL()

Examples

data("CA_NS6")
d = CA_NS6

nptperyear <- 23
INPUT <- check_input(d$t, d$y, d$w, QC_flag = d$QC_flag,
     nptperyear = nptperyear, south = FALSE,
     maxgap = nptperyear/4, alpha = 0.02, wmin = 0.2)
# plot_input(INPUT)

# Rough fitting and growing season dividing
wFUN <- "wTSM"
brks2 <- season_mov(INPUT,
    options = list(
        rFUN = "smooth_wWHIT", wFUN = wFUN,
        r_min = 0.05, ypeak_min = 0.05,
        lambda = 10,
        verbose = FALSE
    ))
# plot_season(INPUT, brks2, d)
# Fine fitting
fits <- curvefits(
    INPUT, brks2,
    options = list(
        methods = c("AG", "Beck", "Elmore", "Zhang"), #,"klos", "Gu"
        wFUN = wFUN,
        nextend = 2, maxExtendMonth = 2, minExtendMonth = 1, minPercValid = 0.2
    )
)

r_param = get_param(fits)
r_pheno = get_pheno(fits)
r_gof = get_GOF(fits)
d_fit = get_fitting(fits)

g <- plot_curvefits(d_fit, brks2)
grid::grid.newpage(); grid::grid.draw(g)

curvefits by local model functions of TIMESAT

Description

Local model functions f_L(t), f_C(t) and f_R(t) describe the VI variation in intervals around the left minima, the central maxima and the right minima.

Local model function are merged into global model function via merge_LocalModels() and Per J\"onsson et al. (2004; their Eq. 12), where cut-off function sharply drop from 1 to 0 in small intervals around (t_L + t_C)/2 and (t_C + t_R)/2.

F(t)= \begin{cases} \alpha(t) f_{L}(t)+[1-\alpha(t)] f_{C}(t), t_{L}<t<t_{C} \\ \beta(t) f_{C}(t)+[1-\beta(t)] f_{R}(t), t_{C}<t<t_{R}\end{cases}

Usage

curvefits_LocalModel(INPUT, brks, options = list(), ...)

merge_LocalModels(fits)

Arguments

options for fitting

• methods (default ⁠c('AG', 'Beck', 'Elmore', 'Zhang')``): Fine curve fitting methods, can be one or more of ⁠c('AG', 'Beck', 'Elmore', 'Zhang', 'Gu', 'Klos')‘. Note that ’Gu' and 'Klos' are very slow.

• iters (default 2): max iterations of fine fitting.

• wFUN (default wTSM): Character or function, weights updating function of fine fitting function.

• wmin (default 0.1): min weights in the weights updating procedure.

• use.rough (default FALSE): Whether to use rough fitting smoothed time-series as input? If false, smoothed VI by rough fitting will be used for Phenological metrics extraction; If true, original input y will be used (rough fitting is used to divide growing seasons and update weights.

• use.y0 (default TRUE): boolean. whether to use original y0 as the input of plot_input, note that not for curve fitting. y0 is the original value before the process of check_input.

• nextend (default 2): Extend curve fitting window, until nextend good or marginal points are found in the previous and subsequent growing season.

• maxExtendMonth (default 1): Search good or marginal good values in previous and subsequent maxExtendMonth period.

• minExtendMonth (default 0.5): Extend period defined by nextend and maxExtendMonth, should be no shorter than minExtendMonth. When all points of the input time-series are good value, then the extending period will be too short. In that situation, we can't make sure the connection between different growing seasons is smoothing.

• minPercValid: (default 0, not use). If the percentage of good- and marginal- quality points is less than minPercValid, curve fiting result is set to NA.

• minT: (not use). If Tn not provided in INPUT, minT will not be used. minT use night temperature Tn to define backgroud value (days with Tn < minT treated as ungrowing season).

References

1. Per J\"onsson, P., Eklundh, L., 2004. TIMESAT - A program for analyzing time-series of satellite sensor data. Comput. Geosci. 30, 833-845. doi:10.1016/j.cageo.2004.05.006.

See Also

curvefits()

Examples

## Not run: 
library(phenofit)

data("CA_NS6")
d = CA_NS6

nptperyear <- 23
INPUT <- check_input(d$t, d$y, d$w, QC_flag = d$QC_flag,
     nptperyear = nptperyear, south = FALSE,
     maxgap = nptperyear/4, alpha = 0.02, wmin = 0.2)
# plot_input(INPUT)

# Rough fitting and growing season dividing
wFUN <- "wTSM"
brks2 <- season_mov(INPUT,
    options = list(
        rFUN = "smooth_wWHIT", wFUN = wFUN,
        r_min = 0.05, ypeak_min = 0.05,
        lambda = 10,
        verbose = FALSE
    ))
# plot_season(INPUT, brks2, d)

# Fine fitting
fits <- curvefits_LocalModel(
    INPUT, brks2,
    options = list(
        methods = c("AG", "Beck", "Elmore", "Zhang", "Gu"), #,"klos", "Gu"
        wFUN = wFUN,
        nextend = 2, maxExtendMonth = 2, minExtendMonth = 1, minPercValid = 0.2
    ),
    constrain = TRUE
)
# merge local model function into global model function
fits_merged = merge_LocalModels(fits) 

## Visualization ---------------------------------------------------------------
l_fitting = map(fits %>% guess_names, get_fitting) #%>% melt_list("period")

d_merged = get_fitting(fits_merged[[2]]) %>% cbind(type = "Merged")
d_raw = l_fitting[2:4] %>% set_names(c("Left", "Central", "Right")) %>%
    melt_list("type")
d_obs = d_raw[, .(t, y, QC_flag)] %>% unique()
d_fit = rbind(d_merged, d_raw)[meth == "Zhang"]

levs = c("Left", "Central", "Right", "Merged")
levs_new = glue("({letters[1:4]}) {levs}") %>% as.character()
d_fit$type %<>% factor(levs, levs_new)

p = ggplot(d_obs, aes(t, y)) +
    geom_point() +
    geom_line(data = d_fit, aes(t, ziter2, color = type)) +
    facet_wrap(~type) +
    labs(x = "Date", y = "EVI") +
    scale_x_date(date_labels = "%b %Y", expand = c(1, 1)*0.08) +
    theme_bw(base_size = 13) +
    theme(legend.position = "none",
          strip.text = element_text(size = 14))
p

## End(Not run)

cutoff

Description

cutoff

Usage

cutoff(n1, n2, k = n1:n2)

Arguments

Value

A numeric vector, in small intervals around (tL + tC)/2 and (tC+tR)/2, respectively, smoothly drop from 1 to 0.

References

1. Per J\"onsson, P., Eklundh, L., 2004. TIMESAT - A program for analyzing time-series of satellite sensor data. Comput. Geosci. 30, 833-845. https://doi.org/10.1016/j.cageo.2004.05.006.

weighted CV

Description

weighted CV

Usage

cv_coef(x, w)

Arguments

Value

Named numeric vector, (mean, sd, cv).

Examples

library(phenofit)
x = rnorm(100)
coefs <- cv_coef(x)

S3 class of fine curve fitting object.

Description

fFIT is returned by optim_pheno().

Format

• tout: Corresponding doy of prediction

• zs: curve fitting values of every iteration

• ws: weight of every iteration

• par: Optimized parameter of fine curve fitting method

• fun: The name of fine curve fitting function.

S3 class of multiple fine curve fittings object.

Description

plot curve fitting VI, gradient (first order difference D1), hessian (D2), curvature (k) and the change rate of curvature(der.k)

Usage

## S3 method for class 'fFITs'
plot(x, method, show.pheno = TRUE, ...)

## S3 method for class 'fFIT'
plot(x, data, ...)

Arguments

See Also

curvature()

Examples

library(phenofit)
# 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)

methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
fit <- curvefit(y, t, tout, methods)

# plot
plot(fit)

Goal function of fine curve fitting methods

Description

Goal function of fine curve fitting methods

Usage

f_goal(par, fun, y, t, pred, w, ylu, ...)

Arguments

Value

RMSE Root Mean Square Error of curve fitting values.

Examples

library(phenofit)

par  = c( mn  = 0.1 , mx  = 0.7 , sos = 50 , rsp = 0.1 , eos = 250, rau = 0.1)
par0 = c( mn  = 0.15, mx  = 0.65, sos = 100, rsp = 0.12, eos = 200, rau = 0.12)

# simulate vegetation time-series
fFUN = doubleLog_Beck
t    <- seq(1, 365, 8)
tout <- seq(1, 365, 1)
y    <- fFUN(par, t)

f_goal(par0, fFUN, y, t)

objective function of double logistics

Description

objective function of double logistics

Usage

f_goal2(par, fun, y, t, pred, w = NULL, ylu = NULL, ...)

Arguments

find_season

Description

find_season

Usage

find_season.peaks(rfit, info_peak, options = list(), ...)

find_season.default(
  ypred,
  t = seq_along(ypred),
  nptperyear = NULL,
  south = NULL,
  options = list(),
  ...
)

Arguments

findpeaks

Description

Find peaks (maxima) in a time series. This function is modified from pracma::findpeaks.

Usage

findpeaks(
  x,
  nups = 1,
  ndowns = nups,
  zero = "0",
  peakpat = NULL,
  minpeakheight = -Inf,
  minpeakdistance = 1,
  h_min = 0,
  h_max = 0,
  npeaks = 0,
  sortstr = FALSE,
  include_gregexpr = FALSE,
  IsPlot = F
)

Arguments

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)
hgt <- c(4, 5, 3, 4, 5, 4.2, 2.1, 4.3, 3.1, 5.1, 4.2)
wdt <- c(0.005, 0.005, 0.006, 0.01, 0.01, 0.03, 0.01, 0.01, 0.005, 0.008, 0.005)
pSignal <- numeric(length(x))
for (i in seq(along=pos)) {
    pSignal <- pSignal + hgt[i]/(1 + abs((x - pos[i])/wdt[i]))^4
}

plot(pSignal, type="l", col="navy"); grid()
x <- findpeaks(pSignal, npeaks=3, h_min=4, sortstr=TRUE)
points(val~pos, x$X, pch=20, col="maroon")

getRealDate

Description

convert MODIS DayOfYear to the exact compositing date.

Usage

getRealDate(date, DayOfYear)

Arguments

Value

A data.table with a new column t, which is the exact compositing date.

Examples

library(phenofit)
data("MOD13A1")

df  <- MOD13A1$dt
df$t <- getRealDate(df$date, df$DayOfYear)

get_GOF

Description

Goodness-of-fitting (GOF) of fine curve fitting results.

Usage

get_GOF(x, ...)

## S3 method for class 'list'
get_GOF(x, ...)

## S3 method for class 'fFITs'
get_GOF(x, ...)

## S3 method for class 'fFIT'
get_GOF(x, data, ...)

Arguments

Value

• meth: The name of fine curve fitting method

• RMSE: Root Mean Square Error

• NSE : Nash-Sutcliffe model efficiency coefficient

• R : Pearson-Correlation

• R2 : determined coefficient

• pvalue: pvalue of R

• n : The number of observations

References

1. https://en.wikipedia.org/wiki/Nash-Sutcliffe_model_efficiency_coefficient
   

2. https://en.wikipedia.org/wiki/Pearson_correlation_coefficient

See Also

curvefit()

Examples

library(phenofit)
# simulate vegetation time-series
FUN = 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 <- FUN(par, t)
methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
fit <- curvefit(y, t, tout, methods) # `fFITs` (fine-fitting) object 
fits <- list(`2001` = fit, `2002` = fit) # multiple years

l_param   <- get_param(fits)
d_GOF     <- get_GOF(fits)
d_fitting <- get_fitting(fits)
l_pheno   <- get_pheno(fits, "AG", IsPlot=TRUE)

getFittings

Description

Get curve fitting data.frame

Usage

get_fitting(x)

## S3 method for class 'list'
get_fitting(x)

## S3 method for class 'fFITs'
get_fitting(x)

Arguments

Examples

library(phenofit)
# simulate vegetation time-series
FUN = 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 <- FUN(par, t)
methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
fit <- curvefit(y, t, tout, methods) # `fFITs` (fine-fitting) object 
fits <- list(`2001` = fit, `2002` = fit) # multiple years

l_param   <- get_param(fits)
d_GOF     <- get_GOF(fits)
d_fitting <- get_fitting(fits)
l_pheno   <- get_pheno(fits, "AG", IsPlot=TRUE)

Get parameters from curve fitting result

Description

Get parameters from curve fitting result

Usage

get_param(x)

## S3 method for class 'list'
get_param(x)

## S3 method for class 'fFITs'
get_param(x)

## S3 method for class 'fFIT'
get_param(x)

Arguments

Value

A list of tibble with the length being equal to the number of methods. Each line of tibble cotains the corresponding parameters of each growing season.

Examples

library(phenofit)
# simulate vegetation time-series
FUN = 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 <- FUN(par, t)
methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
fit <- curvefit(y, t, tout, methods) # `fFITs` (fine-fitting) object 
fits <- list(`2001` = fit, `2002` = fit) # multiple years

l_param   <- get_param(fits)
d_GOF     <- get_GOF(fits)
d_fitting <- get_fitting(fits)
l_pheno   <- get_pheno(fits, "AG", IsPlot=TRUE)

get_pheno

Description

Get yearly vegetation phenological metrics of a curve fitting method

Usage

get_pheno(x, ...)

## S3 method for class 'rfit'
get_pheno(x, TRS = c(0.2, 0.5), asymmetric = TRUE, ...)

## S3 method for class 'list'
get_pheno(
  x,
  method,
  TRS = c(0.2, 0.5, 0.6),
  analytical = FALSE,
  smoothed.spline = FALSE,
  IsPlot = FALSE,
  show.title = TRUE,
  ...
)

## S3 method for class 'fFITs'
get_pheno(
  x,
  method,
  TRS = c(0.2, 0.5),
  analytical = FALSE,
  smoothed.spline = FALSE,
  IsPlot = FALSE,
  title.left = "",
  show.PhenoName = TRUE,
  ...
)

Arguments

Value

List of every year phenology metrics

Examples

library(phenofit)
# simulate vegetation time-series
FUN = 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 <- FUN(par, t)
methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
fit <- curvefit(y, t, tout, methods) # `fFITs` (fine-fitting) object 
fits <- list(`2001` = fit, `2002` = fit) # multiple years

l_param   <- get_param(fits)
d_GOF     <- get_GOF(fits)
d_fitting <- get_fitting(fits)
l_pheno   <- get_pheno(fits, "AG", IsPlot=TRUE)

Initial lambda value of Whittaker smoother

Description

This function is only suitable for 16-day EVI time-series.

Usage

init_lambda(y)

Arguments

Examples

library(phenofit)
data("MOD13A1")

dt <- tidy_MOD13(MOD13A1$dt)
st <- MOD13A1$st

sitename <- dt$site[1]
d      <- dt[site == sitename, ] # get the first site data
lambda <- init_lambda(d$y)

init_param

Description

Initialize parameters of double logistic function

Usage

init_param(y, t, w, type = 1L)

init_Zhang(e, type = 1L, ...)

init_AG(e, type = 1L, ...)

init_AG2(e, type = 1L, ...)

init_Beck(e, type = 1L, ...)

init_Elmore(e, type = 1L, ...)

init_Gu(e, type = 1L, ...)

init_Klos(e, type = 1L, ...)

Arguments

Examples

library(phenofit)
# 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)

l_param <- init_param(y, t)

input object with one growing season per year

Description

Variables in input_single:

• t: date of compositing image

• y: EVI

• w: weights of data point

• ylu: lower and upper boundary

• nptperyear: points per year

• south: boolean, whether in south Hemisphere?

Usage

data('input_single')

Format

An object of class list of length 6.

skewness and kurtosis

Description

Inherit from package e1071

Usage

kurtosis(x, na.rm = FALSE, type = 3)

skewness(x, na.rm = FALSE, type = 3)

Arguments

Examples

x = rnorm(100)
coef_kurtosis <- kurtosis(x)
coef_skewness <- skewness(x)

lambda_vcurve

Description

lambda_vcurve

Usage

lambda_vcurve(
  y,
  w,
  lg_lambdas = seq(0.1, 5, 0.1),
  plot = FALSE,
  verbose = FALSE,
  ...
)

Arguments

Double logistics in Rcpp

Description

Double logistics in Rcpp

Usage

logistic(par, t, pred)

doubleLog_Zhang(par, t, pred)

doubleLog_AG(par, t, pred)

doubleLog_Beck(par, t, pred)

doubleLog_Elmore(par, t, pred)

doubleLog_Gu(par, t, pred)

doubleLog_Klos(par, t, pred)

Arguments

See Also

doubleLog.Beck()

melt_list

Description

melt_list

Usage

melt_list(list, var.name = "variable", na.rm = TRUE, ...)

movmean

Description

NA and Inf values in the yy will be ignored automatically.

Usage

movmean(y, halfwin = 1L, SG_style = FALSE, w = NULL)

Arguments

Examples

x <- 1:100
x[50] <- NA; x[80] <- Inf
s1 <- movmean(x, 2, SG_style = TRUE)
s2 <- movmean(x, 2, SG_style = FALSE)

Unified optimization function

Description

I_optimx is rich of functionality, but with a low computing performance. Some basic optimization functions are unified here, with some input and output format.

• opt_ncminf General-Purpose Unconstrained Non-Linear Optimization, see ucminf::ucminf().

• opt_nlminb Optimization using PORT routines, see stats::nlminb().

• opt_nlm Non-Linear Minimization, stats::nlm().

• opt_optim General-purpose Optimization, see stats::optim().

Usage

opt_ucminf(par0, objective, ...)

opt_nlm(par0, objective, ...)

opt_optim(par0, objective, method = "BFGS", ...)

opt_nlminb(par0, objective, ...)

Arguments

Value

• convcode: An integer code. 0 indicates successful convergence. Various methods may or may not return sufficient information to allow all the codes to be specified. An incomplete list of codes includes

  • 1: indicates that the iteration limit maxit had been reached.

  • 20: indicates that the initial set of parameters is inadmissible, that is, that the function cannot be computed or returns an infinite, NULL, or NA value.

  • 21: indicates that an intermediate set of parameters is inadmissible.

  • 10: indicates degeneracy of the Nelder–Mead simplex.

  • 51: indicates a warning from the "L-BFGS-B" method; see component message for further details.

  • 52: indicates an error from the "L-BFGS-B" method; see component message for further details.

  • 9999: error

• value: The value of fn corresponding to par

• par: The best parameter found

• nitns: the number of iterations

• fevals: The number of calls to objective.

See Also

optim_pheno(), I_optim()

Examples

library(phenofit)
library(ggplot2)
library(magrittr)
library(purrr)
library(data.table)

# 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)
par0 = c( mn  = 0.15, mx  = 0.65, sos = 100, rsp = 0.12, eos = 200, rau = 0.12)

t    <- seq(1, 365, 8)
tout <- seq(1, 365, 1)
y <- fFUN(par, t)

optFUNs <- c("opt_ucminf", "opt_nlminb", "opt_nlm", "opt_optim") %>% set_names(., .)
opts <- lapply(optFUNs, function(optFUN){
    optFUN <- get(optFUN)
    opt    <- optFUN(par0, f_goal, y = y, t = t, fun = fFUN)
    opt$ysim <- fFUN(opt$par, t)
    opt
})

# visualization
df   <- map(opts, "ysim") %>% as.data.table() %>% cbind(t, y, .)
pdat <- data.table::melt(df, c("t", "y"), variable.name = "optFUN")

ggplot(pdat) + 
    geom_point(data = data.frame(t, y), aes(t, y), size = 2) + 
    geom_line(aes(t, value, color = optFUN), linewidth = 0.9)

optim_pheno

Description

Interface of optimization functions for double logistics and other parametric curve fitting functions.

Usage

optim_pheno(
  prior,
  sFUN,
  y,
  t,
  tout,
  method,
  w,
  nptperyear,
  ylu,
  iters = 2,
  wFUN = wTSM,
  lower = -Inf,
  upper = Inf,
  constrain = TRUE,
  verbose = FALSE,
  ...,
  use.cpp = FALSE
)

Arguments

Value

A fFIT() object, with the element of:

• tout: The time of output curve fitting time-series.

• zs : Smoothed vegetation time-series of every iteration.

• ws : Weights of every iteration.

• par : Final optimized parameter of fine fitting.

• fun : The name of fine fitting.

See Also

fFIT(), stats::nlminb()

Examples

# library(magrittr)
# library(purrr)

# simulate vegetation time-series
t    <- seq(1, 365, 8)
tout <- seq(1, 365, 1)

FUN = doubleLog_Beck
par  = c( mn  = 0.1 , mx  = 0.7 , sos = 50 , rsp = 0.1 , eos = 250, rau = 0.1)
par0 = c( mn  = 0.15, mx  = 0.65, sos = 100, rsp = 0.12, eos = 200, rau = 0.12)

y <- FUN(par, t)

methods = c("BFGS", "ucminf", "nlm", "nlminb")
opt1 <- I_optim(par0, doubleLog_Beck, y, t, methods) # "BFGS", "ucminf", "nlm",
# opt2 <- I_optimx(prior, fFUN, y, t, tout, )

sFUN   = "doubleLog.Beck" # doubleLog.Beck
r <- optim_pheno(par0, sFUN, y, t, tout, method = methods[4],
                 nptperyear = 46, iters = 2, wFUN = wTSM, verbose = FALSE, use.julia = FALSE)

plot_curvefits

Description

plot_curvefits

Usage

plot_curvefits(
  d_fit,
  seasons,
  d_obs = NULL,
  title = NULL,
  xlab = "Time",
  ylab = "Vegetation Index",
  yticks = NULL,
  font.size = 14,
  theme = NULL,
  cex = 2,
  shape = "point",
  angle = 30,
  show.legend = TRUE,
  layer_extra = NULL,
  ...
)

Arguments

Examples

data("CA_NS6")
d = CA_NS6

nptperyear <- 23
INPUT <- check_input(d$t, d$y, d$w, QC_flag = d$QC_flag,
     nptperyear = nptperyear, south = FALSE,
     maxgap = nptperyear/4, alpha = 0.02, wmin = 0.2)
# plot_input(INPUT)

# Rough fitting and growing season dividing
wFUN <- "wTSM"
brks2 <- season_mov(INPUT,
    options = list(
        rFUN = "smooth_wWHIT", wFUN = wFUN,
        r_min = 0.05, ypeak_min = 0.05,
        lambda = 10,
        verbose = FALSE
    ))
# plot_season(INPUT, brks2, d)
# Fine fitting
fits <- curvefits(
    INPUT, brks2,
    options = list(
        methods = c("AG", "Beck", "Elmore", "Zhang"), #,"klos", "Gu"
        wFUN = wFUN,
        nextend = 2, maxExtendMonth = 2, minExtendMonth = 1, minPercValid = 0.2
    )
)

r_param = get_param(fits)
r_pheno = get_pheno(fits)
r_gof = get_GOF(fits)
d_fit = get_fitting(fits)

g <- plot_curvefits(d_fit, brks2)
grid::grid.newpage(); grid::grid.draw(g)

Plot INPUT returned by check_input

Description

Plot INPUT returned by check_input

Usage

plot_input(INPUT, wmin = 0.2, show.y0 = TRUE, ylab = "VI", ...)

Arguments

Examples

library(phenofit)
data("CA_NS6"); d = CA_NS6
# global parameter
IsPlot = TRUE
nptperyear = 23
ypeak_min  = 0.05

INPUT    <- check_input(d$t, d$y, d$w, d$QC_flag, nptperyear,
                        maxgap = nptperyear/4, alpha = 0.02, wmin = 0.2)
plot_input(INPUT)

plot_phenofit

Description

plot_phenofit

Usage

plot_phenofit(
  obj,
  type = "all",
  methods,
  title = NULL,
  ylab = "Vegetation Index",
  IsPlot = TRUE,
  show.legend = TRUE,
  newpage = TRUE,
  ...
)

Arguments

plot_season

Description

Plot growing season divding result.

Usage

plot_season(
  INPUT,
  brks,
  plotdat,
  IsPlot.OnlyBad = FALSE,
  show.legend = TRUE,
  ylab = "VI",
  title = NULL,
  show.shade = TRUE,
  margin = 0.35
)

Arguments

Extract Vegetation Phenology at site scale

Description

Extract Vegetation Phenology at site scale

Usage

process_phenofit(
  d_obs,
  nptperyear = 36,
  south = FALSE,
  options_season = list(wFUN = "wTSM", maxExtendMonth = 12, MaxPeaksPerYear = 3,
    MaxTroughsPerYear = 4),
  options_fitting = list(methods = c("AG", "Zhang", "Beck", "Elmore", "Gu"), wFUN =
    "wTSM", maxExtendMonth = 12, minExtendMonth = 0.5, use.y0 = FALSE),
  brks = NULL,
  TRS = c(0.1, 0.2, 0.5, 0.6, 0.8, 0.9),
  run.curvefit = TRUE,
  ...
)

Arguments

divide_seasons

Description

divide_seasons

Usage

process_season(
  d_obs,
  options = list(wFUN = "wTSM", maxExtendMonth = 12, MaxPeaksPerYear = 3,
    MaxTroughsPerYear = 4),
  nptperyear = 36,
  south = FALSE,
  ...
)

Arguments

Note

site-year may be not continuous.

Initial weights according to qc

Description

• getBits: Extract bitcoded QA information from bin value

• qc_summary: Initial weigths based on Quality reliability of VI pixel, suit for MOD13A1, MOD13A2 and MOD13Q1 (SummaryQA band).

• qc_5l: Initial weights based on Quality control of five-level confidence score, suit for MCD15A3H(LAI, FparLai_QC), MOD17A2H(GPP, Psn_QC) and MOD16A2(ET, ET_QC).

• qc_StateQA: Initial weights based on StateQA, suit for MOD09A1, MYD09A1.

• qc_FparLai: For MODIS LAI

• qc_NDVI3g: For AVHRR NDVI3g

• qc_NDVIv4: For AVHRR NDVIv4

Usage

getBits(x, start, end = start)

qc_summary(QA, wmin = 0.2, wmid = 0.5, wmax = 1)

qc_StateQA(QA, wmin = 0.2, wmid = 0.5, wmax = 1)

qc_FparLai(QA, FparLai_QC = NULL, wmin = 0.2, wmid = 0.5, wmax = 1)

qc_5l(QA, wmin = 0.2, wmid = 0.5, wmax = 1)

qc_NDVIv4(QA, wmin = 0.2, wmid = 0.5, wmax = 1)

qc_NDVI3g(QA, wmin = 0.2, wmid = 0.5, wmax = 1)

qc_SPOT(QA, wmin = 0.2, wmid = 0.5, wmax = 1)

Arguments

Details

If FparLai_QC specified, I_margin = SCF_QC >= 2 & SCF_QC <= 3.

Value

A list object with

• weigths: Double vector, initial weights.

• QC_flag: Factor vector, with the level of c("snow", "cloud", "shadow", "aerosol", "marginal", "good")

Note

qc_5l and qc_NDVIv4 only returns weight, without QC_flag.

References

https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD13A1

https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MCD15A3H

Erwin Wolters, Else Swinnen, Carolien Toté, Sindy Sterckx. SPOT-VGT COLLECTION 3 PRODUCTS USER MANUAL V1.2, 2018, P47

See Also

qc_sentinel2()

Examples

set.seed(100)
QA <- as.integer(runif(100, 0, 2^7))

r1 <- qc_summary(QA, wmin = 0.2, wmid = 0.5, wmax = 1)
r2 <- qc_StateQA(QA, wmin = 0.2, wmid = 0.5, wmax = 1)
r_5l <- qc_5l(QA, wmin = 0.2, wmid = 0.5, wmax = 1)
r_NDVI3g <- qc_NDVI3g(QA, wmin = 0.2, wmid = 0.5, wmax = 1)
r_NDVIv4 <- qc_NDVIv4(QA, wmin = 0.2, wmid = 0.5, wmax = 1)

qc level, color and shape

Description

qc level, color and shape

Usage

qc_levels

qc_colors

qc_shapes

Format

An object of class character of length 6.

An object of class character of length 6.

An object of class numeric of length 6.

Initial weights for sentinel2 according to SCL band

Description

Usage

qc_sentinel2(SCL, wmin = 0.2, wmid = 0.5, wmax = 1)

Arguments

References

https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR

Examples

qc_sentinel2(1:11)

season_filter

Description

season_filter

Usage

rcpp_season_filter(d, rm_closed = TRUE, rtrough_max = 0.7, r_min = 0.02)

check_season_dt(dt)

check_season_list(lst_dt)

Arguments

Weighted Savitzky-Golay written in RcppArmadillo

Description

NA and Inf values in the yy has been ignored automatically.

Usage

rcpp_wSG(y, halfwin = 1L, d = 1L, w = NULL)

rcpp_SG(y, halfwin = 1L, d = 1L)

Arguments

Examples

y <- 1:15
w <- seq_along(y)/length(y)

frame = 5
d = 2
s1 <- rcpp_wSG(y, frame, d, w)
s2 <- rcpp_SG(y, frame, d)

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

lubridate

make_date

magrittr

%<>%, %>%

zeallot

%<-%

rm too closed peaks or troughs

Description

rm too closed peaks or troughs

Usage

removeClosedExtreme(pos, ypred, A = NULL, y_min)

Arguments

Value

A data.table with the columns of: ⁠"val", "pos", "left", "right", "type"⁠

Examples

#removeClosedExtreme(pos, ypred, A, y_min)

Rough fitting

Description

Divide growing seasons according to rough fitting (rFUN) result .

For season, rough fitting is applied for whole. For season_mov rough fitting is applied in every year, during which maxExtendMonth is extended.

Usage

roughFit(
  INPUT,
  options = list(),
  frame = floor(INPUT$nptperyear/5) * 2 + 1,
  ...
)

Arguments

Details

Before division growing season, INPUT should be added a year in head and tail first by add_HeadTail.

Finally, use findpeaks() to get local maximum and local minimum values. Two local minimum define a growing season. If two local minimum(maximum) are too closed, then only the smaller(biger) is left.

Value

• fit: A data.table of Rough fitting result, with the columns of (t, y, witer1, ..., witerN, ziter1, ..., ziterN).

• dt: A data.table of growing season division information, with the columns of (beg, peak, end, y_beg, y_peak, y_end, len, year, season, flag).

options for season

(a) Parameters for rough fitting

• rFUN : character (default smooth_wWHIT), the name of rough curve fitting function, can be one of c("smooth_wSG", "smooth_wWHIT", "smooth_wHANTS"), which are corresponding to smooth_wSG(), smooth_wWHIT() and smooth_wHANTS().

• wFUN : character (default wTSM), the name of weights updating functions, can be one of c("wTSM", "wChen", "wBisquare", "wSELF"). See wTSM(), wChen(), wBisquare() and wSELF() for details.

• iters : integer (default 2), the number of rough fitting iterations.

• wmin : double, the minimum weight of bad points (i.e. snow, ice and cloud).

• verbose : logical (default FALSE). If TRUE, options$season will be printed on the console.

• lambda : double (default NULL), the smoothing parameter of smooth_wWHIT().

  • If lambda = NULL, V-curve theory will be employed to find the optimal lambda. See lambda_vcurve() for details.

• frame : integer (default NULL), the parameter of smooth_wSG(), moving window size.

  • If frame = NULL, frame will be reset as floor(nptperyear/5)*2 + 1 (refered by TIMESAT).

• nf : integer (default 4), the number of frequencies in smooth_wHANTS().

• maxExtendMonth: integer (default 12), previous and subsequent maxExtendMonth (in month) data were added to the current year for rough fitting.

• nextend : integer (default NULL), same as maxExtendMonth, but in points.

  • If nextend provided, maxExtendMonth will be ignored.

  • If nextend = NULL, nextend will be reset as ceiling(maxExtendMonth/12*nptperyear)

(b) Parameters for growing season division

• minpeakdistance : double (default NULL), the minimum distance of two peaks (in points). If the distance of two maximum extreme value less than minpeakdistance, only the maximum one will be kept.

  • If minpeakdistance = NULL, it will be reset as nptperyear/6.

• r_max : double (default 0.2; in (0, 1)). r_max and r_min are used to eliminate fake peaks and troughs.

  • The real peaks should satisfy:

    1. max(h_{peak, L}, h_{peak, R}) > r_{max} A

    2. min(h_{peak, L}, h_{peak, R}) > r_{min} A, where h_{peak, L}, h_{peak, R} are height difference from the peak to the left- and right-hand troughs.

  • The troughs should satisfy:

    1. max(h_{trough, L}, h_{trough, R}) > r_{max} A, where h_{trough, L}, h_{trough, R} are height difference from the trough to the left- and right-hand peaks.

• r_min : double (default 0.05; in (0, 1)), see above r_max for details. r_min < r_max.

• rtrough_max : double (default 0.6, in (0, 1)), y_{peak} <= rtrough_max * A + ylu[1].

• ypeak_min : double 0.1 (in VI unit), y_{peak} >= ypeak_min.

• .check_season : logical (default TRUE). check the growing season length according to len_min and len_max. If FALSE, len_min and len_max will lose their effect.

• len_min : integer (default 45), the minimum length (in days) of growing season

• len_max : integer (default 650), the minimum length (in days) of growing season

• adj.param : logical. If TRUE (default), if there are too many or too less peaks and troughs, phenofit will automatically adjust rough curve fitting function parameters. See MaxPeaksPerYear and MaxTroughsPerYear for details.

• MaxPeaksPerYear (optional) : integer (default 2), the max number of peaks per year. If PeaksPerYear > MaxPeaksPerYear, then lambda = lambda*2.

• MaxTroughsPerYear (optional) : integer (default 3), the max number of troughs per year. If TroughsPerYear > MaxTroughsPerYear, then lambda = lambda*2.

• calendarYear : logical (default FALSE). If TRUE, the start and end of a calendar year will be regarded as growing season division (North Hemisphere is from 01 Jan to 31 Dec; South Hemisphere is from 01 Jul to 30 Jun).

• rm.closed : logical (default TRUE). If TRUE, closed peaks (or troughs) will be further tidied. Only the maximum

• is.continuous (not used): logical (default TRUE). This parameter is for fluxnet2015 fluxsite data, where the input might be not continuous.

See Also

findpeaks().

Examples

data("CA_NS6")
d <- CA_NS6

nptperyear <- 23
INPUT <- check_input(d$t, d$y, d$w,
    QC_flag = d$QC_flag,
    nptperyear = nptperyear, south = FALSE,
    maxgap = nptperyear / 4, alpha = 0.02, wmin = 0.2
)
# plot_input(INPUT)

wFUN <- "wTSM"
# all year as a whole
options = list(rFUN = "smooth_wWHIT", wFUN = wFUN, lambda = 10)
brks <- season(INPUT, lambda = 10)
plot_season(INPUT, brks, d)

brks2 = season_input(INPUT, options)
all.equal(brks2, brks)

c(d_fit, info_peak) %<-% roughFit(INPUT)
d_season = find_season.peaks(d_fit, info_peak)

c(t, ypred) %<-% d_fit[, .(t, ziter2)]
d_season = find_season.default(ypred, t)
all.equal(brks$dt, d_season)

# opt <- .options$season
# brks$fit - d_fit # function passed test

# curve fitting by year
brks_mov <- season_mov(INPUT,
    options = list(
        rFUN = "smooth_wWHIT", wFUN = wFUN,
        lambda = 10,
        r_min = 0.05, ypeak_min = 0.05,
        verbose = TRUE
    )
)
plot_season(INPUT, brks_mov)

rfit <- brks2rfit(brks_mov)
r <- get_pheno(rfit)

Growing season division

Description

Divide growing seasons according to rough fitting (rFUN) result .

For season, rough fitting is applied for whole. For season_mov rough fitting is applied in every year, during which maxExtendMonth is extended.

Usage

season(
  INPUT,
  rFUN,
  wFUN,
  iters = 2,
  wmin = 0.1,
  lambda,
  nf = 3,
  frame = floor(INPUT$nptperyear/5) * 2 + 1,
  minpeakdistance = NULL,
  ypeak_min = 0.1,
  r_max = 0.2,
  r_min = 0.05,
  rtrough_max = 0.6,
  MaxPeaksPerYear = 2,
  MaxTroughsPerYear = 3,
  calendarYear = FALSE,
  adj.param = TRUE,
  rm.closed = TRUE,
  is.continuous = TRUE,
  .check_season = TRUE,
  verbose = FALSE,
  ...
)

stat_season(INPUT, d_fit)

Arguments

Details

Before growing season division, INPUT should be added a year in head and tail first by add_HeadTail.

Finally, use findpeaks() to get local maximum and local minimum values. Two local minimum define a growing season. If two local minimum(maximum) are too closed, then only the smaller(biger) is left.

Value

• fit: A data.table of Rough fitting result, with the columns of (t, y, witer1, ..., witerN, ziter1, ..., ziterN).

• dt: A data.table of growing season division information, with the columns of (beg, peak, end, y_beg, y_peak, y_end, len, year, season, flag).

See Also

findpeaks().

Examples

data("CA_NS6")
d <- CA_NS6

nptperyear <- 23
INPUT <- check_input(d$t, d$y, d$w,
    QC_flag = d$QC_flag,
    nptperyear = nptperyear, south = FALSE,
    maxgap = nptperyear / 4, alpha = 0.02, wmin = 0.2
)
# plot_input(INPUT)

wFUN <- "wTSM"
# all year as a whole
options = list(rFUN = "smooth_wWHIT", wFUN = wFUN, lambda = 10)
brks <- season(INPUT, lambda = 10)
plot_season(INPUT, brks, d)

brks2 = season_input(INPUT, options)
all.equal(brks2, brks)

c(d_fit, info_peak) %<-% roughFit(INPUT)
d_season = find_season.peaks(d_fit, info_peak)

c(t, ypred) %<-% d_fit[, .(t, ziter2)]
d_season = find_season.default(ypred, t)
all.equal(brks$dt, d_season)

# opt <- .options$season
# brks$fit - d_fit # function passed test

# curve fitting by year
brks_mov <- season_mov(INPUT,
    options = list(
        rFUN = "smooth_wWHIT", wFUN = wFUN,
        lambda = 10,
        r_min = 0.05, ypeak_min = 0.05,
        verbose = TRUE
    )
)
plot_season(INPUT, brks_mov)

rfit <- brks2rfit(brks_mov)
r <- get_pheno(rfit)

Growing season division (unstable version)

Description

Growing season division (unstable version)

Usage

season_input(INPUT, options = NULL, verbose = FALSE, ...)

Arguments

Moving growing season division

Description

Moving growing season division

Usage

season_mov(INPUT, options = list(), ..., years.run = NULL)

Arguments

options for season

(a) Parameters for rough fitting

• rFUN : character (default smooth_wWHIT), the name of rough curve fitting function, can be one of c("smooth_wSG", "smooth_wWHIT", "smooth_wHANTS"), which are corresponding to smooth_wSG(), smooth_wWHIT() and smooth_wHANTS().

• wFUN : character (default wTSM), the name of weights updating functions, can be one of c("wTSM", "wChen", "wBisquare", "wSELF"). See wTSM(), wChen(), wBisquare() and wSELF() for details.

• iters : integer (default 2), the number of rough fitting iterations.

• wmin : double, the minimum weight of bad points (i.e. snow, ice and cloud).

• verbose : logical (default FALSE). If TRUE, options$season will be printed on the console.

• lambda : double (default NULL), the smoothing parameter of smooth_wWHIT().

  • If lambda = NULL, V-curve theory will be employed to find the optimal lambda. See lambda_vcurve() for details.

• frame : integer (default NULL), the parameter of smooth_wSG(), moving window size.

  • If frame = NULL, frame will be reset as floor(nptperyear/5)*2 + 1 (refered by TIMESAT).

• nf : integer (default 4), the number of frequencies in smooth_wHANTS().

• maxExtendMonth: integer (default 12), previous and subsequent maxExtendMonth (in month) data were added to the current year for rough fitting.

• nextend : integer (default NULL), same as maxExtendMonth, but in points.

  • If nextend provided, maxExtendMonth will be ignored.

  • If nextend = NULL, nextend will be reset as ceiling(maxExtendMonth/12*nptperyear)

(b) Parameters for growing season division

• minpeakdistance : double (default NULL), the minimum distance of two peaks (in points). If the distance of two maximum extreme value less than minpeakdistance, only the maximum one will be kept.

  • If minpeakdistance = NULL, it will be reset as nptperyear/6.

• r_max : double (default 0.2; in (0, 1)). r_max and r_min are used to eliminate fake peaks and troughs.

  • The real peaks should satisfy:

    1. max(h_{peak, L}, h_{peak, R}) > r_{max} A

    2. min(h_{peak, L}, h_{peak, R}) > r_{min} A, where h_{peak, L}, h_{peak, R} are height difference from the peak to the left- and right-hand troughs.

  • The troughs should satisfy:

    1. max(h_{trough, L}, h_{trough, R}) > r_{max} A, where h_{trough, L}, h_{trough, R} are height difference from the trough to the left- and right-hand peaks.

• r_min : double (default 0.05; in (0, 1)), see above r_max for details. r_min < r_max.

• rtrough_max : double (default 0.6, in (0, 1)), y_{peak} <= rtrough_max * A + ylu[1].

• ypeak_min : double 0.1 (in VI unit), y_{peak} >= ypeak_min.

• .check_season : logical (default TRUE). check the growing season length according to len_min and len_max. If FALSE, len_min and len_max will lose their effect.

• len_min : integer (default 45), the minimum length (in days) of growing season

• len_max : integer (default 650), the minimum length (in days) of growing season

• adj.param : logical. If TRUE (default), if there are too many or too less peaks and troughs, phenofit will automatically adjust rough curve fitting function parameters. See MaxPeaksPerYear and MaxTroughsPerYear for details.

• MaxPeaksPerYear (optional) : integer (default 2), the max number of peaks per year. If PeaksPerYear > MaxPeaksPerYear, then lambda = lambda*2.

• MaxTroughsPerYear (optional) : integer (default 3), the max number of troughs per year. If TroughsPerYear > MaxTroughsPerYear, then lambda = lambda*2.

• calendarYear : logical (default FALSE). If TRUE, the start and end of a calendar year will be regarded as growing season division (North Hemisphere is from 01 Jan to 31 Dec; South Hemisphere is from 01 Jul to 30 Jun).

• rm.closed : logical (default TRUE). If TRUE, closed peaks (or troughs) will be further tidied. Only the maximum

• is.continuous (not used): logical (default TRUE). This parameter is for fluxnet2015 fluxsite data, where the input might be not continuous.

References

1. Kong, D., Zhang, Y., Wang, D., Chen, J., & Gu, X. (2020). Photoperiod Explains the Asynchronization Between Vegetation Carbon Phenology and Vegetation Greenness Phenology. Journal of Geophysical Research: Biogeosciences, 125(8), e2020JG005636. https://doi.org/10.1029/2020JG005636

2. Kong, D., Zhang, Y., Gu, X., & Wang, D. (2019). A robust method for reconstructing global MODIS EVI time series on the Google Earth Engine. ISPRS Journal of Photogrammetry and Remote Sensing, 155, 13-24.

See Also

season()

Examples

data("CA_NS6")
d <- CA_NS6

nptperyear <- 23
INPUT <- check_input(d$t, d$y, d$w,
    QC_flag = d$QC_flag,
    nptperyear = nptperyear, south = FALSE,
    maxgap = nptperyear / 4, alpha = 0.02, wmin = 0.2
)

# curve fitting by year
brks_mov <- season_mov(INPUT,
    options = list(
        rFUN = "smooth_wWHIT", wFUN = "wTSM",
        lambda = 10,
        r_min = 0.05, ypeak_min = 0.05,
        verbose = TRUE
    )
)
plot_season(INPUT, brks_mov)

rfit <- brks2rfit(brks_mov)
# Phenological Metrics from rough fitting
r <- get_pheno(rfit)

set and get phenofit option

Description

set and get phenofit option

Usage

set_options(..., options = NULL)

get_options(names = NULL)

Arguments

options for season

(a) Parameters for rough fitting

• rFUN : character (default smooth_wWHIT), the name of rough curve fitting function, can be one of c("smooth_wSG", "smooth_wWHIT", "smooth_wHANTS"), which are corresponding to smooth_wSG(), smooth_wWHIT() and smooth_wHANTS().

• wFUN : character (default wTSM), the name of weights updating functions, can be one of c("wTSM", "wChen", "wBisquare", "wSELF"). See wTSM(), wChen(), wBisquare() and wSELF() for details.

• iters : integer (default 2), the number of rough fitting iterations.

• wmin : double, the minimum weight of bad points (i.e. snow, ice and cloud).

• verbose : logical (default FALSE). If TRUE, options$season will be printed on the console.

• lambda : double (default NULL), the smoothing parameter of smooth_wWHIT().

  • If lambda = NULL, V-curve theory will be employed to find the optimal lambda. See lambda_vcurve() for details.

• frame : integer (default NULL), the parameter of smooth_wSG(), moving window size.

  • If frame = NULL, frame will be reset as floor(nptperyear/5)*2 + 1 (refered by TIMESAT).

• nf : integer (default 4), the number of frequencies in smooth_wHANTS().

• maxExtendMonth: integer (default 12), previous and subsequent maxExtendMonth (in month) data were added to the current year for rough fitting.

• nextend : integer (default NULL), same as maxExtendMonth, but in points.

  • If nextend provided, maxExtendMonth will be ignored.

  • If nextend = NULL, nextend will be reset as ceiling(maxExtendMonth/12*nptperyear)

(b) Parameters for growing season division

• minpeakdistance : double (default NULL), the minimum distance of two peaks (in points). If the distance of two maximum extreme value less than minpeakdistance, only the maximum one will be kept.

  • If minpeakdistance = NULL, it will be reset as nptperyear/6.

• r_max : double (default 0.2; in (0, 1)). r_max and r_min are used to eliminate fake peaks and troughs.

  • The real peaks should satisfy:

    1. max(h_{peak, L}, h_{peak, R}) > r_{max} A

    2. min(h_{peak, L}, h_{peak, R}) > r_{min} A, where h_{peak, L}, h_{peak, R} are height difference from the peak to the left- and right-hand troughs.

  • The troughs should satisfy:

    1. max(h_{trough, L}, h_{trough, R}) > r_{max} A, where h_{trough, L}, h_{trough, R} are height difference from the trough to the left- and right-hand peaks.

• r_min : double (default 0.05; in (0, 1)), see above r_max for details. r_min < r_max.

• rtrough_max : double (default 0.6, in (0, 1)), y_{peak} <= rtrough_max * A + ylu[1].

• ypeak_min : double 0.1 (in VI unit), y_{peak} >= ypeak_min.

• .check_season : logical (default TRUE). check the growing season length according to len_min and len_max. If FALSE, len_min and len_max will lose their effect.

• len_min : integer (default 45), the minimum length (in days) of growing season

• len_max : integer (default 650), the minimum length (in days) of growing season

• adj.param : logical. If TRUE (default), if there are too many or too less peaks and troughs, phenofit will automatically adjust rough curve fitting function parameters. See MaxPeaksPerYear and MaxTroughsPerYear for details.

• MaxPeaksPerYear (optional) : integer (default 2), the max number of peaks per year. If PeaksPerYear > MaxPeaksPerYear, then lambda = lambda*2.

• MaxTroughsPerYear (optional) : integer (default 3), the max number of troughs per year. If TroughsPerYear > MaxTroughsPerYear, then lambda = lambda*2.

• calendarYear : logical (default FALSE). If TRUE, the start and end of a calendar year will be regarded as growing season division (North Hemisphere is from 01 Jan to 31 Dec; South Hemisphere is from 01 Jul to 30 Jun).

• rm.closed : logical (default TRUE). If TRUE, closed peaks (or troughs) will be further tidied. Only the maximum

• is.continuous (not used): logical (default TRUE). This parameter is for fluxnet2015 fluxsite data, where the input might be not continuous.

options for fitting

• methods (default ⁠c('AG', 'Beck', 'Elmore', 'Zhang')``): Fine curve fitting methods, can be one or more of ⁠c('AG', 'Beck', 'Elmore', 'Zhang', 'Gu', 'Klos')‘. Note that ’Gu' and 'Klos' are very slow.

• iters (default 2): max iterations of fine fitting.

• wFUN (default wTSM): Character or function, weights updating function of fine fitting function.

• wmin (default 0.1): min weights in the weights updating procedure.

• use.rough (default FALSE): Whether to use rough fitting smoothed time-series as input? If false, smoothed VI by rough fitting will be used for Phenological metrics extraction; If true, original input y will be used (rough fitting is used to divide growing seasons and update weights.

• use.y0 (default TRUE): boolean. whether to use original y0 as the input of plot_input, note that not for curve fitting. y0 is the original value before the process of check_input.

• nextend (default 2): Extend curve fitting window, until nextend good or marginal points are found in the previous and subsequent growing season.

• maxExtendMonth (default 1): Search good or marginal good values in previous and subsequent maxExtendMonth period.

• minExtendMonth (default 0.5): Extend period defined by nextend and maxExtendMonth, should be no shorter than minExtendMonth. When all points of the input time-series are good value, then the extending period will be too short. In that situation, we can't make sure the connection between different growing seasons is smoothing.

• minPercValid: (default 0, not use). If the percentage of good- and marginal- quality points is less than minPercValid, curve fiting result is set to NA.

• minT: (not use). If Tn not provided in INPUT, minT will not be used. minT use night temperature Tn to define backgroud value (days with Tn < minT treated as ungrowing season).

Examples

set_options(verbose = FALSE) 
get_options("season") %>% str()

Weighted HANTS SMOOTH

Description

Weighted HANTS smoother

Usage

smooth_wHANTS(
  y,
  t,
  w,
  nf = 3,
  ylu,
  periodlen = 365,
  nptperyear,
  wFUN = wTSM,
  iters = 2,
  wmin = 0.1,
  ...
)

Arguments

Value

• ws: weights of every iteration

• zs: curve fittings of every iteration

Author(s)

Wout Verhoef, NLR, Remote Sensing Dept. June 1998 Mohammad Abouali (2011), Converted to MATLAB Dongdong Kong (2018), introduced to R and modified into weighted model.

Examples

library(phenofit)
data("MOD13A1")
dt <- tidy_MOD13(MOD13A1$dt)
d <- dt[site == "AT-Neu", ]

l <- check_input(d$t, d$y, d$w, nptperyear=23)
r_wHANTS <- smooth_wHANTS(l$y, l$t, l$w, ylu = l$ylu, nptperyear = 23, iters = 2)

Weighted Savitzky-Golay

Description

Weighted Savitzky-Golay

Usage

smooth_wSG(
  y,
  w,
  nptperyear,
  ylu,
  wFUN = wTSM,
  iters = 2,
  frame = floor(nptperyear/5) * 2 + 1,
  d = 2,
  ...
)

Arguments

Value

• ws: weights of every iteration

• zs: curve fittings of every iteration

References

1. Chen, J., J\"onsson, P., Tamura, M., Gu, Z., Matsushita, B., Eklundh, L., 2004. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter. Remote Sens. Environ. 91, 332-344. https://doi.org/10.1016/j.rse.2004.03.014.
   

2. https://en.wikipedia.org/wiki/Savitzky%E2%80%93Golay_filter

Examples

library(phenofit)
data("MOD13A1")
dt <- tidy_MOD13(MOD13A1$dt)
d <- dt[site == "AT-Neu", ]

l <- check_input(d$t, d$y, d$w, nptperyear=23)
r_wSG <- smooth_wSG(l$y, l$w, l$ylu, nptperyear = 23, iters = 2)

Weigthed Whittaker Smoother

Description

Weigthed Whittaker Smoother

Usage

smooth_wWHIT(
  y,
  w,
  ylu,
  nptperyear,
  wFUN = wTSM,
  iters = 1,
  lambda = 15,
  second = FALSE,
  ...
)

Arguments

Value

• ws: weights of every iteration

• zs: curve fittings of every iteration

Note

Whittaker smoother of the second order difference is used!

References

1. Eilers, P.H.C., 2003. A perfect smoother. Anal. Chem. doi:10.1021/ac034173t

2. Frasso, G., Eilers, P.H.C., 2015. L- and V-curves for optimal smoothing. Stat. Modelling 15, 91-111. doi:10.1177/1471082X14549288.

See Also

lambda_vcurve()

Examples

data("MOD13A1")
dt <- tidy_MOD13(MOD13A1$dt)
d <- dt[site == "AT-Neu", ]

l <- check_input(d$t, d$y, d$w, nptperyear=23)
r_wWHIT <- smooth_wWHIT(l$y, l$w, l$ylu, nptperyear = 23, iters = 2)

## Optimize `lambda` by V-curve theory
# (a) optimize manually
lambda_vcurve(l$y, l$w, plot = TRUE) 

# (b) optimize automatically by setting `lambda = NULL` in smooth_wWHIT
r_wWHIT2 <- smooth_wWHIT(l$y, l$w, l$ylu, nptperyear = 23, iters = 2, lambda = NULL) # 

tidy_MOD13

Description

Tidy MODIS 'MOD13' VI products' (e.g. MOD13A1, MOD13A2, ...) raw data exported from Google Earth Engine. Tidy contents include:

1. add exact compositing date, see getRealDate().
   

2. Init weigths according SummaryQA, see qc_summary().
   

Usage

tidy_MOD13(infile, outfile, wmin = 0.2)

Arguments

Value

A tidied data.table, with columns of 'site', 'y', 't', 'w', 'date' and 'SummaryQA'.

• site: site name

• y: real value of EVI, ⁠[-1, 1]⁠

• date: image date

• t: exact compositing date constructed from DayOfYear

• w: weights

• SummaryQA: A factor, QA types, one of "good", "margin", "snow/ice" or "cloud".

Examples

library(phenofit)
data("MOD13A1")
dt <- tidy_MOD13(MOD13A1$dt)

tidy_pheno

Description

Tidy for every method with multiple years phenology data

Usage

tidy_pheno(pheno)

date2doy(p_date)

doy2date(p_doy)

Arguments

Examples

library(phenofit)
# 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)

methods <- c("AG", "Beck", "Elmore", "Gu", "Zhang") # "Klos" too slow
fFITs <- curvefit(y, t, tout, methods)

# multiple years
fits <- list(`2001` = fFITs, `2002` = fFITs)
pheno <- get_pheno(fits, "AG", IsPlot = FALSE)

V-curve theory to optimize Whittaker parameter lambda.

Description

V-curve is used to optimize Whittaker parameter lambda. This function is not for users!!!

Update 20180605 add weights updating to whittaker lambda selecting

Usage

v_curve(INPUT, lg_lambdas = seq(0, 5, by = 0.005), plot = FALSE, ...)

Arguments

See Also

lambda_vcurve

Examples

data("CA_NS6"); d = CA_NS6
nptperyear = 23
INPUT <- check_input(d$t, d$y, d$w, nptperyear = nptperyear,
    maxgap = nptperyear/4, alpha = 0.02, wmin = 0.2)

r <- v_curve(INPUT, lg_lambdas = seq(0, 3, 0.1), plot = TRUE)

Weight updating functions

Description

• wSELF weigth are not changed and return the original.

• wTSM weight updating method in TIMESAT.

• wBisquare Bisquare weight update method. wBisquare has been modified to emphasis on upper envelope.

• wBisquare0 Traditional Bisquare weight update method.

• wChen Chen et al., (2004) weight updating method.

• wBeck Beck et al., (2006) weigth updating method. wBeck need sos and eos input. The function parameter is different from others. It is still not finished.

Usage

wSELF(y, yfit, w, ...)

wTSM(y, yfit, w, iter = 2, nptperyear, wfact = 0.5, ...)

wBisquare0(y, yfit, w, ..., wmin = 0.2)

wBisquare(y, yfit, w, ..., wmin = 0.2, .toUpper = TRUE)

wChen(y, yfit, w, ..., wmin = 0.2)

wKong(y, yfit, w, ..., wmin = 0.2)

Arguments

Value

wnew Numeric Vector, adjusted weights.

Author(s)

wTSM is implemented by Per J\"onsson, Malm\"o University, Sweden per.jonsson@ts.mah.se and Lars Eklundh, Lund University, Sweden lars.eklundh@nateko.lu.se. And Translated into Rcpp by Dongdong Kong, 01 May 2018.

References

1. Per J\"onsson, P., Eklundh, L., 2004. TIMESAT - A program for analyzing time-series of satellite sensor data. Comput. Geosci. 30, 833-845. https://doi.org/10.1016/j.cageo.2004.05.006.
   

2. https://au.mathworks.com/help/curvefit/smoothing-data.html#bq_6ys3-3
   

3. Garcia, D., 2010. Robust smoothing of gridded data in one and higher dimensions with missing values. Computational statistics & data analysis, 54(4), pp.1167-1178.
   

4. Chen, J., J\"onsson, P., Tamura, M., Gu, Z., Matsushita, B., Eklundh, L., 2004. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter. Remote Sens. Environ. 91, 332-344. https://doi.org/10.1016/j.rse.2004.03.014.
   

5. Beck, P.S.A., Atzberger, C., Hogda, K.A., Johansen, B., Skidmore, A.K., 2006. Improved monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI. Remote Sens. Environ. https://doi.org/10.1016/j.rse.2005.10.021
   

6. https://github.com/kongdd/phenopix/blob/master/R/FitDoubleLogBeck.R

Weighted Whittaker smoothing with a second order finite difference penalty

Description

This function smoothes signals with a finite difference penalty of order 2. This function is modified from ptw package.

Usage

whit2(y, lambda, w = rep(1, ny))

Arguments

Value

A numeric vector, smoothed signal.

Author(s)

Paul Eilers, Jan Gerretzen

References

1. Eilers, P.H.C. (2004) "Parametric Time Warping", Analytical Chemistry, 76 (2), 404 – 411.
   

2. Eilers, P.H.C. (2003) "A perfect smoother", Analytical Chemistry, 75, 3631 – 3636.

Examples

library(phenofit)
data("MOD13A1")
dt <- tidy_MOD13(MOD13A1$dt)
y <- dt[site == "AT-Neu", ][1:120, y]

plot(y, type = "b")
lines(whit2(y, lambda = 2), col = 2)
lines(whit2(y, lambda = 10), col = 3)
lines(whit2(y, lambda = 100), col = 4)
legend("bottomleft", paste("lambda = ", c(2, 10, 15)), col = 2:4, lty = rep(1, 3))