---
title: "Detailed procedures of phenofit"
output: 
  rmarkdown::html_vignette:
    number_sections: true
vignette: >
  %\VignetteIndexEntry{phenofit-procedures}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=7, fig.height=5
)
```

```{r,include=FALSE}
knitr::opts_chunk$set() 
```

# Example

## Load packages
```{r,message=F}
library(phenofit)
library(data.table)
library(dplyr)
library(ggplot2)
```

## Load data
```{r}
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
```{r}
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.

```{r}
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

```{r}
# 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`

```{r}
# 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)
```