Title: | Seasonal Trend Analysis for Time Series Imagery in R |
---|---|
Description: | Efficiently estimate shape parameters of periodic time series imagery with which a statistical seasonal trend analysis (STA) is subsequently performed. STA output can be exported in conventional raster formats. Methods to visualize STA output are also implemented as well as the calculation of additional basic statistics. STA is based on (R. Eastman, F. Sangermano, B. Ghimire, H. Zhu, H. Chen, N. Neeti, Y. Cai, E. Machado and S. Crema, 2009) <doi:10.1080/01431160902755338>. |
Authors: | Inder Tecuapetla-Gomez [aut, cre] |
Maintainer: | Inder Tecuapetla-Gomez <[email protected]> |
License: | GPL (>= 2) |
Version: | 0.1.7 |
Built: | 2024-12-11 07:25:47 UTC |
Source: | CRAN |
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.
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.
Tecuapetla-Gómez, I. [email protected]
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.
RasterLayer
with no missing values from a Raster*
objectThe 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.
getMaster(x)
getMaster(x)
x |
|
RasterLayer
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.
data(marismas)
data(marismas)
An object of class "numeric"
.
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.
A "RasterStack"
object with 36 rows, 55 columns, 1980 cells and 612 layers.
sta
functionThis function displays some maps of mapview-class
## S3 method for class 'staMatrix' plot(x, significance = NULL, master, ...)
## S3 method for class 'staMatrix' plot(x, significance = NULL, master, ...)
x |
an object of class "staMatrix" |
significance |
numeric indicating significance of each shape parameter trend |
master |
RasterLayer used to transfer STA output to raster layers. |
... |
additional plot parameters |
sta
, getMaster
, matrixToRaster
sta
functionThis function returns a plot
## S3 method for class 'staNumeric' plot(x, significance = NULL, ...)
## S3 method for class 'staNumeric' plot(x, significance = NULL, ...)
x |
an object of class "staNumeric" |
significance |
numeric indicating significance of each shape parameter trend |
... |
additional plot parameters |
RasterStack
Statistical Seasonal Trend Analysis for numeric vector or RasterStack
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 )
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 )
data |
numeric vector, matrix or RasterStack object |
freq |
integer with the number of observations per period. See |
numFreq |
integer with the number of frequencies to employ in harmonic regression fitting.
See |
delta |
numeric (positive) controlling regularization and prevent non-invertible
hat matrix in harmonic regression model. See |
startYear |
numeric, time series initial year |
endYear |
numeric, time series last year |
intraAnnualPeriod |
character indicating seasons (wet or dry) to be considered for additional
statistical analysis. See |
interAnnualPeriod |
numeric vector indicating the number of years to be considered in STA. For instance,
1:5, indicates that the first five years will be
utilized for STA. Similarly, c(2,6,10) indicates that the second, sixth and tenth
years will be utilized for STA. See |
adhocPeriod |
numeric vector with the specific observations to be considered in additional
statistical analysis. See |
significance |
numeric in the interval (0,1) to assess statistical significance of trend analysis.
|
save |
logical, should STA output be saved, default is |
dirToSaveSTA |
character with full path name used to save |
numCores |
integer given the number of cores to use; pertinent when |
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.
When class(data)
is a numeric
, an object of class "staNumeric" containing:
data |
numeric vector |
freq |
integer with the number of observations per period |
startYear |
numeric, time series initial year |
endYear |
numeric, time series last year |
intraAnnualPeriod |
character indicating seasons (wet or dry) |
interAnnualPeriod |
numeric vector indicating the number of years considered in STA |
ts |
time series object; |
fit |
numeric vector with output of |
sta |
a list containing the following elements: |
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
significance |
numeric in the interval (0,1) or |
When class(data)
is a RasterStack
or a matrix
, an object of class
"staMatrix" containing:
freq |
integer with the number of observations per period |
daysUsedFit |
integer vector indicating the indices used to fit harmonic regression. see |
intervalsUsedBasicStats |
numeric vector indicating the indices used on calculation of basic statistics |
sta |
a list containg the following elements: |
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
STA is based on the following ideas. Let denote the value of a periodic time
series at time-point
. It is well-known that this type of observations can be modeled
as:
,
.
This model is known as harmonic regression. denotes the number of observations per period,
is the number of
harmonics included in the fit,
's and
's are called amplitude coefficients and phase angles,
respectively.
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,
,
and
are called mean and annual and semiannual amplitudes, respectively.
Applying the harmonic regression model to observations over
periods of length
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.
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.
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)
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)