# phenofit
**Repository Path**: adv-r/phenofit
## Basic Information
- **Project Name**: phenofit
- **Description**: No description available
- **Primary Language**: Unknown
- **License**: GPL-2.0
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No
## Statistics
- **Stars**: 0
- **Forks**: 0
- **Created**: 2020-03-14
- **Last Updated**: 2020-12-19
## Categories & Tags
**Categories**: Uncategorized
**Tags**: None
## README
# phenofit
[](https://travis-ci.org/kongdd/phenofit)
[](https://ci.appveyor.com/project/kongdd/phenofit)
[](https://codecov.io/gh/kongdd/phenofit)
[](http://www.gnu.org/licenses/gpl-2.0.html)
[](https://cran.r-project.org/package=phenofit)
[](https://www.rpackages.io/package/phenofit)
[](https://www.rpackages.io/package/phenofit)
A state-of-the-art **remote sensing vegetation phenology** extraction
package: `phenofit`
- `phenofit` combine merits of TIMESAT and phenopix
- A simple and stable growing season dividing methods was proposed
- Provide a practical snow elimination method, based on Whittaker
- 7 curve fitting methods and 4 phenology extraction methods
- We add parameters boundary for every curve fitting methods according
to their ecological meaning.
- `optimx` is used to select best optimization method for different
curve fitting methods.
***Task lists***
- [ ] Test the performance of `phenofit` in multiple growing season
regions (e.g. the North China Plain);
- [ ] Uncertainty analysis of curve fitting and phenological metrics;
- [x] shiny app has been moved to
[phenofit.shiny](https://github.com/kongdd/phenofit.shiny);
- [x] Complete script automatic generating module in shinyapp;
- [x] `Rcpp` improve double logistics optimization efficiency by 60%;
- [x] Support spatial analysis;
- [x] Support annual season in curve fitting;
- [x] flexible fine fitting input ( original time-series or smoothed
time-series by rough fitting).
- [x] Asymmetric of Threshold method
*Figure 1. The flowchart of phenology extraction in `phenofit`.*
# Installation
You can install phenofit from github with:
``` r
# install.packages("devtools")
devtools::install_github("kongdd/phenofit")
```
# 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 Initial weights for input data
Load packages.
``` r
suppressMessages({
library(data.table)
library(magrittr)
library(lubridate)
library(purrr)
library(plyr)
library(ggplot2)
library(phenofit)
})
```
Set global parameters for `phenofit`
``` r
# 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
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
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
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
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
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,
FUN = smooth_wWHIT, wFUN = wFUN,
maxExtendMonth = 6, r_min = 0.1,
IsPlot = IsPlot, IsPlot.OnlyBad = FALSE, print = print)
```
## 2.4 Curve fitting
``` r
fit <- curvefits(INPUT, brks2,
methods = c("AG", "Zhang", "Beck", "Elmore"), #,"klos",, 'Gu'
wFUN = wFUN,
nextend = 2, maxExtendMonth = 3, minExtendMonth = 1, minPercValid = 0.2,
print = print, verbose = FALSE)
## check the curve fitting parameters
l_param <- get_param(fit)
print(str(l_param, 1))
# List of 4
# $ AG :Classes 'tbl_df', 'tbl' and 'data.frame': 18 obs. of 8 variables:
# $ Beck :Classes 'tbl_df', 'tbl' and 'data.frame': 18 obs. of 7 variables:
# $ Elmore:Classes 'tbl_df', 'tbl' and 'data.frame': 18 obs. of 8 variables:
# $ Zhang :Classes 'tbl_df', 'tbl' and 'data.frame': 18 obs. of 8 variables:
# NULL
print(l_param$AG)
# # A tibble: 18 x 8
# flag t0 mn mx rsp a3 rau a5
#
# 1 2000_1 201. 0.168 0.407 0.0323 3.53 0.0151 5.83
# 2 2001_1 561. 0.173 0.407 0.0225 4.06 0.0162 6
# 3 2002_1 931. 0.188 0.511 0.0383 2 0.0169 4.21
# 4 2003_1 1275. 0.167 0.434 0.0268 2 0.0122 3.37
# 5 2004_1 1660. 0.175 0.451 0.0363 2 0.0179 4.46
# 6 2005_1 2052. 0.180 0.466 0.0141 6 0.0314 2
# 7 2006_1 2379. 0.174 0.436 0.0226 2.51 0.0131 3.26
# 8 2007_1 2752. 0.165 0.483 0.0208 2 0.0150 2.88
# 9 2008_1 3134. 0.177 0.492 0.0180 3.50 0.0199 6
# 10 2009_1 3525. 0.172 0.480 0.0133 5.40 0.0313 2
# 11 2010_1 3841. 0.181 0.493 0.0238 2 0.0148 2
# 12 2011_1 4206. 0.190 0.464 0.0328 2 0.0135 6
# 13 2012_1 4558. 0.167 0.510 0.0490 2 0.0109 3.64
# 14 2013_1 4966. 0.168 0.484 0.0140 6 0.0190 2
# 15 2014_1 5303. 0.203 0.489 0.0344 2 0.0133 6
# 16 2015_1 5690. 0.215 0.494 0.0182 6 0.0275 4.81
# 17 2016_1 6023. 0.194 0.484 0.0429 2 0.0126 4.99
# 18 2017_1 6406. 0.171 0.449 0.0201 6 0.0129 3.52
d_fit <- get_fitting(fit)
## Get GOF information
d_gof <- get_GOF(fit)
# fit$stat <- stat
print(head(d_gof))
# flag meth RMSE NSE R pvalue n
# 1: 2000_1 AG 0.10058989 0.4952465 0.9016506 1.809700e-09 24
# 2: 2000_1 Beck 0.10059953 0.4951497 0.9036266 1.461349e-09 24
# 3: 2000_1 Elmore 0.10060466 0.4950983 0.9021757 1.710497e-09 24
# 4: 2000_1 Zhang 0.10059866 0.4951585 0.9036627 1.455586e-09 24
# 5: 2001_1 AG 0.08808864 0.5879431 0.9037624 3.439665e-09 23
# 6: 2001_1 Beck 0.08817911 0.5870963 0.9047721 3.093164e-09 23
# print(fit$fits$AG$`2002_1`$ws)
print(fit$`2002_1`$fFIT$AG$ws)
# $iter1
# [1] 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.5 1.0 1.0 1.0 1.0 0.5 1.0 1.0
# [18] 1.0 1.0 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 1.0 1.0
#
# $iter2
# [1] 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000
# [8] 0.2000000 0.2000000 0.5000000 1.0000000 1.0000000 0.7893075 1.0000000
# [15] 0.4984446 0.8873205 1.0000000 1.0000000 1.0000000 0.2000000 0.2000000
# [22] 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000 0.2000000
# [29] 0.2000000 0.2000000 0.2000000 1.0000000 1.0000000
## visualization
g <- plot_phenofit(d_fit, brks2, NULL, title.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
```
``` r
# write_fig(g, "Figure1_phenofit_curve_fitting.pdf", 10, 6)
```
## 2.5 Extract phenology
``` r
# 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)
```
``` r
# print(str(pheno, 1))
head(l_pheno$doy$AG)
# flag origin TRS2.sos TRS2.eos TRS5.sos TRS5.eos TRS6.sos TRS6.eos
# 1: 2000_1 2000-01-01 167 273 175 264 177 261
# 2: 2001_1 2001-01-01 145 263 155 254 158 251
# 3: 2002_1 2002-01-01 168 268 179 256 182 252
# 4: 2003_1 2003-01-01 133 273 149 253 153 246
# 5: 2004_1 2004-01-01 166 262 178 251 181 248
# 6: 2005_1 2005-01-01 150 266 159 252 162 248
# DER.sos DER.pop DER.eos UD SD DD RD Greenup Maturity Senescence
# 1: 174 203 266 164 186 249 280 157 193 243
# 2: 154 196 256 141 170 240 268 133 177 234
# 3: 183 202 257 164 195 238 274 157 201 223
# 4: 153 180 254 128 170 225 283 118 179 196
# 5: 181 201 253 162 193 236 268 154 200 226
# 6: 157 225 248 143 175 233 272 136 181 279
# Dormancy
# 1: 287
# 2: 275
# 3: 282
# 4: 298
# 5: 276
# 6: NA
```
# **References**
> \[1\] Dongdong Kong, R package: A state-of-the-art Vegetation
> Phenology extraction package, `phenofit` version 0.2.2,
>
>
> \[2\] Zhang, Q., Kong, D., Shi, P., Singh, V.P., Sun, P., 2018.
> Vegetation phenology on the Qinghai-Tibetan Plateau and its response
> to climate change (1982–2013). Agric. For. Meteorol. 248, 408–417.
>
# Acknowledgements
Keep in mind that this repository is released under a GPL2 license,
which permits commercial use but requires that the source code (of
derivatives) is always open even if hosted as a web service.