bfastmonitor {bfast}  R Documentation 
Near RealTime Disturbance Detection Based on BFASTType Models
Description
Monitoring disturbances in time series models (with trend/season/regressor terms) at the end of time series (i.e., in near realtime). Based on a model for stable historical behaviour abnormal changes within newly acquired data can be detected. Different models are available for modeling the stable historical behavior. A seasontrend model (with harmonic seasonal pattern) is used as a default in the regresssion modelling.
Usage
bfastmonitor(
data,
start,
formula = response ~ trend + harmon,
order = 3,
lag = NULL,
slag = NULL,
history = c("ROC", "BP", "all"),
type = "OLSMOSUM",
h = 0.25,
end = 10,
level = c(0.05, 0.05),
hpc = "none",
verbose = FALSE,
plot = FALSE,
sbins = 1
)
Arguments
data 
A time series of class 
start 
numeric. The starting date of the monitoring period. Can either
be given as a float (e.g., 
formula 
formula for the regression model. The default is

order 
numeric. Order of the harmonic term, defaulting to 
lag 
numeric. Order of the autoregressive term, by default omitted. 
slag 
numeric. Order of the seasonal autoregressive term, by default omitted. 
history 
specification of the start of the stable history period. Can
either be a character, numeric, or a function. If character, then selection
is possible between reverseordered CUSUM ( 
type 
character specifying the type of monitoring process. By default,
a MOSUM process based on OLS residuals is employed. See

h 
numeric scalar from interval (0,1) specifying the bandwidth relative to the sample size in MOSUM/ME monitoring processes. 
end 
numeric. Maximum time (relative to the history period) that will be monitored (in MOSUM/ME processes). Default is 10 times the history period. 
level 
numeric vector. Significance levels of the monitoring and ROC (if selected) procedure, i.e., probability of type I error. 
hpc 
character specifying the high performance computing support.
Default is 
verbose 
logical. Should information about the monitoring be printed during computation? 
plot 
logical. Should the result be plotted? 
sbins 
numeric. Number of seasonal dummies, passed to

Details
bfastmonitor
provides monitoring of disturbances (or structural
changes) in near realtime based on a wide class of time series regression
models with optional season/trend/autoregressive/covariate terms. See
Verbesselt at al. (2011) for details.
Based on a given time series (typically, but not necessarily, with frequency
greater than 1), the data is first preprocessed for regression modeling.
Trend/season/autoregressive/covariate terms are (optionally) computed using
bfastpp
. Second, the data is split into a history and
monitoring period (starting with start
). Third, a subset of the
history period is determined which is considered to be stable (see also
below). Fourth, a regression model is fitted to the preprocessed data in
the stable history period. Fifth, a monitoring procedure is used to
determine whether the observations in the monitoring period conform with
this stable regression model or whether a change is detected.
The regression model can be specified by the user. The default is to use a
linear trend and a harmonic season: response ~ trend + harmon
.
However, all other terms set up by bfastpp
can also be omitted/added,
e.g., response ~ 1
(just a constant), response ~ season
(seasonal dummies for each period), etc. Further terms precomputed by
bfastpp
can be lag
(autoregressive terms of specified order),
slag
(seasonal autoregressive terms of specified order), xreg
(covariates, if data
has more than one column).
For determining the size of the stable history period, various approaches are available. First, the user can set a start date based on subjectmatter knowledge. Second, datadriven methods can be employed. By default, this is a reverseordered CUSUM test (ROC). Alternatively, breakpoints can be estimated (Bai and Perron method) and only the data after the last breakpoint are employed for the stable history. Finally, the user can also supply a function for his/her own datadriven method.
Value
bfastmonitor
returns an object of class
"bfastmonitor"
, i.e., a list with components as follows.
data 
original 
tspp 
preprocessed

model 
fitted

mefp 
fitted

history 
start and end time of history period, 
monitor 
start and end time of monitoring period, 
breakpoint 
breakpoint detected (if any). 
magnitude 
median of the difference between the data and the model prediction in the monitoring period. 
Author(s)
Achim Zeileis, Jan Verbesselt
References
Verbesselt J, Zeileis A, Herold M (2012). “Near realtime disturbance detection using satellite image time series.” Remote Sensing of Environment, 123, 98–108. ISSN 00344257, doi: 10.1016/j.rse.2012.02.022, https://doi.org/10.1016/j.rse.2012.02.022.
See Also
Examples
NDVIa < as.ts(zoo::zoo(som$NDVI.a, som$Time))
plot(NDVIa)
## apply the bfast monitor function on the data
## start of the monitoring period is c(2010, 13)
## and the ROC method is used as a method to automatically identify a stable history
mona < bfastmonitor(NDVIa, start = c(2010, 13))
mona
plot(mona)
## fitted seasontrend model in history period
summary(mona$model)
## OLSbased MOSUM monitoring process
plot(mona$mefp, functional = NULL)
## the pattern in the running mean of residuals
## this illustrates the empirical fluctuation process
## and the significance of the detected break.
NDVIb < as.ts(zoo(som$NDVI.b, som$Time))
plot(NDVIb)
monb < bfastmonitor(NDVIb, start = c(2010, 13))
monb
plot(monb)
summary(monb$model)
plot(monb$mefp, functional = NULL)
## set the stable history period manually and use a 4th order harmonic model
bfastmonitor(NDVIb, start = c(2010, 13),
history = c(2008, 7), order = 4, plot = TRUE)
## just use a 6th order harmonic model without trend
mon < bfastmonitor(NDVIb, formula = response ~ harmon,
start = c(2010, 13), order = 6, plot = TRUE)
summary(mon$model)
AIC(mon$model)
## use a custom number of seasonal dummies (11/yr) instead of harmonics
mon < bfastmonitor(NDVIb, formula = response ~ season,
start = c(2010, 13), sbins = 11, plot = TRUE)
summary(mon$model)
AIC(mon$model)
## Example for processing raster bricks (satellite image time series of 16day NDVI images)
f < system.file("extdata/modisraster.grd", package = "bfast")
library("raster")
modisbrick < raster::brick(f)
data < as.vector(modisbrick[1])
ndvi < bfastts(data, dates, type = c("16day"))
plot(ndvi/10000)
## derive median NDVI of a NDVI raster brick
medianNDVI < raster::calc(modisbrick, fun = function(x) median(x, na.rm = TRUE))
raster::plot(medianNDVI)
## helper function to be used with the calc() function
xbfastmonitor < function(x, timestamps = dates) {
ndvi < bfastts(x, timestamps, type = c("16day"))
ndvi < window(ndvi, end = c(2011, 14))/10000
## delete end of the time to obtain a dataset similar to RSE paper (Verbesselt et al.,2012)
bfm < bfastmonitor(data = ndvi, start = c(2010, 12), history = c("ROC"))
return(c(breakpoint = bfm$breakpoint, magnitude = bfm$magnitude))
}
## apply on one pixel for testing
ndvi < bfastts(as.numeric(modisbrick[1])/10000, dates, type = c("16day"))
plot(ndvi)
bfm < bfastmonitor(data = ndvi, start = c(2010, 12), history = c("ROC"))
bfm$magnitude
plot(bfm)
xbfastmonitor(modisbrick[1], dates) ## helper function applied on one pixel
## apply the bfastmonitor function onto a raster brick
timeofbreak < raster::calc(modisbrick, fun=xbfastmonitor)
raster::plot(timeofbreak) ## time of break and magnitude of change
raster::plot(timeofbreak,2) ## magnitude of change