Version:	1.7.0
Title:	Breaks for Additive Season and Trend
Description:	Decomposition of time series into trend, seasonal, and remainder components with methods for detecting and characterizing abrupt changes within the trend and seasonal components. 'BFAST' can be used to analyze different types of satellite image time series and can be applied to other disciplines dealing with seasonal or non-seasonal time series, such as hydrology, climatology, and econometrics. The algorithm can be extended to label detected changes with information on the parameters of the fitted piecewise linear models. 'BFAST' monitoring functionality is described in Verbesselt et al. (2010) <doi:10.1016/j.rse.2009.08.014>. 'BFAST monitor' provides functionality to detect disturbance in near real-time based on 'BFAST'- type models, and is described in Verbesselt et al. (2012) <doi:10.1016/j.rse.2012.02.022>. 'BFAST Lite' approach is a flexible approach that handles missing data without interpolation, and will be described in an upcoming paper. Furthermore, different models can now be used to fit the time series data and detect structural changes (breaks).
Depends:	R (≥ 3.0.0), strucchangeRcpp (≥ 1.5-4)
Imports:	graphics, stats, zoo, forecast, Rcpp (≥ 0.12.7), Rdpack (≥ 0.7)
Suggests:	MASS, sfsmisc, stlplus, terra
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
URL:	https://bfast2.github.io/
BugReports:	https://github.com/bfast2/bfast/issues
LazyLoad:	yes
LazyData:	yes
LinkingTo:	Rcpp
RoxygenNote:	7.3.2
RdMacros:	Rdpack
Encoding:	UTF-8
NeedsCompilation:	yes
Packaged:	2024-10-22 08:51:57 UTC; dainius
Author:	Jan Verbesselt [aut], Dainius Masiliūnas [aut, cre], Achim Zeileis [aut], Rob Hyndman [ctb], Marius Appel [aut], Martin Jung [ctb], Andrei Mîrț [ctb], Paulo Negri Bernardino [ctb], Dongdong Kong [ctb]
Maintainer:	Dainius Masiliūnas <pastas4@gmail.com>
Repository:	CRAN
Date/Publication:	2024-10-22 10:10:02 UTC

Breaks For Additive Season and Trend (BFAST)

Description

BFAST integrates the decomposition of time series into trend, seasonal, and remainder components with methods for detecting and characterizing abrupt changes within the trend and seasonal components. BFAST can be used to analyze different types of satellite image time series and can be applied to other disciplines dealing with seasonal or non-seasonal time series,such as hydrology, climatology, and econometrics. The algorithm can be extended to label detected changes with information on the parameters of the fitted piecewise linear models.

Additionally monitoring disturbances in BFAST-type models at the end of time series (i.e., in near real-time) is available: 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 season-trend model (with harmonic seasonal pattern) is used as a default in the regresssion modelling.

Details

The package contains:

• bfast(): Main function for iterative decomposition and break detection as described in Verbesselt et al (2010a,b).

• bfastlite(): lightweight and fast detection of all breaks in a time series using a single iteration with all components at once.

• bfastmonitor(): Monitoring approach for detecting disturbances in near real-time (see Verbesselt et al. 2012).

• bfastpp(): Data pre-processing for BFAST-type modeling.

• Functions for plotting and printing, see bfast().

• simts: Artificial example data set.

• harvest: NDVI time series of a P. radiata plantation that is harvested.

• som: NDVI time series of locations in the south of Somalia to illustrate the near real-time disturbance approach

Package options

bfast uses the following options to modify the default behaviour:

• bfast.prefer_matrix_methods: logical value defining whether methods should try to use the design matrix instead of the formula and a dataframe whenever possible. This can avoid expensive repeated calls of model.matrix and model.frame and make model fitting faster using lm.fit.

• bfast.use_bfastts_modifications: logical value defining whether a faster version of bfastts() should be used.

• strucchange.use_armadillo: logical value defining whether to use C++ optimised code paths in strucchangeRcpp.

By default, all three are enabled. See set_fallback_options() for a convenient interface for setting them all off for debugging purposes.

References

Verbesselt J, Zeileis A, Herold M (2012). “Near real-time disturbance detection using satellite image time series.” Remote Sensing of Environment, 123, 98–108. ISSN 0034-4257, doi:10.1016/j.rse.2012.02.022.

Verbesselt J, Hyndman R, Newnham G, Culvenor D (2010). “Detecting trend and seasonal changes in satellite image time series.” Remote Sensing of Environment, 114(1), 106–115. ISSN 0034-4257, doi:10.1016/j.rse.2009.08.014.

Verbesselt J, Hyndman R, Zeileis A, Culvenor D (2010). “Phenological change detection while accounting for abrupt and gradual trends in satellite image time series.” Remote Sensing of Environment, 114(12), 2970–2980. ISSN 0034-4257, doi:10.1016/j.rse.2010.08.003.

For all elements of a vector a, find the closest elements in a vector B and returns resulting indexes

Description

For all elements of a vector a, find the closest elements in a vector B and returns resulting indexes

Usage

.bfast_cpp_closestfrom(a, b, twosided)

Arguments

Value

integer vector of the same size as a with elements represnting indexes pointing to closest values in b

Break Detection in the Seasonal and Trend Component of a Univariate Time Series

Description

Iterative break detection in seasonal and trend component of a time series. Seasonal breaks is a function that combines the iterative decomposition of time series into trend, seasonal and remainder components with significant break detection in the decomposed components of the time series.

Usage

bfast(
  Yt,
  h = 0.15,
  season = c("dummy", "harmonic", "none"),
  max.iter = 10,
  breaks = NULL,
  hpc = "none",
  level = 0.05,
  decomp = c("stl", "stlplus"),
  type = "OLS-MOSUM",
  ...
)

Arguments

Details

The algorithm decomposes the input time series Yt into three components: trend, seasonality and remainder, using the function defined by the decomp parameter. Then each component is checked for at least one significant break using strucchangeRcpp::efp(), and if there is one, strucchangeRcpp::breakpoints() is run on the component. The result allows differentiating between breaks in trend and seasonality.

Value

An object of the class "bfast" is a list with the following elements:

Author(s)

Jan Verbesselt

References

Verbesselt J, Hyndman R, Newnham G, Culvenor D (2010). “Detecting trend and seasonal changes in satellite image time series.” Remote Sensing of Environment, 114(1), 106–115. ISSN 0034-4257, doi:10.1016/j.rse.2009.08.014.

Verbesselt J, Hyndman R, Zeileis A, Culvenor D (2010). “Phenological change detection while accounting for abrupt and gradual trends in satellite image time series.” Remote Sensing of Environment, 114(12), 2970–2980. ISSN 0034-4257, doi:10.1016/j.rse.2010.08.003.

See Also

plot.bfast for plotting of bfast() results.
breakpoints for more examples and background information about estimation of breakpoints in time series.

Examples

## Simulated Data
plot(simts) # stl object containing simulated NDVI time series
datats <- ts(rowSums(simts$time.series))
# sum of all the components (season,abrupt,remainder)
tsp(datats) <- tsp(simts$time.series) # assign correct time series attributes
plot(datats)

fit <- bfast(datats, h = 0.15, season = "dummy", max.iter = 1)
plot(fit, sim = simts)
fit
# prints out whether breakpoints are detected
# in the seasonal and trend component
 
## Real data
## The data should be a regular ts() object without NA's
## See Fig. 8 b in reference
plot(harvest, ylab = "NDVI")
# MODIS 16-day cleaned and interpolated NDVI time series

(rdist <- 10/length(harvest))
# ratio of distance between breaks (time steps) and length of the time series

fit <- bfast(harvest, h = rdist, season = "harmonic", max.iter = 1, breaks = 2)
plot(fit)
## plot anova and slope of the trend identified trend segments
plot(fit, ANOVA = TRUE)
## plot the trend component and identify the break with
## the largest magnitude of change
plot(fit, type = "trend", largest = TRUE)

## plot all the different available plots
plot(fit, type = "all")

## output
niter <- length(fit$output) # nr of iterations
out <- fit$output[[niter]]
# output of results of the final fitted seasonal and trend models and
## #nr of breakpoints in both.

## running bfast on yearly data
t <- ts(as.numeric(harvest), frequency = 1, start = 2006)
fit <- bfast(t, h = 0.23, season = "none", max.iter = 1)
plot(fit)
fit

## handling missing values with stlplus
(NDVIa <- as.ts(zoo::zoo(som$NDVI.a, som$Time)))
fit <- bfast(NDVIa, season = "harmonic", max.iter = 1, decomp = "stlplus")
plot(fit)
fit

Deprecated functions in package bfast.

Description

The functions listed below are deprecated and will be removed in the future. Please use the alternative functions as listed below. Help pages for deprecated functions are available via help("<function>-deprecated").

Usage

create16dayts(data, dates)

create16dayts

Use bfastts with the parameter type="16-day".

Checking for one major break in the time series

Description

A function to select a suitable model for the data by choosing either a model with 0 or with 1 breakpoint.

Usage

bfast01(
  data,
  formula = NULL,
  test = "OLS-MOSUM",
  level = 0.05,
  aggregate = all,
  trim = NULL,
  bandwidth = 0.15,
  functional = "max",
  order = 3,
  lag = NULL,
  slag = NULL,
  na.action = na.omit,
  reg = c("lm", "rlm"),
  stl = "none",
  sbins = 1
)

Arguments

Details

bfast01 tries to select a suitable model for the data by choosing either a model with 0 or with 1 breakpoint. It proceeds in the following steps:

1. The data is preprocessed with bfastpp using the arguments order/lag/slag/na.action/stl/sbins.

2. A linear model with the given formula is fitted. By default a suitable formula is guessed based on the preprocessing parameters.

3. The model with 1 breakpoint is estimated as well where the breakpoint is chosen to minimize the segmented residual sum of squares.

4. A sequence of tests for the null hypothesis of zero breaks is performed. Each test results in a decision for FALSE (no breaks) or TRUE (structural break(s)). The test decisions are then aggregated to a single decision (by default using all() but any() or some other function could also be used).

Available methods for the object returned include standard methods for linear models (coef, fitted, residuals, predict, AIC, BIC, logLik, deviance, nobs, model.matrix, model.frame), standard methods for breakpoints (breakpoints, breakdates), coercion to a zoo series with the decomposed components (as.zoo), and a plot method which plots such a zoo series along with the confidence interval (if the 1-break model is visualized). All methods take a 'breaks' argument which can either be 0 or 1. By default the value chosen based on the 'test' decisions is used.

Note that the different tests supported have power for different types of alternatives. Some tests (such as supLM/supF or BIC) assess changes in all coefficients of the model while residual-based tests (e.g., OLS-CUSUM or OLS-MOSUM) assess changes in the conditional mean. See Zeileis (2005) for a unifying view.

Value

bfast01 returns a list of class "bfast01" with the following elements:

Author(s)

Achim Zeileis, Jan Verbesselt

References

De Jong R, Verbesselt J, Zeileis A, Schaepman ME (2013). “Shifts in Global Vegetation Activity Trends.” Remote Sensing, 5(3), 1117–1133. ISSN 2072-4292, doi:10.3390/rs5031117.

Zeileis A (2005). “A Unified Approach to Structural Change Tests Based on ML Scores, F Statistics, and OLS Residuals.” Econometric Reviews, 24(4), 445–466. ISSN 0747-4938, doi:10.1080/07474930500406053.

See Also

bfastmonitor, breakpoints

Examples

library(zoo)
## define a regular time series
ndvi <- as.ts(zoo(som$NDVI.a, som$Time))

## fit variations
bf1 <- bfast01(ndvi)
bf2 <- bfast01(ndvi, test = c("BIC", "OLS-MOSUM", "supLM"), aggregate = any)
bf3 <- bfast01(ndvi, test = c("OLS-MOSUM", "supLM"), aggregate = any, bandwidth = 0.11) 

## inspect test decisions
bf1$test
bf1$breaks
bf2$test
bf2$breaks
bf3$test
bf3$breaks

## look at coefficients
coef(bf1)
coef(bf1, breaks = 0)
coef(bf1, breaks = 1) 

## zoo series with all components
plot(as.zoo(ndvi))
plot(as.zoo(bf1, breaks = 1))
plot(as.zoo(bf2))
plot(as.zoo(bf3))

## leveraged by plot method
plot(bf1, regular = TRUE)
plot(bf2)
plot(bf2, plot.type = "multiple",
     which = c("response", "trend", "season"), screens = c(1, 1, 2))
plot(bf3)

Change type analysis of the bfast01 function

Description

A function to determine the change type

Usage

bfast01classify(
  object,
  alpha = 0.05,
  pct_stable = NULL,
  typology = c("standard", "drylands")
)

Arguments

Details

bfast01classify

Value

bfast01classify returns a data.frame with the following elements:

and also significance and percentage of both segments before and after the potentially detected break: "p_segment1", "p_segment2", "pct_segment1", "pct_segment2".

Author(s)

Rogier de Jong, Jan Verbesselt

References

Bernardino PN, De Keersmaecker W, Fensholt R, Verbesselt J, Somers B, Horion S (2020). “Global-scale characterization of turning points in arid and semi-arid ecosystem functioning.” Global Ecology and Biogeography, 29(7), 1230–1245. doi:10.1111/geb.13099.

De Jong R, Verbesselt J, Zeileis A, Schaepman ME (2013). “Shifts in Global Vegetation Activity Trends.” Remote Sensing, 5(3), 1117–1133. ISSN 2072-4292, doi:10.3390/rs5031117.

See Also

bfast01

Examples

library(zoo)
## define a regular time series
ndvi <- as.ts(zoo(som$NDVI.a, som$Time))
## fit variations
bf1 <- bfast01(ndvi)
bfast01classify(bf1, pct_stable = 0.25)

Detect multiple breaks in a time series

Description

A combination of bfastpp and breakpoints to do light-weight detection of multiple breaks in a time series while also being able to deal with NA values by excluding them via bfastpp.

Usage

bfastlite(
  data,
  formula = response ~ trend + harmon,
  order = 3,
  breaks = "LWZ",
  lag = NULL,
  slag = NULL,
  na.action = na.omit,
  stl = c("none", "trend", "seasonal", "both"),
  decomp = c("stl", "stlplus"),
  sbins = 1,
  h = 0.15,
  level = 0,
  type = "OLS-MOSUM",
  ...
)

bfast0n(
  data,
  formula = response ~ trend + harmon,
  order = 3,
  breaks = "LWZ",
  lag = NULL,
  slag = NULL,
  na.action = na.omit,
  stl = c("none", "trend", "seasonal", "both"),
  decomp = c("stl", "stlplus"),
  sbins = 1,
  h = 0.15,
  level = 0,
  type = "OLS-MOSUM",
  ...
)

Arguments

Value

An object of class bfastlite, with three elements:

Author(s)

Dainius Masiliunas, Jan Verbesselt

References

Masiliūnas D, Tsendbazar N, Herold M, Verbesselt J (2021). “BFAST Lite: A Lightweight Break Detection Method for Time Series Analysis.” Remote Sensing, 13(16), 3308. doi:10.3390/rs13163308.

Examples

plot(simts) # stl object containing simulated NDVI time series
datats <- ts(rowSums(simts$time.series))
# sum of all the components (season,abrupt,remainder)
tsp(datats) <- tsp(simts$time.series) # assign correct time series attributes
plot(datats)

# Detect breaks
bp = bfastlite(datats)

# Default method of estimating breakpoints
bp[["breakpoints"]][["breakpoints"]]

# Plot
plot(bp)

# Custom method of estimating number of breaks (request 2 breaks)
strucchangeRcpp::breakpoints(bp[["breakpoints"]], breaks = 2)

# Plot including magnitude based on RMSD for the cos1 component of harmonics
plot(bp, magstat = "RMSD", magcomp = "harmoncos1", breaks = 2)

# Try with a structural change test
bp <- bfastlite(datats, level=0.05)
print(bp)
plot(bp)

# Details of the structural change test with the type RE
bfastlite(datats, level=0.05, type="RE")$sctest


## Run bfastlite() on a raster
f <- system.file("extdata/modisraster.tif", package="bfast")
modisbrick <- terra::rast(f)

# Run on pixel 10
data <- unlist(modisbrick[10])
ndvi <- bfastts(data, dates, type = c("16-day"))
bfl <- bfastlite(ndvi, breaks = "BIC")
# Get max magnitude by RMSD
max(magnitude(bfl[["breakpoints"]])$Mag[,"RMSD"])

# Wrapper function
bflSpatial <- function(pixels)
{
    ts <- bfastts(pixels, dates, type = c("16-day"))
    bfl <- bfastlite(ts, breaks="BIC")
    bp <- bfl[["breakpoints"]]
    # Return number of breakpoints and max breakpoint magnitude
    if (length(bp[["breakpoints"]]) == 1 && is.na(bp[["breakpoints"]]))
        return(c(0, 0))
    
    return(c(length(bp[["breakpoints"]]), max(magnitude(bp)$Mag[,"RMSD"])))
}

# Run function on each raster pixel
rastbfl <- terra::app(modisbrick, bflSpatial)
terra::plot(rastbfl)

Near Real-Time Disturbance Detection Based on BFAST-Type Models

Description

Monitoring disturbances in time series models (with trend/season/regressor terms) at the end of time series (i.e., in near real-time). 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 season-trend 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 = "OLS-MOSUM",
  h = 0.25,
  end = 10,
  level = c(0.05, 0.05),
  hpc = "none",
  verbose = FALSE,
  plot = FALSE,
  sbins = 1
)

Arguments

Details

bfastmonitor provides monitoring of disturbances (or structural changes) in near real-time 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 subject-matter knowledge. Second, data-driven methods can be employed. By default, this is a reverse-ordered 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 data-driven method.

Value

bfastmonitor returns an object of class "bfastmonitor", i.e., a list with components as follows.

Author(s)

Achim Zeileis, Jan Verbesselt

References

Verbesselt J, Zeileis A, Herold M (2012). “Near real-time disturbance detection using satellite image time series.” Remote Sensing of Environment, 123, 98–108. ISSN 0034-4257, doi:10.1016/j.rse.2012.02.022.

See Also

monitor, mefp, breakpoints

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 season-trend model in history period
summary(mona$model)
## OLS-based 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 16-day NDVI images)
f <- system.file("extdata/modisraster.tif", package="bfast")
modisbrick <- terra::rast(f)
data <- unlist(modisbrick[1])
ndvi <- bfastts(data, dates, type = c("16-day"))
plot(ndvi/10000)

## derive median NDVI of a NDVI raster brick
medianNDVI <- terra::app(modisbrick, fun = "median")
terra::plot(medianNDVI)

## helper function to be used with the app() function
xbfastmonitor <- function(x, timestamps = dates) {
	ndvi <- bfastts(x, timestamps, type = c("16-day"))
	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
bfm <- bfastmonitor(data = ndvi, start = c(2010, 12), history = c("ROC"))
bfm$magnitude
plot(bfm)
xbfastmonitor(data, dates) ## helper function applied on one pixel

## apply the bfastmonitor function onto a raster brick
timeofbreak <- terra::app(modisbrick, fun=xbfastmonitor)

terra::plot(timeofbreak) ## time of break and magnitude of change
terra::plot(timeofbreak,2) ## magnitude of change

Time Series Preprocessing for BFAST-Type Models

Description

Time series preprocessing for subsequent regression modeling. Based on a (seasonal) time series, a data frame with the response, seasonal terms, a trend term, (seasonal) autoregressive terms, and covariates is computed. This can subsequently be employed in regression models.

Usage

bfastpp(
  data,
  order = 3,
  lag = NULL,
  slag = NULL,
  na.action = na.omit,
  stl = c("none", "trend", "seasonal", "both"),
  decomp = c("stl", "stlplus"),
  sbins = 1
)

Arguments

Details

To facilitate (linear) regression models of time series data, bfastpp facilitates preprocessing and setting up regressor terms. It returns a data.frame containing the first column of the data as the response while further columns (if any) are used as covariates xreg. Additionally, a linear trend, seasonal dummies, harmonic seasonal terms, and (seasonal) autoregressive terms are provided.

Optionally, each column of data can be seasonally adjusted and/or trend-adjusted via STL (season-trend decomposition via LOESS smoothing) prior to preprocessing. The idea would be to capture season and/or trend nonparametrically prior to regression modelling.

Value

If no formula is provided, bfastpp returns a "data.frame" with the following variables (some of which may be matrices).

If a formula is given, bfastpp returns a list with components X, y, and t, where X is the design matrix of the model, y is the response vector, and t represents the time of observations. X will only contain variables that occur in the formula. Columns of X have names as decribed above.

Author(s)

Achim Zeileis

References

Verbesselt J, Zeileis A, Herold M (2012). “Near real-time disturbance detection using satellite image time series.” Remote Sensing of Environment, 123, 98–108. ISSN 0034-4257, doi:10.1016/j.rse.2012.02.022.

See Also

bfastmonitor

Examples

## set up time series
ndvi <- as.ts(zoo::zoo(cbind(a = som$NDVI.a, b = som$NDVI.b), som$Time))
ndvi <- window(ndvi, start = c(2006, 1), end = c(2009, 23))

## parametric season-trend model
d1 <- bfastpp(ndvi, order = 2)
d1lm <- lm(response ~ trend + harmon, data = d1)
summary(d1lm)
# plot visually (except season, as it's a factor)
plot(zoo::read.zoo(d1)[,-3],
  # Avoid clipping plots for pretty output
  ylim = list(c(min(d1[,2]), max(d1[,2])),
            c(min(d1[,3]), max(d1[,3])),
            c(-1, 1), c(-1, 1), c(-1, 1), c(-1, 1),
            c(min(d1[,6]), max(d1[,6]))
       ))

## autoregressive model (after nonparametric season-trend adjustment)
d2 <- bfastpp(ndvi, stl = "both", lag = 1:2)
d2lm <- lm(response ~ lag, data = d2)
summary(d2lm)

## use the lower level lm.fit function
d3 <- bfastpp(ndvi, stl = "both", lag = 1:2)
d3mm <- model.matrix(response ~ lag, d3)
d3lm <- lm.fit(d3mm, d3$response)
d3lm$coefficients

Create a regular time series object by combining data and date information

Description

Create a regular time series object by combining measurements (data) and time (dates) information.

Usage

bfastts(data, dates, type = c("irregular", "16-day", "10-day"))

Arguments

Details

bfastts create a regular time series

Value

bfastts returns an object of class "ts", i.e., a list with components as follows.

Author(s)

Achim Zeileis, Jan Verbesselt

See Also

monitor, mefp, breakpoints

Examples

# 16-day time series (i.e. MODIS)
timedf <- data.frame(y = som$NDVI.b, dates = dates[1:nrow(som)])
bfastts(timedf$y, timedf$dates, type = "16-day")

# Irregular
head(bfastts(timedf$y, timedf$dates, type = "irregular"), 50)


# Example of use with a raster
f <- system.file("extdata/modisraster.tif", package="bfast")
modisbrick <- terra::rast(f)
ndvi <- bfastts(unlist(modisbrick[1]), dates, type = c("16-day")) ## data of pixel 1
plot(ndvi/10000) 

# Time series of 4 pixels
modis_ts = t(modisbrick[1:4])
# Data with multiple columns, 2-4 are external regressors
ndvi <- bfastts(modis_ts, dates, type = c("16-day"))
plot(ndvi/10000)

A helper function to create time series

Description

A deprecated alias to bfastts. Please use bfastts(type="16-day") instead.

Usage

create16dayts(data, dates)

Arguments

Author(s)

Achim Zeileis, Jan Verbesselt

See Also

bfastmonitor

bfast-deprecated

A vector with date information (a Datum type) to be linked with each NDVI layer within the modis raster datacube (modisraster data set)

Description

dates is an object of class "Date" and contains the "Date" information to create a 16-day time series object.

Source

Verbesselt J, Zeileis A, Herold M (2012). “Near real-time disturbance detection using satellite image time series.” Remote Sensing of Environment, 123, 98–108. ISSN 0034-4257, doi:10.1016/j.rse.2012.02.022.

Examples

## see ?bfastmonitor for examples

16-day NDVI time series for a Pinus radiata plantation.

Description

A univariate time series object of class "ts". Frequency is set to 23 – the approximate number of observations per year.

Source

Verbesselt J, Hyndman R, Newnham G, Culvenor D (2010). “Detecting trend and seasonal changes in satellite image time series.” Remote Sensing of Environment, 114(1), 106–115. ISSN 0034-4257, doi:10.1016/j.rse.2009.08.014.

Examples

plot(harvest,ylab='NDVI')

A raster datacube of 16-day satellite image NDVI time series for a small subset in south eastern Somalia.

Description

A Cloud-Optimised GeoTIFF containing 16-day NDVI satellite images (MOD13C1 product).

Source

Verbesselt J, Zeileis A, Herold M (2012). “Near real-time disturbance detection using satellite image time series.” Remote Sensing of Environment, 123, 98–108. ISSN 0034-4257, doi:10.1016/j.rse.2012.02.022.

Examples

## see ?bfastmonitor

A random NDVI time series

Description

A univariate time series object of class "ts". Frequency is set to 24.

Examples

plot(ndvi)

Methods for objects of class "bfast".

Description

Plot methods for objects of class "bfast".

Usage

## S3 method for class 'bfast'
plot(
  x,
  type = c("components", "all", "data", "seasonal", "trend", "noise"),
  sim = NULL,
  largest = FALSE,
  main,
  ANOVA = FALSE,
  ...
)

Arguments

Details

This function creates various plots to demonstrate the results of a bfast decomposition. The type of plot shown depends on the value of type.

• components Shows the final estimated components with breakpoints.

• all Plots the estimated components and breakpoints from all iterations.

• data Just plots the original time series data.

• seasonal Shows the seasonal component including breakpoints.

• trend Shows the trend component including breakpoints.

• noise Plots the noise component along with its acf and pacf.

If sim is not NULL, the components used in simulation are also shown on each graph.

Value

No return value, called for side effects.

Author(s)

Jan Verbesselt, Rob Hyndman and Rogier De Jong

Examples

## See \code{\link[bfast]{bfast}} for examples.

Plot the time series and results of BFAST Lite

Description

The black line represents the original input data, the green line is the fitted model, the blue lines are the detected breaks, and the whiskers denote the magnitude (if magstat is specified).

Usage

## S3 method for class 'bfastlite'
plot(x, breaks = NULL, magstat = NULL, magcomp = "trend", ...)

Arguments

Value

Nothing, called for side effects.

Set package options with regard to computation times

Description

These functions set options of the bfast and strucchangeRcpp packages to enable faster computations. By default (set_default_options), these optimizations are enabled. Notice that only some functions of the bfast package make use of these options. set_fast_options is an alias for set_default_options.

Usage

set_default_options()

set_fast_options()

set_fallback_options()

Value

A list of modified options and their new values.

Examples

# run bfastmonitor with different options and compare computation times
library(zoo)
NDVIa <- as.ts(zoo(som$NDVI.a, som$Time))

set_default_options()
## Not run: 
system.time(replicate(100,  bfastmonitor(NDVIa, start = c(2010, 13))))

## End(Not run)

set_fallback_options()
## Not run: 
system.time(replicate(100,  bfastmonitor(NDVIa, start = c(2010, 13))))

## End(Not run)

Simulated seasonal 16-day NDVI time series

Description

simts is an object of class "stl" and consists of seasonal, trend (equal to 0) and noise components. The simulated noise is typical for remotely sensed satellite data.

Source

Verbesselt J, Hyndman R, Newnham G, Culvenor D (2010). “Detecting trend and seasonal changes in satellite image time series.” Remote Sensing of Environment, 114(1), 106–115. ISSN 0034-4257, doi:10.1016/j.rse.2009.08.014.

Examples

plot(simts)

Two 16-day NDVI time series from the south of Somalia

Description

som is a dataframe containing time and two NDVI time series to illlustrate how the monitoring approach works.

Source

Verbesselt J, Zeileis A, Herold M (2012). “Near real-time disturbance detection using satellite image time series.” Remote Sensing of Environment, 123, 98–108. ISSN 0034-4257, doi:10.1016/j.rse.2012.02.022.

Examples

## first define the data as a regular time series (i.e. ts object)
library(zoo)
NDVI <- as.ts(zoo(som$NDVI.b,som$Time))
plot(NDVI)