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.

library(data.table)
library(lubridate)
library(purrr)
library(ggplot2)
library(phenofit)

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`.

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)]

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(FUN = 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))
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

# 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)


Try the phenofit package in your browser

Any scripts or data that you put into this service are public.

phenofit documentation built on Oct. 15, 2021, 5:09 p.m.