Metrics of land surface phenology and greenness are often derived from vegetation index time series to map the spatial variability of vegetation phenology or to analyse temporal changes in vegetation. Phenology and greenness metrics (PGMs) are for example the start and end of growing season (SOS and EOS) or the greenness during the peak of the growing season (PEAK). A description of these metrics and a discussion on the caveats and challenges of phenology detection from remote sensing-derived vegetation index time series can be found in Forkel et al. (2015):
In the greenbrown package, phenology detection is implemented as a step-wise procedure with prepocessing and analysis steps as described in Forkel et al. (2015). The steps are:
The corresponding functions to these steps can be used separately or the entire step-wise procedure can be applied just by using the main function Phenology for time series and PhenologyRaster for multi-temporal raster data.
The function Seasonality is used during several steps to check if the time series has seasonality. Otherwise, the calculation of most phenology metrics like SOS and EOS is not useful.
The main functions for phenology detection that cover these three steps are Phenology for time series and PhenologyRaster for multi-temporal raster data. These two functions encapsulate several other functions that are needed during the different steps and that are explained below.
> # load example NDVI time series
> # compute phenology metrics
> spl.trs <- Phenology(ndvi, tsgf="TSGFspline", approach="White")
--- Phenology ----------------------------
Calculate phenology metrics on time series
- smoothing and gap filling : TSGFspline
- summary approach : White
Mean and standard deviation of annual metrics:
SOS : 86 +- 11
EOS : 280 +- 15
LOS : 190 +- 21
MSP : 0.27 +- 0.019
MAU : 0.27 +- 0.019
RSP : NaN +- NA
RAU : NaN +- NA
POP : 180 +- 10
POT : 0.37 +- 1.3
MGS : 0.3 +- 0.021
PEAK : 0.32 +- 0.022
TROUGH : 0.23 +- 0.018
Fig. 1: Annual time series of the the start of growing season (sos), end of growing season (eos) and position of peak day (pop).
The phenology analysis can be easily applied to raster datasets with PhenologyRaster.
> # load example NDVI time series
> phenmap <- PhenologyRaster(ndvimap, start=c(1982, 1), freq=12, tsgf="TSGFspline", approach="Deriv")
> plot(phenmap, c(grep("SOS.1982", names(phenmap)), grep("EOS.1982", names(phenmap))))
Fig. 2: Maps of start and end of the growing season in 1982 (in day of year).
Permanent gaps in time series are missing values that occur regulary every year or in many years. Permanent gaps in satellite-derived vegetation index time series are for example the winter months in northern regions where reliable observation are not possible because of high sun zenit angles or snow cover. Such gaps can be identified using the function IsPermanentGap and these gaps can be filled with FillPermanentGaps.
> # set NA values into winter months to simulate gaps
> winter <- (1:length(ndvi))[cycle(ndvi) == 1 | cycle(ndvi) == 2 | cycle(ndvi) == 12]
> gaps <- sample(winter, length(winter)*0.3)
> ndvi2 <- ndvi
> ndvi2[gaps] <- NA
> # check and fill permanent gaps
> fill <- FillPermanentGaps(ndvi2) # default fills with the minimum value
> plot(fill, col="red"); lines(ndvi)
Several functions are implemented to smooth, to gap-fill and interpolate time series to daily time steps:
The functions of steps 1 and 2 are both encapsulated by the function TsPP (time series preprocessing).
> # introduce random gaps
> gaps <- ndvi
> gaps[runif(100, 1, length(ndvi))] <- NA
> # smoothing and gap filling
> tsgf <- TSGFspline(gaps)
> plot(gaps, main="Temporal smoothing and gap-filling using splines")
> lines(tsgf, col="red")
Fig. 3: Original time series with gaps (black) and smoothed and interpolated time series (red).
In the last step, phenology events like the start of growing season (SOS) and end of growing season (EOS) are estimated from the daily interpolated, gap-filled and smoothed time series. This can be done either by using thresholds (PhenoTrs) or by extreme values of the derivative of the seasonal cycle (PhenoDeriv). PhenoTrs implements thereby two approaches to estimate phenology events.
> # time series pre-processing
> x <- TsPP(ndvi, interpolate=TRUE)[1:365]
> # calculate phenology metrics for first year
> PhenoDeriv(x, plot=TRUE)
Fig. 4: Seasonal cycle of a single year with estimated PGMs (phenology and greenness metrics).
greenbrown, Matthias Forkel, 2016-11-17