Extract vegetation phenology from MOD13A1 EVI

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.

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_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

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.

d$Tn %<>% zoo::na.approx(maxgap = 4)
plot(d$Tn, type = "l"); abline(a = 5, b = 0, col = "red")

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 -167. 0.169 0.399 0.0195  4.34 0.0175  6   
##  2 2000_1  202. 0.168 0.407 0.0320  3.56 0.0152  5.79
##  3 2001_1  564. 0.168 0.402 0.0204  3.92 0.0169  6   
##  4 2002_1  934. 0.168 0.530 0.0327  2    0.0156  2   
##  5 2003_1 1275. 0.167 0.433 0.0268  2    0.0124  3.64
##  6 2004_1 1689. 0.169 0.445 0.0166  5.09 0.0372  2   
##  7 2005_1 2025. 0.172 0.450 0.0223  6    0.0162  3.62
##  8 2006_1 2370. 0.166 0.420 0.0280  2.01 0.0112  3.07
##  9 2007_1 2752. 0.166 0.483 0.0211  2    0.0150  2.90
## 10 2008_1 3153. 0.170 0.481 0.0132  5.09 0.0301  6   
## 11 2009_1 3523. 0.168 0.479 0.0141  6    0.0282  2   
## 12 2010_1 3840. 0.168 0.488 0.0231  2    0.0129  2   
## 13 2011_1 4209. 0.175 0.468 0.0283  2    0.0135  4.92
## 14 2012_1 4586. 0.166 0.523 0.0201  4.61 0.0156  2   
## 15 2013_1 4966. 0.167 0.482 0.0142  6    0.0188  2   
## 16 2014_1 5306. 0.169 0.507 0.0281  2    0.0127  2.70
## 17 2015_1 5700. 0.183 0.485 0.0146  5.40 0.0277  2   
## 18 2016_1 6024. 0.175 0.487 0.0363  2    0.0124  4.38
## 19 2017_1 6405. 0.169 0.448 0.0203  5.99 0.0128  3.61
## 20 2018_1 6789. 0.165 0.450 0.0121  4.11 0.0208  2
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.
## Get GOF information
d_gof <- get_GOF(fit)
# fit$stat <- stat
print(head(d_gof))
##      flag   meth        R2       NSE         R       RMSE       pvalue n_sim
##    <char> <char>     <num>     <num>     <num>      <num>        <num> <num>
## 1: 1999_1     AG 0.8233472 0.6145234 0.9073848 0.08520015 2.337043e-09    23
## 2: 1999_1  Zhang 0.8282894 0.6183996 0.9101041 0.08477070 1.730224e-09    23
## 3: 1999_1   Beck 0.8286265 0.6186920 0.9102893 0.08473822 1.694577e-09    23
## 4: 1999_1 Elmore 0.8334801 0.6051443 0.9129513 0.08623043 1.249992e-09    23
## 5: 2000_1     AG 0.8189371 0.5535112 0.9049514 0.09770633 7.301265e-09    22
## 6: 2000_1  Zhang 0.8263049 0.5556140 0.9090131 0.09747597 4.799998e-09    22
# 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      167      275      175      264      177      261
## 3: 2001_1 2001-01-01      143      263      154      255      157      252
## 4: 2002_1 2002-01-01      165      286      178      258      182      250
## 5: 2003_1 2003-01-01      132      271      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     199     254   136   168   239   267     129      175        235
## 2:     174     203     266   163   186   249   281     156      193        243
## 3:     153     199     256   138   171   241   269     130      179        236
## 4:     182     204     249   161   196   222   295     153      204        220
## 5:     153     180     254   127   171   227   282     118      180        197
## 6:     171     229     248   157   189   235   268     150      196        274
##    Dormancy
##       <num>
## 1:      272
## 2:      287
## 3:      275
## 4:      310
## 5:      294
## 6:      375