The irg package opts for a tabular calculation of the instantaneous rate of green-up (IRG) as opposed to a raster based approach. Sampling MODIS imagery is left up to the user and a prerequisite for all functions. The main input (DT) for all functions is a data.table of an NDVI time series. The sampling unit (id) is flexible (a decision for the user) though we would anticipate points or polygons, or maybe a pixel. All functions leverage the speed of data.table to efficiently filter, scale and model NDVI time series, and calculate IRG.

Setup

Install the latest version with remotes.

remotes::install_gitlab('robit.a/irg')

Packages

irg depends on two packages (and stats):

No external dependencies.

Input data

irg requires an NDVI time series in a data.table.

Though names can be different and specified at input, the default names and required columns are:

  • id: sampling unit identification (point id, polygon id, …)
  • yr: year of sample
  • DayOfYear: julian day
  • NDVI: sampled NDVI value
  • SummaryQA: Quality reliability of pixel

SummaryQA details:

  • 0: Good data, use with confidence
  • 1: Marginal data, useful but look at detailed QA for more information
  • 2: Pixel covered with snow/ice
  • 3: Pixel is cloudy

Let’s take a look at the example data.

library(irg)
library(data.table)

ndvi <- fread(system.file('extdata', 'ndvi.csv', package = 'irg'))

# or look at the help page
?ndvi
id yr DayOfYear NDVI SummaryQA
5 2002 284 6495 0
1 2002 299 5517 1
2 2002 297 6352 1
3 2002 289 5725 0
4 2002 299 5990 1
5 2002 299 5950 1

data.frame

If your data is a data.frame, convert it by reference:

# Pretend
DF <- as.data.frame(ndvi)

# Convert by reference
setDT(DF)

Sampling NDVI

Though irg is not involved in the sampling step, it is important that the input data matches the package’s expectations.

We used the incredible Google Earth Engine to sample MODIS NDVI (MOD13Q1.006). There are also R packages specific to MODIS (MODIStsp) and general purpose raster operations (raster), and others (let us know)…

Temporal extent

Filtering steps in irg use a baseline ‘winterNDVI’ and upper quantile as described by Bischoff et al. (2012). These steps require multiple years of sampled NDVI for each id. In addition, make sure to include all samples throughout the year, leaving the filtering for irg.

Workflow

There are 5 filtering functions, 2 scaling functions, 3 modeling functions and 2 IRG functions.

The irg::irg function is a wrapper for all steps - filtering, scaling, modeling and calculating IRG in one step. At this point, only defaults. Here’s 5 rows from the result.

For options, head to the steps below.

out <- irg(ndvi)
#> Warning in calc_irg(model_ndvi(DT, observed = FALSE)): NAs found in DT, IRG will
#> be set to NA.
#> Warning in min(irg, na.rm = TRUE): no non-missing arguments to min; returning
#> Inf
#> Warning in max(irg, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(irg, na.rm = TRUE): no non-missing arguments to min; returning
#> Inf

#> Warning in min(irg, na.rm = TRUE): no non-missing arguments to min; returning
#> Inf
#> Warning in max(irg, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(irg, na.rm = TRUE): no non-missing arguments to min; returning
#> Inf

#> Warning in min(irg, na.rm = TRUE): no non-missing arguments to min; returning
#> Inf
#> Warning in max(irg, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(irg, na.rm = TRUE): no non-missing arguments to min; returning
#> Inf

#> Warning in min(irg, na.rm = TRUE): no non-missing arguments to min; returning
#> Inf
#> Warning in max(irg, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(irg, na.rm = TRUE): no non-missing arguments to min; returning
#> Inf

#> Warning in min(irg, na.rm = TRUE): no non-missing arguments to min; returning
#> Inf
#> Warning in max(irg, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(irg, na.rm = TRUE): no non-missing arguments to min; returning
#> Inf

#> Warning in min(irg, na.rm = TRUE): no non-missing arguments to min; returning
#> Inf
#> Warning in max(irg, na.rm = TRUE): no non-missing arguments to max; returning
#> -Inf
#> Warning in min(irg, na.rm = TRUE): no non-missing arguments to min; returning
#> Inf
id yr t fitted irg
2 2003 0.4000000 0.1504764 0.5111129
2 2003 0.4027397 0.1568717 0.5288405
2 2003 0.4054795 0.1634865 0.5468322
2 2003 0.4082192 0.1703238 0.5650611
2 2003 0.4109589 0.1773864 0.5834979

Filtering

There are 5 filtering functions.

functions arguments
filter_ndvi DT
filter_qa DT, qa, good
filter_roll DT, window, id, method
filter_top DT, probs, id
filter_winter DT, probs, limits, doy, id
# Load data.table
library(data.table)
library(irg)

# Read in example data
ndvi <- fread(system.file('extdata', 'ndvi.csv', package = 'irg'))

# Filter NDVI time series
filter_qa(ndvi, qa = 'SummaryQA', good = c(0, 1))

filter_winter(ndvi, probs = 0.025, limits = c(60L, 300L),
                            doy = 'DayOfYear', id = 'id')

filter_roll(ndvi, window = 3L, id = 'id', method = 'median')

filter_top(ndvi, probs = 0.925, id = 'id')

Scaling

Two scaling functions are use to scale the day of year column and filtered NDVI time series between 0-1.

# Scale variables
scale_doy(ndvi, doy = 'DayOfYear')
scale_ndvi(ndvi)

Modeling

Three functions are used to model the NDVI times series to a double logistic curve, as described by Bischoff et al. (2012).

\[fitted = \frac{1}{1 + e^ \frac{xmidS - t}{scalS}} - \frac{1}{1 + e^ \frac{xmidA - t}{scalA}}\]

Two options from this point are available: fitting NDVI and calculating IRG for observed data only, or for the full year.

To calculate for every day of every year, specify returns = 'models' in model_params, observed = FALSE in model_ndvi and assign the output of model_ndvi.

# Guess starting parameters
model_start(ndvi, id = 'id', year = 'yr')

# Double logistic model parameters given starting parameters for nls
mods <- model_params(
  ndvi,
  returns = 'models',
  id = 'id', year = 'yr',
  xmidS = 'xmidS_start', xmidA = 'xmidA_start',
  scalS = 0.05,
  scalA = 0.01
)

# Fit double log to NDVI
fit <- model_ndvi(mods, observed = FALSE)

Alternatively, to calculate for the observed data only, specify returns = 'columns' in model_params and observed = TRUE in model_ndvi.

# Guess starting parameters
model_start(ndvi, id = 'id', year = 'yr')

# Double logistic model parameters given starting parameters for nls
model_params(
  ndvi,
  returns = 'columns',
  id = 'id', year = 'yr',
  xmidS = 'xmidS_start', xmidA = 'xmidA_start',
  scalS = 0.05,
  scalA = 0.01
)

# Fit double log to NDVI
 model_ndvi(ndvi, observed = TRUE)

IRG

\[IRG = \frac{e ^ \frac{t + xmidS}{scalS}}{2 scalS e ^ \frac{t + xmidS}{scalS} + scalS e ^ \frac{2t}{scalS} + scalS e ^ \frac{2midS}{scalS}}\] Finally, calculate IRG:

# Calculate IRG for each day of the year
calc_irg(fit)

# or for observed data
calc_irg(ndvi)