Breaks for Additive Season and Trend (BFAST) Walkthrough

Overview

Breaks for Additive Season and Trend (BFAST) is a parametric time-series analysis method that applies an additive decomposition model to separate each pixel’s time series into three components:

Trend
Seasonal
Remainder (noise or residuals)

The primary objectives of BFAST are:

1. Detect multiple abrupt changes in the trend and seasonal components.

2. Characterize gradual and abrupt ecosystem changes by estimating the timing, magnitude, and direction of change.

Statistical significance of trends is assessed using a Student’s t-test to determine whether the slope significantly differs from zero. Based on slope sign and significance, vegetation change is classified as:

Positive
Negative
Stable

Data Acquisition

Data were acquired from the NASA AppEEARS website.

Parameters Used

Boundary: Ukraine shapefile
Start Date: 01-01-2001
End Date: 12-31-2023
Dataset: MOD13Q1
Resolution: 250 m
Temporal Frequency: 16-day
Product: EVI
Format: GeoTIFF
Projection: GCS WGS84

Rescaling EVI Values

MODIS EVI values range from -2000 to 10,000.

To convert to standard EVI values (-1 to 1), multiply by 0.0001.

EVI_rescaled <- EVI * 0.0001

Fishnet Creation

Fishnet tiles were created across Ukraine to reduce computational demand.

Workflow

1. Convert MODIS EVI imagery to Integer data type.

2. Use Raster to Points to extract center points.

3. Use Extract Multi Values to Points to add EVI values.

4. Export attribute table to .csv.

BFAST Analysis

Step 1: Handle Missing Values

BFAST will not run if NA values are present.

Ensure all missing values are filled or removed before analysis.

Step 2: Parallel BFAST Processing

Parallel processing was performed to reduce computation time. Each tile generates multiple output files containing trend and seasonal break information.

library("zoo")
library("sandwich")
library("MASS")
library("quadprog")
library("tseries")
library("strucchange")
library("fracdiff")
library("forecast")
library("iterators")
library("codetools")
library("foreach")
library("bfast")
#library("mailR")
library("parallel")

bfastFunction <- function(count) { #Function that gets run through the parallel process
   pid=Sys.getpid()          # PID for the forked process
   path=paste0(dir,"/out-tile",tileName,"-pid",pid,"/")
   system( paste0( "mkdir -p ", path) )
    poly_id<-tpdata[count,1]  #highlighted number identifies field that will used to ID pixel
    ndvi<-tpdata[count,2:vmax[2]]  #highlighted number identifies first column with NDVI data
    #plot(ndvi)
    tsdata<-ts(ndvi,frequency=23,start=c(2001,1))
    dim(tsdata)<-NULL
    rdist<-23/length(tsdata)
    fits<-bfast(tsdata,h=rdist,season="harmonic",max.iter=1)
    #plot(fits)
    fits2<-fits$Time
    ts_trend_break_time<-t(fits2[1])
    fits3<-fits$Magnitude
    ts_trend_break_magnitude<-t(fits3[1])
    fits4<-fits$output
    fits4a<-fits4[[1]]$Vt.bp
    fits4adata<-as.matrix(fits4a)
    fits4amax<-dim(fits4adata)
    ts_trend_nbbreak<-t(fits4amax[1])
    results1<-ts_trend_break_time
    aLine<-t(c(poly_id,results1))
    write.table(aLine, file=paste0(path,tileName, "_trend_breaks_time.txt"), append=TRUE,quote=FALSE,sep=",", eol="\n",na="NA", dec=".",row.names=FALSE,col.names=FALSE,qmethod=c("escape","double"))
    results2<-ts_trend_break_magnitude
    aLine<-t(c(poly_id,results2))
    write.table(aLine,file=paste0(path,tileName, "_trend_breaks_magnitude.txt"), append=TRUE,quote=FALSE, sep=",",eol="\n", na="NA", dec=".",row.names=FALSE,col.names=FALSE,qmethod=c("escape","double"))
    results3<-ts_trend_nbbreak
    aLine<-t(c(poly_id,results3))
    write.table(aLine,file=paste0(path,tileName, "_trend_nbbreaks.txt"),append=TRUE,quote=FALSE,sep=",", eol="\n",na="NA", de=".", ,row.names=FALSE,col.names=FALSE,qmethod=c("escape","double"))
    fits4b<-fits4[[1]]$Tt
    results4<-fits4b
    aLine<-t(c(poly_id,results4))
    write.table(aLine,file=paste0(path,tileName, "_trend_bfast.txt"),append=TRUE,quote=FALSE, sep=",",eol="\n",na="NA",dec=".", row.names=FALSE,col.names=FALSE,qmethod=c("escape","double"))
    fits4c<-fits4[[1]]$Wt.bp
    fits4cdata<-as.matrix(fits4c)
    fits4cmax<-dim(fits4cdata)
    ts_season_nbbreak<-t(fits4cmax[1])
    results5<-ts_season_nbbreak
    aLine<-t(c(poly_id,results5))
    write.table(aLine,file=paste0(path,tileName, "_season_nbbreaks.txt"),append=TRUE,quote=FALSE, sep=",",eol="\n", na="NA", de=".", row.names=FALSE,col.names=FALSE,qmethod=c("escape","double"))
    ts_season_breaks_time<-t(fits4cdata)
    results6<- ts_season_breaks_time
    aLine<- t(c(poly_id,results6))
    write.table(aLine,file=paste0(path,tileName, "_season_breaks_time.txt"),append=TRUE,quote=FALSE,sep=",", eol="\n",na="NA", de=".",row.names=FALSE,col.names=FALSE,qmethod=c("escape","double"))
    fits4d<-fits4[[1]]$St
    results7<-fits4d
    aLine<-t(c(poly_id,results7))
    write.table(aLine,file=paste0(path,tileName, "_season_bfast.txt"),append=TRUE,quote=FALSE,sep=",",eol="\n", na="NA",dec=".", row.names=FALSE,col.names=FALSE,qmethod=c("escape","double"))
}

args <- commandArgs(trailingOnly = TRUE)
no_cores <- strtoi(args[1]) #Get the number of cores as the first command line argument
no_cores
filename <- args[2] #Get the file to run as the second command line argument
filename
tileName <- strsplit(filename, "_")[[1]][2] #Get the tile number by splitting up the filename string

inidata<-read.table(filename, header=TRUE, sep = ",", dec = ".") #use only with small files; modify if no labels are in the input 
mdata<-as.matrix(inidata)
tpdata<-mdata
vmax<-dim(mdata)
vmax[1]
vmax[2]

dir="/tmp"
#dir="Output"

arr <- seq(1, vmax[1]) #List of numbers from 1 to vmax[1], incrementing by 1
a <- mclapply(arr, bfastFunction, mc.cores = no_cores) #parallelize running bfastFunction with each number in arr as a count parameter

# Now recombine all individual files of each type

filename=paste0(tileName, "_trend_breaks_time.txt")
system( paste0( "cat ",dir,"/out-tile",tileName,"-pid*/", filename, " > Output/",filename) )

filename=paste0(tileName, "_trend_breaks_magnitude.txt")
system( paste0( "cat ",dir,"/out-tile",tileName,"-pid*/", filename, " > Output/",filename) )

filename=paste0(tileName, "_trend_nbbreaks.txt")
system( paste0( "cat ",dir,"/out-tile",tileName,"-pid*/", filename, " > Output/",filename) )

filename=paste0(tileName, "_trend_bfast.txt")
system( paste0( "cat ",dir,"/out-tile",tileName,"-pid*/", filename, " > Output/",filename) )

filename=paste0(tileName, "_season_nbbreaks.txt")
system( paste0( "cat ",dir,"/out-tile",tileName,"-pid*/", filename, " > Output/",filename) )

filename=paste0(tileName, "_season_breaks_time.txt")
system( paste0( "cat ",dir,"/out-tile",tileName,"-pid*/", filename, " > Output/",filename) )

filename=paste0(tileName, "_season_bfast.txt")
system( paste0( "cat ",dir,"/out-tile",tileName,"-pid*/", filename, " > Output/",filename) )

# Remove all the individual files
system( paste0( "rm -fr ",dir,"/out-tile",tileName,"-pid*") )


rm(list=ls())

Step 3: Rebuild Raster Outputs

The BFAST outputs include:

Trend component
Seasonal component
Number of breaks
Break timing
Break magnitude
Slope

Rebuilding Workflow

1. Join BFAST results to the original point shapefile.

2. Use Point to Raster to create rasters for:

Trend
Seasonal component
Number of breaks
Slope

3. Mosaic tile rasters into a single raster for each variable.

Final outputs include country-wide rasters representing vegetation change dynamics.

BFAST Analysis

Michael Pegg

2026-03-04

Breaks for Additive Season and Trend (BFAST) Walkthrough

Overview

Breaks for Additive Season and Trend (BFAST) is a parametric time-series analysis method that applies an additive decomposition model to separate each pixel’s time series into three components:

The primary objectives of BFAST are:

1. Detect multiple abrupt changes in the trend and seasonal components.

2. Characterize gradual and abrupt ecosystem changes by estimating the timing, magnitude, and direction of change.

Statistical significance of trends is assessed using a Student’s t-test to determine whether the slope significantly differs from zero. Based on slope sign and significance, vegetation change is classified as:

Data Acquisition

Data were acquired from the NASA AppEEARS website.

Parameters Used

Rescaling EVI Values

MODIS EVI values range from -2000 to 10,000.

To convert to standard EVI values (-1 to 1), multiply by 0.0001.

Fishnet Creation

Fishnet tiles were created across Ukraine to reduce computational demand.

Workflow

1. Convert MODIS EVI imagery to Integer data type.

2. Use Raster to Points to extract center points.

3. Use Extract Multi Values to Points to add EVI values.

4. Export attribute table to .csv.

BFAST Analysis

Step 1: Handle Missing Values

BFAST will not run if NA values are present.

Ensure all missing values are filled or removed before analysis.

Step 2: Parallel BFAST Processing

Parallel processing was performed to reduce computation time. Each tile generates multiple output files containing trend and seasonal break information.

Step 3: Rebuild Raster Outputs

The BFAST outputs include:

Rebuilding Workflow

1. Join BFAST results to the original point shapefile.

2. Use Point to Raster to create rasters for:

3. Mosaic tile rasters into a single raster for each variable.

Final outputs include country-wide rasters representing vegetation change dynamics.