2

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.

0

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")
  • 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 Sep 19 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. – Mikkel Lydholm Rasmussen Sep 20 at 8:31
  • Sorry for bothering you again! please see my answer with revised code. I couldn't still run the code. Can you please notice my mistake. Thanks again – SA Khan Sep 22 at 20:12
0

As commented, my rasterstack of 597 ndvi images is:

chitral_ROI_4

create a regular daily time series object by combining data and date information

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

s.d.linear <- round(na.approx(s),0)  

Create time series to data frame function

time.series.to.dataframe <- function(time_series, source) {

  s.df           <- data.frame(as.numeric(time_series))
  colnames(s.df) <- "NDVI"
  s.df$Time      <- as.Date(date_decimal(as.numeric(time(time_series))))
  s.df$Type      <- source
  return(s.df)

}

Apply function time series to dataframe

s.d.linear.df <- time.series.to.dataframe(s.d.linear, "Linear Interpolation")

s.d.periodic <- round(na.interp(s),0)
plot(s.d.periodic)

Now aggregated from daily to monthly time series

aggregate.daily.to.monthly <- function(daily.ts) {

  s.month <- round(aggregate(as.zoo(daily.ts), as.yearmon, median), 0)   
  s.month <- as.ts(s.month)

  return(s.month)

}

s.m.linear   <- aggregate.daily.to.monthly(s.d.linear)
s.m.periodic <- aggregate.daily.to.monthly(s.d.periodic) 

Now from here onwards, I adapt to @Mikkel's code by taking out ts() in xbfst function as I have s.m.periodic as a regular monthly time series:

xbfast <- function(s.m.periodic) {  
  result <- bfastts(s.m.periodic, 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)
}

Now I replace vegevi with my own rasterstack and apply bfast_output which results in the following error.

bfast_output <- calc(chitral_ROI_4, fun=bfastfun) 
Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) : 
  cannot use this function

What could be wrong here? Any ideas please.

  • You main issue lies in the way you treat your rasters - as.vector() is not a good thing to use on a rasterStack or rasterBrick. You can exemplify this by checking the type of 's.m.periodic'. My example code expects a rasterBrick, but it is being provided with something else - likely a ts-object. Additionally, you seem to have combined bfastts and bfast into one command, which doesn't work. – Mikkel Lydholm Rasmussen Sep 23 at 14:11
  • You are right. I am providing it with a ts-object which is not working. – SA Khan Sep 25 at 17:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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