¿Puede aplicarse el algoritmo BFAST a toda una zona de estudio/referencia?

Question

Estoy trabajando en la identificación de puntos de ruptura en las series temporales de NDVI basadas en Landsat (1990-2018). Estoy utilizando el algoritmo Breaks for Additive Seasonal and Trend (BFAST), que es un método de descomposición de series temporales basado en píxeles para detectar y caracterizar los cambios bruscos. Me pregunto si sería posible aplicarlo para una región de interés o área de estudio en su conjunto en lugar de aplicarlo en píxeles seleccionados.

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

Así, ndvi_ROI_1 es un área de estudio o región de interés. Ahora selecciono el píxel:

selected_pixel <- 100

He creado una serie temporal irregular para el píxel seleccionado (100 en este caso) utilizando bfastts() del paquete BFAST.

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

La serie temporal se convierte en una serie temporal regular utilizando los consejos de https://philippgaertner.github.io/2018/04/bfast-preparation/ .

También he probado a no seleccionar un píxel concreto y ejecutar bfastts() en toda la región de estudio y parece que funciona (aunque no estoy seguro de la validez de los resultados). He utilizado el siguiente código:

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

Todavía no estoy seguro de si es técnicamente correcto aplicar el método bfast a toda una región de estudio en lugar de hacerlo píxel a píxel. Y si lo es, ¿la codificación es correcta o tengo que introducir algún otro paso?

Espero los amables consejos.

Answer 1

Utilicé BFAST para mi tesis hace algunos años, y este es el código que utilicé. Dada la antigüedad de este código, las cosas pueden haber cambiado significativamente, para hacer que elementos específicos ya no funcionen del todo, pero el proceso de crear una función para aplicar BFAST y luego aplicar esa función en el raster sigue siendo una buena manera de abordar el problema. Tenga en cuenta que el tiempo de procesamiento puede llegar a ser largo.
Hay que tener en cuenta que preprocesé mis datos en una serie temporal regular y suave antes de utilizar R.

Una cosa clave que parece haber cambiado puede ser que la función que solía llamarse ts es ahora bfastts pero no he comprobado si la sintaxis es la misma que antes.

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")

Answer 2

Como se ha comentado, mi rasterstack de 597 imágenes ndvi es:

chitral_ROI_4

crear un objeto regular de serie temporal diaria combinando información de datos y fechas

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

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

Crear series de tiempo a la función de marco de datos

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)

}

Aplicar la función series de tiempo al marco de datos

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)

Ahora se agregan las series temporales diarias a las mensuales

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) 

Ahora a partir de aquí, me adapto al código de @Mikkel quitando ts() en la función xbfst ya que tengo s.m.periodic como serie temporal mensual regular:

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)
}

Ahora reemplazo vegevi con mi propio rasterstack y aplico bfast_output lo que resulta en el siguiente error.

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

¿Qué podría estar mal aquí? Cualquier idea, por favor.