Detailed procedures of phenofit"
In phenofit: Extract Remote Sensing Vegetation Phenology

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=7, fig.height=5
)

knitr::opts_chunk$set()

Example

Load packages

library(phenofit)
library(data.table)
library(dplyr)
library(ggplot2)

Load data

d = MOD13A1$dt %>% subset(site == "CA-NS6" & date >= "2010-01-01" & date <= "2016-12-31") %>%
    .[, .(date, y = EVI/1e4, DayOfYear, QC = SummaryQA)]
d %<>% mutate(t = getRealDate(date, DayOfYear)) %>%
    cbind(d[, as.list(qc_summary(QC, wmin = 0.2, wmid = 0.5, wmax = 0.8))]) %>%
    .[, .(date, t, y, QC_flag, w)]
print(d)

date: image date t : composite date

QC_flag and date are optional.

phenofit parameters

lambda         <- 8
nptperyear     <- 23
minExtendMonth <- 0.5
maxExtendMonth <- 1
minPercValid   <- 0
wFUN           <- wTSM # wBisquare
wmin           <- 0.2
methods_fine <- c("AG", "Zhang", "Beck", "Elmore", "Gu")

growing season division and fine-curve fitting

Simply treating calendar year as a complete growing season will induce a considerable error for phenology extraction. A simple growing season dividing method was proposed in phenofit.

The growing season dividing method rely on heavily in Whittaker smoother.

Procedures of initial weight, growing season dividing, curve fitting, and phenology extraction are conducted separately.

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

brks <- season_mov(INPUT,
    list(FUN = "smooth_wWHIT", wFUN = wFUN,
        maxExtendMonth = 3,
        wmin = wmin, r_min = 0.1
    ))
# plot_season(INPUT, brks)

## 2.4 Curve fitting
fit <- curvefits(INPUT, brks,
    list(
        methods = methods_fine, # ,"klos",, 'Gu'
        wFUN = wFUN,
        iters = 2,
        wmin = wmin,
        # constrain = FALSE,
        nextend = 2,
        maxExtendMonth = maxExtendMonth, minExtendMonth = minExtendMonth,
        minPercValid = minPercValid
    ))

## check the curve fitting parameters
l_param <- get_param(fit)
print(l_param$Beck)

dfit <- get_fitting(fit)
print(dfit)

## 2.5 Extract phenology
TRS <- c(0.1, 0.2, 0.5)
l_pheno <- get_pheno(fit, TRS = TRS, IsPlot = FALSE) # %>% map(~melt_list(., "meth"))
print(l_pheno$doy$Beck)

pheno <- l_pheno$doy %>% melt_list("meth")

Visualization

# growing season dividing
plot_season(INPUT, brks, ylab = "EVI")
# Ipaper::write_fig({  }, "Figure4_seasons.pdf", 9, 4)

# fine curvefitting
g <- plot_curvefits(dfit, brks, title = NULL, cex = 1.5, ylab = "EVI", angle = 0)
grid::grid.newpage()
grid::grid.draw(g)
# Ipaper::write_fig(g, "Figure5_curvefitting.pdf", 8, 6, show = TRUE)

# extract phenology metrics, only the first 3 year showed at here
# write_fig({
l_pheno <- get_pheno(fit[1:7], method = "AG", TRS = TRS, IsPlot = TRUE, show.title = FALSE)
# }, "Figure6_phenology_metrics.pdf", 8, 6, show = TRUE)

Comparison with `TIMESAT` and `phenopix`

# library(ggplot2)
# library(ggnewscale)

# # on the top of `Figure7_predata...`
# d_comp = fread("data-raw/dat_Figure7_comparison_with_others-CA-NS6.csv")
# d_comp = merge(d[, .(date, t)], d_comp[, .(date, TIMESAT, phenopix)]) %>%
#     merge(dfit[meth == "Beck", .(t, phenofit = ziter2)], by = "t") %>%
#     melt(c("date", "t"), variable.name = "meth")

# labels = c("good", "marginal", "snow", "cloud")
# theme_set(theme_grey(base_size = 16))
# cols_line = c(phenofit = "red", TIMESAT = "blue", phenopix = "black")
# p <- ggplot(dfit, aes(t, y)) +
#     geom_point(aes(color = QC_flag, fill = QC_flag, shape = QC_flag), size = 3) +
#     scale_shape_manual(values = qc_shapes[labels], guide = guide_legend(order = 1)) +
#     scale_color_manual(values = qc_colors[labels], guide = guide_legend(order = 1)) +
#     scale_fill_manual(values = qc_colors[labels], guide = guide_legend(order = 1)) +
#     new_scale_color() +
#     geom_line(data = d_comp, aes(t, value, color = meth)) +
#     # geom_line(data = d_comp[meth == "phenofit"], aes(t, value),
#     #           size = 1, show.legend = FALSE, color = "red") +
#     scale_color_manual(values = cols_line, guide = guide_legend(order = 2)) +
#     labs(x = "Time", y = "EVI") +
#     theme(
#         axis.title.x = element_text(margin = margin(t = 0, unit='cm')),
#         # plot.margin = margin(t = 0, unit='cm'),
#         legend.text = element_text(size = 13),
#         legend.position = "bottom",
#             legend.title  = element_blank(),
#             legend.margin = margin(t = -0.3, unit='cm'))
# p
# # write_fig(p, "Figure7_comparison_with_others.pdf", 10, 4, show = TRUE)