3

I am working on identification of breakpoints in Landsat based NDVI time series (1990-2018). I am using The Breaks for Additive Seasonal and Trend (BFAST) algorithm which is a pixel based time series decomposition method to detect and characterise abrupt changes. I wonder would it be possible to apply it for a region of interest or study area as a whole instead of applying it on selected pixels?

ndvi_stack <- stack(ndvi_list)
ndvi_ROI_1 <- crop(ndvi_stack, extent(ROI_1))
ndvi_ROI_1

So, ndvi_ROI_1 is a study area or region of interest. Now I select the pixel:

selected_pixel <- 100

I created an irregular time series for the selected pixel (100 in this case) using bfastts() from BFAST package.

(s <- bfastts(as.vector(ndvi_ROI_1[selected_pixel]), dates, type = c("irregular")))

The time series is converted to a regular time series using advice from https://philippgaertner.github.io/2018/04/bfast-preparation/.

I also tried not to select to a specific pixel and run bfastts() on whole study region and it seems to work (although I am not sure about the validity of results). I used the following code:

(s <- bfastts(as.vector(ndvi_ROI_1), dates, type = c("irregular")))

I am still not sure if this is technically right to apply bfast method to a whole study region instead of pixel by pixel basis? And if it is, is the coding correct or I have to introduce any other steps?

Looking forward to the kind advice.

1 Answer 1

1

I used BFAST for my thesis some years ago, and this is the code I used. Given the age of this code, things may have changed significantly, to make specific elements not function entirely any more, but the process of creating a function to apply BFAST and then applying that function on the raster is still a good way to approach the problem. Do note that the processing time can become lengthy.
It should be noted that I preprocessed my data into a regular smooth time-series before utilizing R.

One key thing that seems to have changed may be that the function that used to be named ts is now bfastts, but I haven't checked if the syntax is the same as it was.

library(bfast)
library(raster)
library(gtools)
vegevi <- brick("evi_subset_brick.grd") # loads my regular EVI data into a raster-brick. EVI is just a different vegetation index, that is more suited for areas with a high vegetation density than NDVI.

xbfast <- function(data) {  
  ndvi <- ts(data, frequency=46, start=1) # adjust the frequency to match yours
  result <- bfast(ndvi, season="harmonic", max.iter=2, breaks=1)
  niter <- length(result$output)
  out <- result$output[[niter]]
  bp <- out$Wt.bp #breakpoint of the seasonality component
  st <- out$St #the seasonality component
  st_a <- st[1:bp] #seasonality untill the breakpoint
  st_b <- st[bp:460] #hard coded end-point - should likely be changed to something better
  st_amin <- min(st_a)
  st_amax <- max(st_a)
  st_bmin <- min(st_b)
  st_bmax <- max(st_b)
  st_adif <- st_amax - st_amin
  st_bdif <- st_bmax - st_bmin
  st_dif <- st_bdif - st_adif
  Magni<-result$Magnitude #magnitude of the biggest change detected in the trend component
  Timing<-result$Time #timing of the biggest change detected in the trend component
  return(c(st_dif,bp,Magni,Timing)) 
}

bfastfun <- function(y) {
  percNA <- apply(y, 1, FUN=function(x) (sum(is.na(x))/length(x)) )
  i <- (percNA<0.2)
  res <- matrix(NA, length(i), 4) #
  if (sum(i) > 0) {
    res <- t(apply(y[i,], 1, xbfast))
  }
  return(res)
}

bfast_output <- calc(vegevi, fun=bfastfun)

diff <- subset(bfast_output, 1)
time <- subset(bfast_output, 2)
magn <- subset(bfast_output, 3)
magn_time <- subset(bfast_output, 4)

writeRaster(diff, filename="evi_diff.tif", format="GTiff")
writeRaster(time, filename="evi_diff_time.tif", format="GTiff")
writeRaster(magn, filename="evi_magni.tif", format="GTiff")
writeRaster(magn_time, filename="evi_magni_timing.tif", format="GTiff")
2
  • Dear Mikkel, thank you very much for sharing your code. My time series data is irregular therefore I can't use ts() function. I have to convert irregular time series to a regular first and then aggregate to regular monthly time series. My aggregated monthly time series is: s.m.periodic <- aggregate.daily.to.monthly(s.d.periodic) I am not sure how can I amend your code to complete remaining steps. Should I run xbfast & bfastfun functions as it is & only replace vegevi in bfast_output with my raster brick? Appreciate your help.
    – SA Khan
    Commented Sep 19, 2019 at 20:23
  • @SAKhan - yes, if you have a monthly time-series, then you can go right ahead and replace my 'vegevi' with your own raster brick. Commented Sep 20, 2019 at 8:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.