Extract vegetation phenology from MOD13A1 EVI
Source:vignettes/phenofit_CA-NS6.Rmd
phenofit_CA-NS6.RmdExample
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.
Set global parameters for phenofit
# 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 |
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_HeadTailis 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
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.
2.2 Check input data
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.
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
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))## List of 4
## $ AG : tibble [20 × 8] (S3: tbl_df/tbl/data.frame)
## $ Zhang : tibble [20 × 8] (S3: tbl_df/tbl/data.frame)
## $ Beck : tibble [20 × 7] (S3: tbl_df/tbl/data.frame)
## $ Elmore: tibble [20 × 8] (S3: tbl_df/tbl/data.frame)
## NULL
print(l_param$AG)## # A tibble: 20 × 8
## flag t0 mn mx rsp a3 rau a5
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1999_1 -166. 0.169 0.398 0.0193 4.43 0.0177 6
## 2 2000_1 209. 0.167 0.407 0.0255 4.36 0.0171 4.84
## 3 2001_1 567. 0.169 0.398 0.0193 4.43 0.0177 6
## 4 2002_1 934. 0.168 0.528 0.0329 2 0.0161 2
## 5 2003_1 1275. 0.167 0.432 0.0268 2 0.0123 3.73
## 6 2004_1 1689. 0.168 0.446 0.0167 5.04 0.0369 2
## 7 2005_1 2025. 0.172 0.449 0.0224 6 0.0161 3.62
## 8 2006_1 2370. 0.166 0.420 0.0280 2.02 0.0112 3.12
## 9 2007_1 2752. 0.166 0.483 0.0210 2 0.0150 2.90
## 10 2008_1 3123. 0.170 0.481 0.0215 2.74 0.0165 6
## 11 2009_1 3522. 0.168 0.479 0.0142 6 0.0277 2
## 12 2010_1 3840. 0.167 0.489 0.0228 2 0.0131 2
## 13 2011_1 4207. 0.175 0.467 0.0290 2 0.0133 5.04
## 14 2012_1 4558. 0.166 0.506 0.0483 2 0.0109 4.15
## 15 2013_1 4965. 0.166 0.483 0.0144 6 0.0190 2.47
## 16 2014_1 5306. 0.169 0.507 0.0281 2 0.0127 2.70
## 17 2015_1 5701. 0.184 0.486 0.0145 5.62 0.0280 2
## 18 2016_1 6024. 0.176 0.484 0.0370 2 0.0126 5.14
## 19 2017_1 6384. 0.168 0.445 0.0367 3.09 0.0101 5.37
## 20 2018_1 6774. 0.165 0.451 0.0155 3.22 0.0164 3.07
d_fit <- get_fitting(fit)## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.
## Warning in merge.data.table(x$data[Ix], df, id = "t"): Unknown argument 'id'
## has been passed.## flag meth R2 NSE R RMSE pvalue n_sim
## <char> <char> <num> <num> <num> <num> <num> <num>
## 1: 1999_1 AG 0.8235796 0.6150465 0.9075129 0.08514232 2.304668e-09 23
## 2: 1999_1 Zhang 0.8285788 0.6189771 0.9102630 0.08470653 1.699585e-09 23
## 3: 1999_1 Beck 0.8289473 0.6192976 0.9104654 0.08467090 1.661281e-09 23
## 4: 1999_1 Elmore 0.8334679 0.6054137 0.9129446 0.08620101 1.250963e-09 23
## 5: 2000_1 AG 0.8252408 0.5496414 0.9084277 0.09660186 2.084845e-09 23
## 6: 2000_1 Zhang 0.8271460 0.5499378 0.9094757 0.09657007 1.856258e-09 23
# print(fit$fits$AG$`2002_1`$ws)
print(fit$`2002_1`$fFIT$AG$ws)## NULL
## visualization
g <- plot_curvefits(d_fit, brks2, NULL, ylab = "NDVI", "Time",
theme = coord_cartesian(xlim = c(ymd("2000-04-01"), ymd("2017-07-31"))))## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_point()`).
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
# 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)## flag origin TRS2.sos TRS2.eos TRS5.sos TRS5.eos TRS6.sos TRS6.eos
## <char> <Date> <num> <num> <num> <num> <num> <num>
## 1: 1999_1 1999-01-01 142 261 152 253 155 250
## 2: 2000_1 2000-01-01 166 275 174 264 177 261
## 3: 2001_1 2001-01-01 143 263 154 255 157 252
## 4: 2002_1 2002-01-01 164 284 178 256 182 248
## 5: 2003_1 2003-01-01 132 273 149 253 153 247
## 6: 2004_1 2004-01-01 162 264 173 252 176 248
## DER.sos DER.pos DER.eos UD SD DD RD Greenup Maturity Senescence
## <num> <num> <num> <num> <num> <num> <num> <num> <num> <num>
## 1: 151 200 254 136 167 240 267 129 175 235
## 2: 173 210 266 162 186 248 281 156 192 241
## 3: 153 202 256 138 169 242 269 131 177 237
## 4: 182 204 248 161 196 222 292 153 204 219
## 5: 153 180 254 127 170 228 282 118 180 197
## 6: 172 228 248 157 189 235 268 150 196 274
## Dormancy
## <num>
## 1: 272
## 2: 288
## 3: 274
## 4: 307
## 5: 294
## 6: 377