Detailed procedures of phenofit"

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)


Try the phenofit package in your browser

Any scripts or data that you put into this service are public.

phenofit documentation built on Oct. 15, 2021, 5:09 p.m.