Statistical Trend Analysis (STA) for Time Series of Satellite Imagery

Description

STA applies the Mann-Kendall test for trend to the so-called shape parameters of periodic time series. STA estimates shape parameters via harmonic regression. STA can handle numeric time series and RasterStack of satellite images.

Details

Shape parameters is the term used in vegetation monitoring to refer to the amplitudes and phase angles resulting from fitting a harmonic regression model to time series of vegetation indices derived from satellite images. Regardless of its origin, STA can be applied to any periodic time series which makes this package potentially useful to other disciplines such as hydrology, climatology and econometrics.

With sta (the main function of this package) it is possible to perform the Mann-Kendall test for trend on time series of the three most commonly used shape parameters: mean, annual and semiannual. These parameters are the estimated amplitude coefficients of the aforementioned harmonic regresion model. This function allows parallel processing to handle large satellite time series imagery.

STA includes the following graphical methods:

• plot.staNumeric: generic plot displaying sta's output for numeric time series.

• plot.staMatrix: maps of mapview-class displaying sta's output for RasterStack.

STA include the following datasets:

• marismas: numeric vector containing 10-day Composite NDMI values from 2000 to 2018.

• ndmi: RasterStack containing 612 spatial subsets of 10-day Composite NDMI images acquired from 2001 to 2017.

Author(s)

Tecuapetla-Gómez, I. itecuapetla@conabio.gob.mx

References

Eastman, R., Sangermano, F., Ghimire, B., Zhu, H., Chen, H., Neeti, N., Cai, Y., Machado, E., Crema, S. (2009). Seasonal trend analysis of image time series, International Journal of Remote Sensing 30(10), 2721–2726.

Get a RasterLayer with no missing values from a Raster* object

Description

The term master refers to a raster layer whose extent and coordinate reference system are used as a reference to rasterize further objects, e.g. matrices. To rasterize, master must be free of missing values.

Usage

getMaster(x)

Arguments

Value

RasterLayer

See Also

matrixToRaster

10-day Composite NDMI Time Series

Description

A numeric vector of length 684 containing 10-day composite values of the Normalized Difference Moisture Index (NDMI) from 2000 to 2018. NDMI = (NIR-MIR)/(NIR+MIR) where NIR and MIR are the Near Infrared and Mid-Infrared bands of the MCD43A4 MODIS product, respectively. This numeric vector was taken from a RasterStack covering the Natural Protected Area Reserva de la Biosfera Marismas Nacionales at Nayarit, Mexico.

Usage

data(marismas)

Format

An object of class "numeric".

10-day Composite NDMI RasterStack

Description

A RasterStack containing 612 layers of the Normalized Difference Moisture Index (NDMI) from 2001 to 2017. NDMI = (NIR-MIR)/(NIR+MIR) where NIR and MIR are the Near Infrared and Mid-Infrared bands of the MCD43A4 MODIS product, respectively. This RasterStack is a spatial subset of a larger RasterStack covering the Natural Protected Area Reserva de la Biosfera Marismas Nacionales at Nayarit, Mexico.

ndmi.tif

A "RasterStack" object with 36 rows, 55 columns, 1980 cells and 612 layers.

Plot method for sta function

Description

This function displays some maps of mapview-class

Usage

## S3 method for class 'staMatrix'
plot(x, significance = NULL, master, ...)

Arguments

See Also

sta, getMaster, matrixToRaster

Plot method for sta function

Description

This function returns a plot

Usage

## S3 method for class 'staNumeric'
plot(x, significance = NULL, ...)

Arguments

See Also

sta

Statistical Seasonal Trend Analysis for numeric vector or RasterStack

Description

Statistical Seasonal Trend Analysis for numeric vector or RasterStack

Usage

sta(
  data,
  freq,
  numFreq = 4,
  delta = 0,
  startYear = 2000,
  endYear = 2018,
  intraAnnualPeriod = c("wetSeason", "drySeason"),
  interAnnualPeriod,
  adhocPeriod = NULL,
  significance = NULL,
  save = FALSE,
  dirToSaveSTA = NULL,
  numCores = 20
)

Arguments

Details

When the input is a matrix, its first two columns must correspond to geographic coordinates. For instance, the matrix resulting from applying rasterToPoints to a RasterStack has this format.

freq must be either 12 (monthly observations), 23 (Landsat annual scale) or 36 (10-day composite) as this version implements STA for time series with these frequencies.

This version sets intraAnnualPeriod to either the wetSeason or the drySeason of Mexico. Empirical evidence suggests that while wet season runs from May to October, dry season runs from November to April. Should a desired STA require specific months/days, these must be provided through adhocPeriod.

When interAnnualPeriod is not specified and class(data)=numeric, interAnnualPeriod = 1:(length(data)/freq); when class(data) is either RasterStack or matrix, interAnnualPeriod = 1:((ncol(data)-2)/freq).

Since adhocPeriod defines an inter annual period "ad-hoc", the specific days of this ad-hoc season must be known in advance and consequently, the specific time-points (with respect to the time series under consideration) must be provided in a numeric vector.

When save=T, a valid dirToSaveSTA must be provided, that is, this folder should have been created previously. In this case, sta's output is saved on dirToSaveSTA. This version saves arrays of STA of the mean, annual and semi-annual parameters (along with their corresponding basic statistics) in the file sta_matrix_output.RData inside dirToSaveSTA. Also, in the same directory, the file sta_progress.txt records the progress of the STA process.

save=T, dirToSaveSTA, numCores and master are required when data is either a RasterStack or a matrix. The aforementioned basic statistics are: mean and standard deviation of the time series of annual maximum and minimum as well as the global minima and maxima.

Value

When class(data) is a numeric, an object of class "staNumeric" containing:

• mean a list of 12 elements with STA output for shape parameter mean

• annual a list of 12 elements with STA output for shape parameter annual

• semiannual a list of 12 elements with STA output for shape parameter semiannual

When class(data) is a RasterStack or a matrix, an object of class "staMatrix" containing:

• mean a matrix of 4 columns with STA output for shape parameter mean. First two columns provide geolocation of analyzed pixels, third and fourth columns show p-value and slope estimate of STA, respectively

• mean_basicStats a matrix of 10 columns with basic statistics for shape parameter mean. First two columns provide geolocation of analyzed pixels, from third to tenth columns show mean, standard deviation, global minimum, and maximum of minimum values, as well as mean, standard deviation, global minimum, and maximum of maximum values, respectively

• annual a matrix of 4 columns with STA output for shape parameter annual. First two columns provide geolocation of analyzed pixels, third and fourth columns show p-value and slope estimate of STA, respectively

• annual_basicStats a matrix of 10 columns with basic statistics for shape parameter annual. First two columns provide geolocation of analyzed pixels, from third to tenth columns show mean, standard deviation, global minimum, and maximum of minimum values, as well as mean, standard deviation, global minimum, and maximum of maximum values, respectively

• semiannual a matrix of 4 columns with STA output for shape parameter semiannual. First two columns provide geolocation of analyzed pixels, third and fourth columns show p-value and slope estimate of STA, respectively

• semiannual_basicStats a matrix of 10 columns with basic statistics for shape parameter semiannual. First two columns provide geolocation of analyzed pixels, from third to tenth columns show mean, standard deviation, global minimum, and maximum of minimum values, as well as mean, standard deviation, global minimum, and maximum of maximum values, respectively

Note

STA is based on the following ideas. Let y_t denote the value of a periodic time series at time-point t. It is well-known that this type of observations can be modeled as:

y_t = a_0 + a_1 cos( (2 \pi t)/L - \phi_1) + ... + a_K cos( (2 \pi K t)/L - \phi_K) + \varepsilon_t, t=1,\ldots,L.

This model is known as harmonic regression. L denotes the number of observations per period, K is the number of harmonics included in the fit, a_i's and \phi_i's are called amplitude coefficients and phase angles, respectively. K can be known empirically. Amplitudes and phase angle can be obtained as the solution of a least-squares problem.

In vegetation monitoring, amplitudes and phase angles are known as shape parameters. In particular, a_0, a_1 and a_2 are called mean and annual and semiannual amplitudes, respectively. Applying the harmonic regression model to observations over P periods of length L each, results in estimates of shape parameters for each period. Thus, focusing only on amplitudes, sta yields time series of mean, annual and semiannual parameters. A subsequent Mann-Kendall test for trend is performed on each of these series.

References

Eastman, R., Sangermano, F., Ghimine, B., Zhu, H., Chen, H., Neeti, N., Cai, Y., Machado E., Crema, S. (2009). Seasonal trend analysis of image time series, International Journal of Remote Sensing, 30(10), 2721–2726.

Examples

sta_marismas <- sta(data=marismas, freq=36)
str(sta_marismas)
plot(sta_marismas)
plot(sta_marismas, significance=0.09)

# Use of interAnnualPeriod
sta_21016 <- sta(data = marismas, freq = 36, interAnnualPeriod = c(2, 10, 16))
plot(sta_21016)

# Use of intraAnnualPeriod
sta_drySeason_218 <- sta(data = marismas, freq = 36,
                     interAnnualPeriod = 2:18, intraAnnualPeriod = "drySeason")
plot(sta_drySeason_218)

# Use of adhocPeriod and significance
adhoc <- list()
beginPeriod <- (1:17) * 36
endPeriod <- 2:18 * 36 
adhoc$partial <- c( sapply(1:length(beginPeriod), 
                 function(s) c(beginPeriod[s]+1, endPeriod[s]) ) )
adhoc$full <- c( sapply(1:length(beginPeriod), 
              function(s) (beginPeriod[s]+1):endPeriod[s]) )
sta_adhoc_218 <- sta(data = marismas, freq = 36, interAnnualPeriod = 2:18,
                 startYear = 2000, endYear = 2018, adhocPeriod = adhoc, significance=0.05)
plot(sta_adhoc_218)

# Use of ndmi RasterStack

ndmi_path = system.file("extdata", "ndmi.tif", package = "sta")
ndmiSTACK <- stack(ndmi_path)
dir.create(path=paste0(system.file("extdata", package="sta"), "/output_ndmi"),
          showWarnings=FALSE)
outputDIR = paste0(system.file("extdata", package="sta"), "/output_ndmi")

sta_ndmi_21016 <- sta(data = ndmiSTACK, freq = 36,
                  numFreq = 4, delta = 0.2, intraAnnualPeriod = "wetSeason",
                  startYear = 2000, endYear = 2018, interAnnualPeriod = c(2,10,16),
                  save = TRUE, numCores = 2L, dirToSaveSTA = outputDIR)