forked from eco-hydro/phenofit
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathREADME.Rmd
249 lines (203 loc) · 9.85 KB
/
README.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
---
output: github_document
---
```{r, echo = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#",
fig.width = 10, fig.height = 5,
fig.align = "center",
fig.path = "man/Figure/",
dev = 'svg'
)
```
# 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
*<u>Figure 1. The flowchart of phenology extraction in `phenofit`.</u>*
# Installation
You can install phenofit from github with:
```{r gh-installation, eval = FALSE}
# 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 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}
suppressMessages({
library(data.table)
library(magrittr)
library(lubridate)
library(purrr)
library(plyr)
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,
FUN = smooth_wWHIT, wFUN = wFUN,
maxExtendMonth = 6, r_min = 0.1,
IsPlot = IsPlot, IsPlot.OnlyBad = FALSE, print = print)
```
## 2.4 Curve fitting
```{r curve fitting, fig.height=7, fig.align="center"}
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))
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_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
# 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)
```
# **References**
> [1\] Dongdong Kong, R package: A state-of-the-art Vegetation Phenology extraction package, `phenofit` version 0.2.2, <https://github.com/kongdd/phenofit>
>
> [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. <https://doi.org/10.1016/j.agrformet.2017.10.026>
# 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.