Movatterモバイル変換

Version:

1.7.1

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:

2025-08-07 09:14:37 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:

2025-08-18 10:20:14 UTC

Breaks For Additive Season and Trend (BFAST)

Description

BFAST integrates the decomposition of time series into trend, seasonal, andremainder components with methods for detecting and characterizing abruptchanges within the trend and seasonal components. BFAST can be used toanalyze different types of satellite image time series and can be applied toother disciplines dealing with seasonal or non-seasonal time series,such ashydrology, climatology, and econometrics. The algorithm can be extended tolabel detected changes with information on the parameters of the fittedpiecewise linear models.

Additionally monitoring disturbances in BFAST-type models at the end of timeseries (i.e., in near real-time) is available: Based on a model for stablehistorical behaviour abnormal changes within newly acquired data can bedetected. Different models are available for modeling the stable historicalbehavior. A season-trend model (with harmonic seasonal pattern) is used as adefault in the regresssion modelling.

Details

The package contains:

bfast(): Main function for iterative decomposition and break detection as described inVerbesselt et al (2010a,b).
bfastlite(): lightweight and fast detection of all breaks in a time seriesusing a single iteration with all components at once.
bfastmonitor():Monitoring approach for detecting disturbances in near real-time (seeVerbesselt et al. 2012).
bfastpp(): Data pre-processing for BFAST-type modeling.
Functions for plotting and printing, seebfast().
simts: Artificial example data set.
harvest: NDVI time series of a P. radiata plantationthat is harvested.
som: NDVI time series oflocations in the south of Somalia to illustrate the near real-timedisturbance approach

Package options

bfast uses the following options to modify thedefault behaviour:

bfast.prefer_matrix_methods:logical value defining whether methods should try touse the design matrix instead of the formula and a dataframe wheneverpossible. This can avoid expensive repeated calls ofmodel.matrix andmodel.frame and make model fitting faster usinglm.fit.
bfast.use_bfastts_modifications:logical value defining whether a faster version ofbfastts() should be used.
strucchange.use_armadillo:logical value defining whether to use C++ optimised code paths in strucchangeRcpp.

By default, all three are enabled.Seeset_fallback_options() for a convenient interface for setting them all offfor 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

a

numeric vector, ordered

b

numeric vector, ordered

twosided

logical value, if false, indexes will always point to elements in b that are less than or equal to elements in a but not greater than.

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 TimeSeries

Description

Iterative break detection in seasonal and trend component of a time series.Seasonal breaks is a function that combines the iterative decomposition oftime series into trend, seasonal and remainder components with significantbreak 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

Yt

univariate time series to be analyzed. This should be an object ofclass "ts" with a frequency greater than one.

h

minimal segment size between potentially detected breaks in thetrend model given as fraction relative to the sample size (i.e. the minimalnumber of observations in each segment divided by the total length of thetimeseries).

season

the seasonal model used to fit the seasonal component anddetect seasonal breaks (i.e. significant phenological change). There arethree options: "dummy", "harmonic", or "none" where "dummy" is the modelproposed in the first Remote Sensing of Environment paper and "harmonic" isthe model used in the second Remote Sensing of Environment paper (See paperfor more details) and where "none" indicates that no seasonal model will befitted (i.e. St = 0 ). If there is no seasonal cycle (e.g. frequency of thetime series is 1) "none" can be selected to avoid fitting a seasonal model.

max.iter

maximum amount of iterations allowed for estimation ofbreakpoints in seasonal and trend component.

breaks

either an integer specifying the number of breaks to compute,or a string specifying a statistic with which to computethe optimal number of breakpoints (seestrucchangeRcpp::breakpoints() formore information). If one value is given, it is used for both the trend andseasonal components, and if two are given, the first one is considered to befor the trend and the second for the seasonal component.

hpc

A character specifying the high performance computing support.Default is "none", can be set to "foreach". Install the "foreach" packagefor hpc support.

level

numeric; threshold value for thesctest.efptest; if a length 2 vector is passed, the first value is used for the trend,the second for the seasonality

decomp

"stlplus" or "stl": the function to use for decomposition.stl can handle sparse time series (1 < frequency < 4),stlpluscan handleNA values in the time series.

type

character, indicating the type argument toefp

...

additional arguments passed tostlplus::stlplus(), ifits usage has been enabled, otherwise ignored.

Details

The algorithm decomposes the input time seriesYt into three components:trend, seasonality and remainder, using the function defined by thedecompparameter. Then each component is checked for at least one significantbreak usingstrucchangeRcpp::efp(), and if there is one,strucchangeRcpp::breakpoints()is run on the component. The result allows differentiating between breaks intrend and seasonality.

Value

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

Yt

equals the Yt used as input.

output

is a listwith the following elements (for each iteration):

`Tt`	the fitted trend component
`St`	the fitted seasonalcomponent
`Nt`	the noise or remainder component
`Vt`	equals the deseasonalized data`Yt - St` for each iteration
`bp.Vt`	output of the`breakpoints`function for the trend model. Note that the output breakpoints are indexnumbers of`na.omit(as.numeric(Vt))`.
`ci.Vt`	output of the`breakpoints` confint function for the trendmodel
`Wt`	equals the detrended data`Yt - Tt` for eachiteration
`bp.Wt`	output of the`breakpoints` function for the seasonal model.Note that the output breakpoints are index numbers of`na.omit(as.numeric(Wt))`.
`ci.Wt`	output of the`breakpoints`confint function for the seasonal model

nobp

is a list with thefollowing elements:

`nobp.Vt`	logical, TRUE if thereare no breakpoints detected
`nobp.Wt`	logical, TRUE if thereare no breakpoints detected

Magnitude

magnitude of the biggestchange detected in the trend component

Time

timing of the biggestchange detected in the trend component

Author(s)

Jan Verbesselt

References

Examples

## Simulated Dataplot(simts) # stl object containing simulated NDVI time seriesdatats <- ts(rowSums(simts$time.series))# sum of all the components (season,abrupt,remainder)tsp(datats) <- tsp(simts$time.series) # assign correct time series attributesplot(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 referenceplot(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 seriesfit <- bfast(harvest, h = rdist, season = "harmonic", max.iter = 1, breaks = 2)plot(fit)## plot anova and slope of the trend identified trend segmentsplot(fit, ANOVA = TRUE)## plot the trend component and identify the break with## the largest magnitude of changeplot(fit, type = "trend", largest = TRUE)## plot all the different available plotsplot(fit, type = "all")## outputniter <- length(fit$output) # nr of iterationsout <- fit$output[[niter]]# output of results of the final fitted seasonal and trend models and## #nr of breakpoints in both.## running bfast on yearly datat <- 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 packagebfast.

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 viahelp("<function>-deprecated").

Usage

create16dayts(data, dates)

`create16dayts`

Usebfastts with the parametertype="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 amodel 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

data

A time series of classts, or anotherobject that can be coerced to such. The time series is processed bybfastpp. A time series of classtscan be prepared by a convenience functionbfastts incase of daily, 10 or 16-daily time series.

formula

formula for the regression model. The default isintelligently guessed based on the arguments order/lag/slag i.e.response ~ trend + harmon, i.e., a linear trend and a harmonic seasoncomponent. Other specifications are possible using all terms set up bybfastpp, i.e.,season (seasonal pattern withdummy variables),lag (autoregressive terms),slag (seasonalautoregressiv terms), orxreg (further covariates). Seebfastpp for details.

test

character specifying the type of test(s) performed. Can be oneor more of BIC, supLM, supF, OLS-MOSUM, ..., or any other test supported bysctest.formula

level

numeric. Significance for thesctest.formula performed.

aggregate

function that aggregates a logical vector to a singlevalue. This is used for aggregating the individual test decisions fromtest to a single one.

trim

numeric. The mimimal segment size passed to thefromargument of theFstats function.

bandwidth

numeric scalar from interval (0,1), functional. Thebandwidth argument is passed to theh argument of thesctest.formula.

functional

arguments passed on tosctest.formula

order

numeric. Order of the harmonic term, defaulting to3.

lag

numeric. Order of the autoregressive term, by default omitted.

slag

numeric. Order of the seasonal autoregressive term, by defaultomitted.

na.action

arguments passed on tobfastpp

reg

whether to use OLS regressionlm() orrobust regressionMASS::rlm().

stl

argument passed on tobfastpp

sbins

argument passed on tobfastpp

Details

bfast01 tries to select a suitable model for the data by choosingeither a model with 0 or with 1 breakpoint. It proceeds in the followingsteps:

The data is preprocessed with bfastpp using the argumentsorder/lag/slag/na.action/stl/sbins.
A linear model with the given formula is fitted. By default a suitableformula is guessed based on the preprocessing parameters.
The model with 1 breakpoint is estimated as well where the breakpoint ischosen to minimize the segmented residual sum of squares.
A sequence of tests for the null hypothesis of zero breaks is performed. Eachtest results in a decision for FALSE (no breaks) or TRUE (structuralbreak(s)). The test decisions are then aggregated to a single decision (bydefault using all() but any() or some other function could also be used).

Available methods for the object returned include standard methods forlinear 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 decomposedcomponents (as.zoo), and a plot method which plots such a zoo series alongwith the confidence interval (if the 1-break model is visualized). Allmethods take a 'breaks' argument which can either be 0 or 1. By default thevalue chosen based on the 'test' decisions is used.

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

Value

bfast01 returns a list of class"bfast01" with thefollowing elements:

call

the original function call.

data

thedata preprocessed by"bfastpp".

formula

the model formulae.

breaks

the number of breaks chosen based on thetest decision(either 0 or 1).

test

the individual test decisions.

breakpoints

the optimal breakpoint for the model with 1 break.

model

A list of two 'lm' objects with no and one breaks,respectively.

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.

Examples

library(zoo)## define a regular time seriesndvi <- as.ts(zoo(som$NDVI.a, som$Time))## fit variationsbf1 <- 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 decisionsbf1$testbf1$breaksbf2$testbf2$breaksbf3$testbf3$breaks## look at coefficientscoef(bf1)coef(bf1, breaks = 0)coef(bf1, breaks = 1) ## zoo series with all componentsplot(as.zoo(ndvi))plot(as.zoo(bf1, breaks = 1))plot(as.zoo(bf2))plot(as.zoo(bf3))## leveraged by plot methodplot(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

object

bfast01 object, i.e. the output of thebfast01 function.

alpha

threshold for significance tests, default 0.05

pct_stable

threshold for segment stability, unit: percent change perunit time (0-100), default NULL

typology

classification legend to use:standard refers to theoriginal legend as used in De Jong et al. (2013),drylands refers to the legend used in Bernardino et al. (2020).

Details

bfast01classify

Value

bfast01classify returns a data.frame with the followingelements:

flag_type

Type of shift: (1) monotonic increase, (2)monotonic decrease, (3) monotonic increase (with positive break), (4)monotonic decrease (with negative break), (5) interruption: increase withnegative break, (6) interruption: decrease with positive break, (7)reversal: increase to decrease, (8) reversal: decrease to increase

flag_significance

SIGNIFICANCE FLAG: (0) both segments significant(or no break and significant), (1) only first segment significant, (2) only2nd segment significant, (3) both segments insignificant (or no break andnot significant)

flag_pct_stable

STABILITY FLAG: (0) change in bothsegments is substantial (or no break and substantial), (1) only firstsegment substantial, (2) only 2nd segment substantial (3) both segments arestable (or no break and stable)

and also significance and percentage ofboth 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.

Examples

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

Detect multiple breaks in a time series

Description

A combination ofbfastpp andbreakpointsto do light-weight detection of multiple breaks in a time serieswhile also being able to deal with NA values by excluding themviabfastpp.

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

data

A time series of classts, or anotherobject that can be coerced to such. For seasonal components, a frequencygreater than 1 is required.

formula

a symbolic description for the model in which breakpointswill be estimated.

order

numeric. Order of the harmonic term, defaulting to3.

breaks

either a positive integer specifying the maximal number of breaks to be calculated,or a string specifying the information criterion to use to automatically determinethe optimal number of breaks (see alsologLik)."all" means the maximal number allowed byh is used.NULL is treated as the default of thebreakpoints function (i.e. BIC).

lag

numeric. Orders of the autoregressive term, by default omitted.

slag

numeric. Orders of the seasonal autoregressive term, by defaultomitted.

na.action

function for handlingNAs in the data (after allother preprocessing).

stl

character. Prior to all other preprocessing, STL (season-trenddecomposition via LOESS smoothing) can be employed for trend-adjustmentand/or season-adjustment. The"trend" or"seasonal" componentor both fromstl are removed from each column indata. By default ("none"), no STL adjustment is used.

decomp

"stlplus" or "stl": use the NA-tolerant decomposition packageor the reference package (which can make use of time series with 2-3observations per year)

sbins

numeric. Controls the number of seasonal dummies. If integer > 1,sets the number of seasonal dummies to use per year.If <= 1, treated as a multiplier to the number of observations per year, i.e.ndummies = nobs/year * sbins.

h

minimal segment size either given as fraction relative to thesample size or as an integer giving the minimal number of observationsin each segment.

level

numeric; threshold value for thesctest.efptest for at least one break in a time series, to save processing time forno-change areas. The test is skipped if set to <= 0.

type

character, indicating the type argument toefp.

...

Additional arguments tobreakpoints.

Value

An object of classbfastlite, with three elements:

breakpoints

output frombreakpoints,containing information about the estimated breakpoints.

data_pp

preprocessed data as output bybfastpp.

sctest

output fromsctest.efp,contatining information about the likelihood of the time serieshaving at least one break.

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 seriesdatats <- ts(rowSums(simts$time.series))# sum of all the components (season,abrupt,remainder)tsp(datats) <- tsp(simts$time.series) # assign correct time series attributesplot(datats)# Detect breaksbp = bfastlite(datats)# Default method of estimating breakpointsbp[["breakpoints"]][["breakpoints"]]# Plotplot(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 harmonicsplot(bp, magstat = "RMSD", magcomp = "harmoncos1", breaks = 2)# Try with a structural change testbp <- bfastlite(datats, level=0.05)print(bp)plot(bp)# Details of the structural change test with the type REbfastlite(datats, level=0.05, type="RE")$sctest## Run bfastlite() on a rasterf <- system.file("extdata/modisraster.tif", package="bfast")modisbrick <- terra::rast(f)# Run on pixel 10data <- unlist(modisbrick[10])ndvi <- bfastts(data, dates, type = c("16-day"))bfl <- bfastlite(ndvi, breaks = "BIC")# Get max magnitude by RMSDmax(magnitude(bfl[["breakpoints"]])$Mag[,"RMSD"])# Wrapper functionbflSpatial <- 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 pixelrastbfl <- 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/regressorterms) at the end of time series (i.e., in near real-time). Based on a modelfor stable historical behaviour abnormal changes within newly acquired datacan be detected. Different models are available for modeling the stablehistorical 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

data

A time series of classts, or anotherobject that can be coerced to such. For seasonal components, a frequencygreater than 1 is required.

start

numeric. The starting date of the monitoring period. Can eitherbe given as a float (e.g.,2000.5) or a vector giving period/cycle(e.g.,c(2000, 7)).

formula

formula for the regression model. The default isresponse ~ trend + harmon, i.e., a linear trend and a harmonic seasoncomponent. Other specifications are possible using all terms set up bybfastpp, i.e.,season (seasonal pattern withdummy variables),lag (autoregressive terms),slag (seasonalautoregressive terms), orxreg (further covariates). Seebfastpp for details.

order

numeric. Order of the harmonic term, defaulting to3.

lag

numeric. Order of the autoregressive term, by default omitted.

slag

numeric. Order of the seasonal autoregressive term, by defaultomitted.

history

specification of the start of the stable history period. Caneither be a character, numeric, or a function. If character, then selectionis possible between reverse-ordered CUSUM ("ROC", default), Bai andPerron breakpoint estimation ("BP"), or all available observations("all"). If numeric, the start date can be specified in the same formasstart. If a function is supplied it is called ashistory(formula, data) to compute a numeric start date.

type

character specifying the type of monitoring process. By default,a MOSUM process based on OLS residuals is employed. Seemefp for alternatives.

h

numeric scalar from interval (0,1) specifying the bandwidthrelative to the sample size in MOSUM/ME monitoring processes.

end

numeric. Maximum time (relative to the history period) that willbe monitored (in MOSUM/ME processes). Default is 10 times the historyperiod.

level

numeric vector. Significance levels of the monitoring and ROC (ifselected) procedure, i.e., probability of type I error.

hpc

character specifying the high performance computing support.Default is"none", can be set to"foreach". Seebreakpoints for more details.

verbose

logical. Should information about the monitoring be printedduring computation?

plot

logical. Should the result be plotted?

sbins

numeric. Number of seasonal dummies, passed tobfastpp.

Details

bfastmonitor provides monitoring of disturbances (or structuralchanges) in near real-time based on a wide class of time series regressionmodels with optional season/trend/autoregressive/covariate terms. SeeVerbesselt at al. (2011) for details.

Based on a given time series (typically, but not necessarily, with frequencygreater than 1), the data is first preprocessed for regression modeling.Trend/season/autoregressive/covariate terms are (optionally) computed usingbfastpp. Second, the data is split into a history andmonitoring period (starting withstart). Third, a subset of thehistory period is determined which is considered to be stable (see alsobelow). Fourth, a regression model is fitted to the preprocessed data inthe stable history period. Fifth, a monitoring procedure is used todetermine whether the observations in the monitoring period conform withthis stable regression model or whether a change is detected.

The regression model can be specified by the user. The default is to use alinear trend and a harmonic season:response ~ trend + harmon.However, all other terms set up bybfastpp can also be omitted/added,e.g.,response ~ 1 (just a constant),response ~ season(seasonal dummies for each period), etc. Further terms precomputed bybfastpp can belag (autoregressive terms of specified order),slag (seasonal autoregressive terms of specified order),xreg(covariates, ifdata has more than one column).

For determining the size of the stable history period, various approachesare available. First, the user can set a start date based on subject-matterknowledge. Second, data-driven methods can be employed. By default, this isa reverse-ordered CUSUM test (ROC). Alternatively, breakpoints can beestimated (Bai and Perron method) and only the data after the lastbreakpoint are employed for the stable history. Finally, the user can alsosupply 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.

data

original"ts" time series,

tspp

preprocessed"data.frame" for regression modeling,

model

fitted"lm" model for the stable history period,

mefp

fitted"mefp" process for the monitoring period,

history

start andend time of history period,

monitor

start and end time of monitoringperiod,

breakpoint

breakpoint detected (if any).

magnitude

median of the difference between the data and the modelprediction in the monitoring period.

Author(s)

Achim Zeileis, Jan Verbesselt

References

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 historymona <- bfastmonitor(NDVIa, start = c(2010, 13))monaplot(mona)## fitted season-trend model in history periodsummary(mona$model)## OLS-based MOSUM monitoring processplot(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))monbplot(monb)summary(monb$model)plot(monb$mefp, functional = NULL)## set the stable history period manually and use a 4th order harmonic modelbfastmonitor(NDVIb, start = c(2010, 13),  history = c(2008, 7), order = 4, plot = TRUE)## just use a 6th order harmonic model without trendmon <- 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 harmonicsmon <- 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 brickmedianNDVI <- terra::app(modisbrick, fun = "median")terra::plot(medianNDVI)## helper function to be used with the app() functionxbfastmonitor <- 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 testingbfm <- bfastmonitor(data = ndvi, start = c(2010, 12), history = c("ROC"))bfm$magnitudeplot(bfm)xbfastmonitor(data, dates) ## helper function applied on one pixel## apply the bfastmonitor function onto a raster bricktimeofbreak <- terra::app(modisbrick, fun=xbfastmonitor)terra::plot(timeofbreak) ## time of break and magnitude of changeterra::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, atrend 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

data

A time series of classts, or anotherobject that can be coerced to such. For seasonal components, a frequencygreater than 1 is required.

order

numeric. Order of the harmonic term, defaulting to3.

lag

numeric. Orders of the autoregressive term, by default omitted.

slag

numeric. Orders of the seasonal autoregressive term, by defaultomitted.

na.action

function for handlingNAs in the data (after allother preprocessing).

stl

decomp

"stlplus" or "stl": use the NA-tolerant decomposition packageor the reference package (which can make use of time series with 2-3observations per year)

sbins

Details

To facilitate (linear) regression models of time series data,bfastppfacilitates preprocessing and setting up regressor terms. It returns adata.frame containing the first column of thedata as theresponse while further columns (if any) are used as covariatesxreg. Additionally, a linear trend, seasonal dummies, harmonicseasonal terms, and (seasonal) autoregressive terms are provided.

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

Value

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

time

numeric vector of time stamps,

response

response vector (first column ofdata),

trend

linear time trend (running from 1 to number of observations),

season

factor indicating season period,

harmon

harmonicseasonal terms (of specifiedorder),

lag

autoregressive terms(or orderslag, if any),

slag

seasonal autoregressive terms(or ordersslag, if any),

xreg

covariate regressor (allcolumns ofdata except the first, if any).

If a formula is given,bfastpp returns alist with componentsX,y, andt, whereX is the design matrix of themodel,y is the response vector, andt represents the time ofobservations.X will only contain variables that occur in theformula. Columns ofX have names as decribed above.

Author(s)

Achim Zeileis

References

Examples

## set up time seriesndvi <- 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 modeld1 <- 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 functiond3 <- 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) andtime (dates) information.

Usage

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

Arguments

data

A data vector or matrix where columns represent variables

dates

Optional input of dates for each measurement in the 'data'variable. In case the data is a irregular time series, a vector with 'dates'for each measurement can be supplied using this 'dates' variable. Theirregular data will be linked with the dates vector to create daily regulartime series with a frequency = 365. Extra days in leap years might causeproblems. Please be careful using this option as it is experimental.Feedback is welcome.

type

("irregular") indicates that the data is collected atirregular dates and as such will be converted to a daily time series.("16-day") indicates that data is collected at a regular time interval(every 16-days e.g. like the MODIS 16-day data products). ("10-day")indicates that data is collected at a 10-day time interval of the SPOT VEGETATION(S10) product. Warning: Only use this function for the SPOT VEGETATION S10 time series,as for other 10-day time series a different approach might be required.

Details

bfastts create a regular time series

Value

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

zz

a regular"ts" time serieswith a frequency equal to 365 or 23 i.e. 16-day time series.

Author(s)

Achim Zeileis, Jan Verbesselt

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")# Irregularhead(bfastts(timedf$y, timedf$dates, type = "irregular"), 50)# Example of use with a rasterf <- system.file("extdata/modisraster.tif", package="bfast")modisbrick <- terra::rast(f)ndvi <- bfastts(unlist(modisbrick[1]), dates, type = c("16-day")) ## data of pixel 1plot(ndvi/10000) # Time series of 4 pixelsmodis_ts = t(modisbrick[1:4])# Data with multiple columns, 2-4 are external regressorsndvi <- 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 usebfastts(type="16-day") instead.

Usage

create16dayts(data, dates)

Arguments

data

Passed to bfastts.

dates

Passed to bfastts.

Author(s)

Achim Zeileis, Jan Verbesselt

A vector with date information (a Datum type) to be linked with each NDVIlayer 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

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 – theapproximate number of observations per year.

Source

Examples

plot(harvest,ylab='NDVI')

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

Description

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

Source

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

x

bfast object

type

Indicates the type of plot. See details.

sim

Optionalstl object containing the originalcomponents used when simulatingx.

largest

If TRUE, show the largest jump in the trend component.

main

an overall title for the plot.

ANOVA

if TRUE Derive Slope and Significance values for eachidentified trend segment

...

further arguments passed to theplotfunction.

Details

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

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.

Ifsim is notNULL, the components used in simulation are alsoshown 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, andthe whiskers denote the magnitude (ifmagstat is specified).

Usage

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

Arguments

x

bfastlite object frombfastlite()

breaks

number of breaks or optimal break selection method, seestrucchangeRcpp::breakpoints()

magstat

name of the magnitude column to plot (e.g.RMSD,MAD,diff), see theMag component ofstrucchangeRcpp::magnitude.breakpointsfull()

magcomp

name of the component (i.e. column inx$data_pp) to plot magnitudes of

...

other parameters to pass toplot()

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 enablefaster computations. By default (set_default_options), these optimizations areenabled. Notice that only some functions of thebfastpackage make use of these options.set_fast_options is an alias forset_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 timeslibrary(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 forremotely sensed satellite data.

Source

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 toilllustrate how the monitoring approach works.

Source

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)

Movatterモバイル変換

Breaks For Additive Season and Trend (BFAST)

Description

Details

Package options

References

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

Description

Usage

Arguments

Value

Break Detection in the Seasonal and Trend Component of a Univariate TimeSeries

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

Deprecated functions in packagebfast.

Description

Usage

create16dayts

Checking for one major break in the time series

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

Change type analysis of the bfast01 function

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

Detect multiple breaks in a time series

Description

Usage

Arguments

Value

Author(s)

References

Examples

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

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

Time Series Preprocessing for BFAST-Type Models

Description

Usage

Arguments

Details

Value

Author(s)

References

See Also

Examples

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

Description

Usage

Arguments

Details

Value

Author(s)

`create16dayts`