---
title: "Extract vegetation phenology from MOD13A1 EVI"
output: 
  rmarkdown::html_vignette:
    number_sections: true
vignette: >
  %\VignetteIndexEntry{phenofit-MODIS}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

# Example

Here, we illustrate how to use `phenofit` to extract vegetation phenology from 
MOD13A1 in the sampled points. Regional analysis also can be conducted in the 
similar way.

<!-- ## 1.1 Download MOD13A1 data
Upload point shapefile into GEE, clip MOD13A1 and download vegetation index
data. [Here](https://code.earthengine.google.com/ee3ec39cf3061374dab435c358d008a3) is the corresponding GEE script. 
 -->

## 1.1 Initial weights for input data

Load packages.
```{r load pkg, message=FALSE, echo=FALSE}
library(data.table)
library(lubridate)
library(purrr)
library(ggplot2)
library(phenofit)
```
Set global parameters for `phenofit`
```{r phenofit_parameters}
# lambda   <- 5    # non-parameter Whittaker, only suit for 16-day. Other time-scale
# should assign a lambda.
ymax_min   <- 0.1  # the maximum ymax shoud be greater than `ymax_min` 
rymin_less <- 0.8  # trough < ymin + A*rymin_less
nptperyear <- 23   # How many points for a single year
wFUN       <- wBisquare #wTSM #wBisquare # Weights updating function, could be one of `wTSM`, 'wBisquare', `wChen` and `wSELF`.
```

- Add date according to composite day of the year (DayOfYear), other than image date.
- Add weights according to `SummaryQA`. 

For MOD13A1, Weights can by initialed by `SummaryQA` band (also suit for 
MOD13A2 and MOD13Q1). There is already a `QC` function for `SummaryQA`, i.e. `qc_summary`.

SummaryQA      | Pixel reliability summary QA | weight
---------------| ---------------------------- | ------
-1 Fill/No data| Not processed                | `wmin`
0 Good data    | Use with confidence          | 1
1 Marginal data| Useful but look at detailed QA for more information | 0.5
2 Snow/ice     | Pixel covered with snow/ice  | `wmin`
3 Cloudy       | Pixel is cloudy              | `wmin`

```{r tidy MOD13A1}
data('MOD13A1')
df <- MOD13A1$dt 
st <- MOD13A1$st

df[, `:=`(date = ymd(date), year = year(date), doy = as.integer(yday(date)))]
df[is.na(DayOfYear), DayOfYear := doy] # If DayOfYear is missing
    
# In case of last scene of a year, doy of last scene could in the next year
df[abs(DayOfYear - doy) >= 300, t := as.Date(sprintf("%d-%03d", year+1, DayOfYear), "%Y-%j")] # last scene
df[abs(DayOfYear - doy) <  300, t := as.Date(sprintf("%d-%03d", year  , DayOfYear), "%Y-%j")]

df <- df[!duplicated(df[, .(site, t)]), ]
# # MCD12Q1.006 land cover 1-17, IGBP scheme
# IGBPnames_006 <- c("ENF", "EBF", "DNF", "DBF", "MF" , "CSH", 
#               "OSH", "WSA", "SAV", "GRA", "WET", "CRO", 
#               "URB", "CNV", "SNOW", "BSV", "water", "UNC")
# Initial weights
df[, c("QC_flag", "w") := qc_summary(SummaryQA)]
df <- df[, .(site, y = EVI/1e4, t, date, w, QC_flag)]
```

- `add_HeadTail` is used to deal with such situation, e.g. MOD13A2 begins from 2000-02-08. 
We need to construct a series with complete year, which begins from 01-01 for NH, and 07-01 for SH. 
For example, the input data period is 20000218 ~ 20171219. After adding one year in head and 
tail, it becomes 19990101 ~ 20181219. 

## 2.1 load site data
```{r load_data}
sites        <- unique(df$site)
sitename     <- sites[3]
d            <- df[site == sitename] # get the first site data
sp           <- st[site == sitename]

south      <- sp$lat < 0
print      <- FALSE # whether print progress
IsPlot     <- TRUE  # for brks

prefix_fig <- "phenofit"
titlestr   <- with(sp, sprintf('[%03d,%s] %s, lat = %5.2f, lon = %6.2f',
                                     ID, site, IGBPname, lat, lon))
file_pdf   <- sprintf('Figure/%s_[%03d]_%s.pdf', prefix_fig, sp$ID[1], sp$site[1])
```

If need night temperature (Tn) to constrain ungrowing season backgroud value, NA 
values in Tn should be filled.
```{r interp Tn, eval=F}
d$Tn %<>% zoo::na.approx(maxgap = 4)
plot(d$Tn, type = "l"); abline(a = 5, b = 0, col = "red")
```

## 2.2 Check input data
```{r check_input}    
dnew  <- add_HeadTail(d, south, nptperyear = 23) # add additional one year in head and tail
INPUT <- check_input(dnew$t, dnew$y, dnew$w, dnew$QC_flag,
                     nptperyear, south, 
                     maxgap = nptperyear/4, alpha = 0.02, wmin = 0.2)
```

## 2.3 Divide growing seasons

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 divide growing season}
par(mar = c(3, 2, 2, 1), mgp = c(3, 0.6, 0))
lambda <- init_lambda(INPUT$y)
# The detailed information of those parameters can be seen in `season`.
# brks   <- season(INPUT, nptperyear,
#                FUN = smooth_wWHIT, wFUN = wFUN, iters = 2,
#                lambda = lambda,
#                IsPlot = IsPlot, plotdat = d,
#                south = d$lat[1] < 0,
#                rymin_less = 0.6, ymax_min = ymax_min,
#                max_MaxPeaksperyear =2.5, max_MinPeaksperyear = 3.5) #, ...
# get growing season breaks in a 3-year moving window
brks2 <- season_mov(INPUT, 
    list(rFUN = "smooth_wWHIT", wFUN = wFUN, maxExtendMonth = 6, r_min = 0.1))
plot_season(INPUT, brks2)
```

## 2.4 Curve fitting
```{r curve fitting, fig.height=7, fig.align="center"}
fit  <- curvefits(INPUT, brks2,
    options = list(
        methods = c("AG", "Zhang", "Beck", "Elmore"), #,"klos",, 'Gu'
        wFUN = wFUN,
        nextend = 2, maxExtendMonth = 3, minExtendMonth = 1, minPercValid = 0.2
    ))

## check the curve fitting parameters
l_param <- get_param(fit)
print(str(l_param, 1))
print(l_param$AG)

d_fit <- get_fitting(fit)
## Get GOF information
d_gof <- get_GOF(fit)
# fit$stat <- stat
print(head(d_gof))

# print(fit$fits$AG$`2002_1`$ws)
print(fit$`2002_1`$fFIT$AG$ws)
## visualization
g <- plot_curvefits(d_fit, brks2, NULL, ylab = "NDVI", "Time",
                   theme = coord_cartesian(xlim = c(ymd("2000-04-01"), ymd("2017-07-31"))))
grid::grid.newpage(); grid::grid.draw(g)# plot to check the curve fitting
# write_fig(g, "Figure1_phenofit_curve_fitting.pdf", 10, 6)
```

## 2.5 Extract phenology
```{r Extract phenology, fig.height=5, fig.width=8, fig.align="center"}
# pheno: list(p_date, p_doy)
l_pheno <- get_pheno(fit, IsPlot = F) #%>% map(~melt_list(., "meth"))

# ratio = 1.15
# file <- "Figure5_Phenology_Extraction_temp.pdf"
# cairo_pdf(file, 8*ratio, 6*ratio)
# temp <- get_pheno(fit$fits$ELMORE[2:6], IsPlot = T)
# dev.off()
# file.show(file)

## check the extracted phenology
pheno <- get_pheno(fit[1:6], "Elmore", IsPlot = T)
# print(str(pheno, 1))
head(l_pheno$doy$AG)
```

