MODIS based time series analysis using BFAST Monitor and BFAST Lite

Introduction to MODIS data and online analysis

The MODIS satellite data that we will download for this time series analysis exercise is available from the MODIS Subsetting Website. MODIS data is made available for specific locations.

More specifically, we will look at the MODIS product called MOD13Q1 which are global 16-day images at a spatial resolution of 250 m. Each image contains several bands; i.e. blue, red, and near-infrared reflectances, centered at 469 nanometers, 645 nanometers, and 858 nanometers, respectively. These bands are used to determine the MODIS vegetation indices.

The MODIS Normalized Difference Vegetation Index (NDVI) complements NOAA’s Advanced Very High Resolution Radiometer (AVHRR) NDVI products and provides continuity for applications requiring long time series. MODIS also includes a new Enhanced Vegetation Index (EVI) that minimizes canopy background variations and maintains sensitivity over dense vegetation conditions. The EVI also uses the blue band to remove residual atmosphere contamination caused by smoke, sub-pixel and thin clouds. The MODIS NDVI and EVI products are computed from atmospherically corrected surface reflectances that have been masked for water, clouds, heavy aerosols, and cloud shadows.

Vegetation indices are used for global monitoring of vegetation conditions and are used in products displaying land cover and its changes. We will work with the MODIS NDVI band of the MOD13Q1 product. More information about this MODIS data set can be found via the MODIS product table.

Now go to the MODIS Land Products:

• Search for ‘Netherlands’ and select the ‘Gelderland Loobos’ Site.
• Select the MODIS 13Q1 product by clicking on it.There is no immediate need to create an account unless you want to manually download the MODIS data via this website. In the below explanation, we will download MODIS NDVI data via an R script.
• Under ‘Download data’ tab and table, via the website, you can see a small link named Pixel Numbering Scheme, Click it to see that the blue pixel (i.e. pixel number 545) is the “site pixel” (i.e. the middle of the spatial subset). This is where a flux tower is located.

Preprocessing: automated MODIS NDVI download and analysis via R

We will download the MODIS data for the Loobos Site via R and process the data for one location to detect changes within the time series.

getwd() # check if your working directory is correctly set

# pkgTest is a helper function to load packages and install packages only when they are not installed yet.
pkgTest <- function(x)
{
  if (x %in% rownames(installed.packages()) == FALSE) {
    install.packages(x, dependencies= TRUE)
  }
  library(x, character.only = TRUE)
}
neededPackages <- c("zoo", "bfast", "terra", "raster", "leaflet", "MODISTools")
for (package in neededPackages){pkgTest(package)}

# Utility function to create time series object from a numeric vector
# val_array: data array for one single pixel (length is number of time steps)
# time_array: array with Dates at which raster data is recorded (same length as val_array)
timeser <- function(val_array, time_array) {
    z <- zoo(val_array, time_array) # create zoo object
    yr <- as.numeric(format(time(z), "%Y")) # extract the year numbers
    jul <- as.numeric(format(time(z), "%j")) # extract the day numbers (1-365)
    delta <- min(unlist(tapply(jul, yr, diff))) # calculate minimum time difference (days) between observations
    zz <- aggregate(z, yr + (jul - 1) / delta / 23) # aggregate into decimal year timestamps
    (tso <- as.ts(zz)) # convert into timeseries object
    return(tso)
  }

# Downloading the NDVI data, starting from 2000-01-01
VI <- mt_subset(product = "MOD13Q1",
                site_id = "nl_gelderland_loobos",
                band = "250m_16_days_NDVI",
                start = "2000-01-01",
                end = "2022-03-22",
                km_lr = 2,
                km_ab = 2,
                site_name = "testsite",
                internal = TRUE,
                progress = TRUE)

## Downloading chunks:

# Downloading the pixel reliability data, starting from 2000-01-01
QA <- mt_subset(product = "MOD13Q1",
                site_id = "nl_gelderland_loobos",
                band = "250m_16_days_pixel_reliability",
                start = "2000-01-01",
                end = "2022-03-22",
                km_lr = 2,
                km_ab = 2,
                site_name = "testsite",
                internal = TRUE,
                progress = TRUE)

## Downloading chunks:

# convert df to raster
VI_r <- mt_to_terra(df = VI)
QA_r <- mt_to_terra(df = QA)

## clean the data
# create mask on pixel reliability flag set all values <0 or >1 NA
m <- QA_r
m[(QA_r < 0 | QA_r > 1)] <- NA # continue working with QA 0 (good data), and 1 (marginal data)

# apply the mask to the NDVI raster
VI_m <- mask(VI_r, m, maskvalue=NA, updatevalue=NA)

# plot the first image
plot(m,1) # plot mask

plot(VI_m,1) # plot cleaned NDVI raster

# extract data from the cleaned raster for selected pixels
click(VI_m, id=TRUE, xy=TRUE, cell=TRUE, n=1)

library(leaflet)
r <- raster(VI_m[[1]]) # Select only the first layer (as a RasterLayer)
pal <- colorNumeric(c("#ffffff", "#4dff88", "#004d1a"), values(r),
  na.color = "transparent")

map <- leaflet() %>% addTiles() %>%
  addRasterImage(r, colors = pal, opacity = 0.8) %>%
  addLegend(pal = pal, values = values(r),
    title = "NDVI")
map

## check VI data at a certain pixel e.g. 1 row, complete left hand site:
## the dimensions of the raster are: 33x33

px <- 78 # pixel number; adjust this number to select the center pixel
tspx <- timeser(unlist(VI_m[px]),as.Date(names(VI_m), "%Y-%m-%d")) # convert pixel "px" to a time series
plot(tspx, main = 'NDVI') # NDVI time series cleaned using the "reliability information"

Option 1: detect break at the end of the time series with BFAST Monitor

Now we apply the bfastmonitor function using a trend + harmon model with order 3 for the harmonics (i.e. seasonality modelling):

bfm1 <- bfastmonitor(tspx, response ~ trend + harmon, order = 3, start = c(2018,3)) # Note: the third observation in 2018 marks the transition from 'history' to 'monitoring'
plot(bfm1)

Let’s run the bfastmonitor code on the full raster time series spatially using the app function:

dates <- as.Date(names(VI_m), "%Y-%m-%d")

# here we define the function that we will apply across the brick using the `app` function:
bfmRaster = function(pixels)
{
    tspx <- timeser(pixels, dates) # create a timeseries of all pixels
    bfm <- bfastmonitor(tspx, response ~ trend + harmon, order = 3, start = c(2019,1)) # run bfast on all pixels
    return(c(bfm$breakpoint, bfm$magnitude)) 
}

# apply the function to each raster cell (time series)
# Optionally you can supply an argument cores=n (where n is the number of cores on your computer) for a potential speed boost
bfmR <- app(VI_m, bfmRaster)
names(bfmR) <- c('time of break', 'magnitude of change')
plot(bfmR) # resulting time and magnitude of change

Here is example R code to get the pixel number. When you run click(), click in the plot on a pixel with a break (i.e. where an estimated time of break is available).

plot(bfmR,1)
click(VI_m, id=FALSE, xy=FALSE, cell=TRUE, n=1)

Here we selected one pixel, and do the bfastmonitor analysis for that pixel:

px <- 460 # pixel number so adjust this number to select the center pixel
tspx <- timeser(unlist(VI_m[px]),as.Date(names(VI_m), "%Y-%m-%d")) # convert pixel 1 to a time series
plot(tspx, main = 'NDVI') # NDVI time series cleaned using the "reliability information"

tspx[tspx < 0] <- NA
bfm <- bfastmonitor(tspx, response ~ trend + harmon, order = 3, start = c(2019,2))
plot(bfm)

Option 2: detecting all breaks in the middle of a time series with BFAST Lite

BFAST Monitor is made to detect the first break at the end of a time series. If you need to detect more than one break, then you need to use a different algorithm: either BFAST or BFAST Lite. In this example, we will run BFAST Lite (a lightweight version of BFAST that can handle NA values) on the same data as we did above.

Let’s run the function bfastlite() on our data from the previous steps. Setting the parameter breaks to BIC results in more liberal detection of breaks compared to the default, LWZ.

breaks <- bfastlite(tspx, response ~ trend + harmon, order = 3, breaks = "BIC")
breaks

## 
##   Optimal 1-segment partition: 
## 
## Call:
## breakpoints.matrix(obj = X, y = y, h = h, breaks = breaks, hpc = hpc)
## 
## Breakpoints at observation number:
## NA 
## 
## Corresponding to breakdates:
## NA

What if we try another pixel:

px <- 82
tspx <- timeser(unlist(VI_m[px]),as.Date(names(VI_m), "%Y-%m-%d"))
breaks <- bfastlite(tspx, response ~ trend + harmon, order = 3, breaks = "BIC")
breaks

## 
##   Optimal 2-segment partition: 
## 
## Call:
## breakpoints.matrix(obj = X, y = y, h = h, breaks = breaks, hpc = hpc)
## 
## Breakpoints at observation number:
## 342 
## 
## Corresponding to breakdates:
## 0.753304

Which date is the break on?

dates.no.na <- as.numeric(time(tspx))
dates.no.na[is.na(tspx)] <- NA
dates.no.na <- na.omit(dates.no.na)
dates.no.na[breaks$breakpoints$breakpoints[1]]

## [1] 2016.696

Plot the model and the break(s):

plot(breaks)

Why did the model decide to place one break? Let’s take a look at the statistics:

summary(breaks$breakpoints)

## 
##   Optimal (m+1)-segment partition: 
## 
## Call:
## breakpoints.matrix(obj = X, y = y, h = h, breaks = breaks, hpc = hpc)
## 
## Breakpoints at observation number:
##                           
## m = 1                  342
## m = 2              257 342
## m = 3      142 221     342
## m = 4   83 151 221     342
## m = 5   68 136 204 274 342
## 
## Corresponding to breakdates:
##                                                                                
## m = 1                                                                          
## m = 2                                                         0.566079295154185
## m = 3                     0.312775330396476 0.486784140969163                  
## m = 4   0.182819383259912 0.332599118942731 0.486784140969163                  
## m = 5   0.149779735682819 0.299559471365639 0.449339207048458 0.60352422907489 
##                          
## m = 1   0.753303964757709
## m = 2   0.753303964757709
## m = 3   0.753303964757709
## m = 4   0.753303964757709
## m = 5   0.753303964757709
## 
## Fit:
##                                                                       
## m    0          1          2          3          4          5         
## RSS      2.9606     1.6347     1.5375     1.4525     1.4070     1.4232
## BIC   -941.3943 -1155.9667 -1128.7523 -1099.5009 -1058.8795  -998.6316
## LWZ   -875.7292 -1024.6364  -931.7569  -836.8403  -730.5538  -604.6408
## R.sq     0.2865     0.6060     0.6295     0.6500     0.6609     0.6570

The above shows that the model tried putting up to 5 breaks in the time series at the particular observation numbers, and selected by minimum BIC that there should be one break. To see that information visually, plot the breakpoints component:

plot(breaks$breakpoints)

The residual sum of squares keeps decreasing with more breakpoints, as that makes the model more flexible and thus fit the data better, but BIC and other information criteria apply a penalty for adding more breaks than really necessary, thus limiting the number of false detections. In this case, both BIC and LWZ agree that 1 is the optimal number of breaks.

We can visually see that the identified break is indeed quite significant. We can also look at it from a statistical point of view by using the function magnitude():

magnitude(breaks$breakpoints)

## $Mag
##         before     after       diff      RMSD      MAD        MD
## [1,] 0.6793469 0.4604267 -0.2189201 0.2209773 0.219694 -0.219694
## 
## $m.x
## [1] 342 342
## 
## $m.y
##    before     after 
## 0.6793469 0.4604267 
## 
## $Magnitude
##       diff 
## -0.2189201 
## 
## $Time
## [1] 342
## 
## attr(,"class")
## [1] "magnitude"

This shows that the difference between the actual observations and the predictions of the models on both sides of the break, when extrapolated to the other side, is fairly high (0.22 RMSD and MAD, thus at the magnitude of 0.22 NDVI units).

Let’s try to run BFAST Lite on the whole raster now. There is one problem with that: since the number of breaks is variable, we don’t have a variable number of layers. Thus, let’s plot only the break that is the most significant, so that we get two layers as output.

# The code for getting a date from above, in a function
# index is which breakpoint to list, tspx is the original time series
IndexToDate <- function(index, tspx, breaks) {
  dates.no.na <- as.numeric(time(tspx))
  dates.no.na[is.na(tspx)] <- NA
  dates.no.na <- na.omit(dates.no.na)
  dates.no.na[breaks$breakpoints$breakpoints[index]]
}

bflRaster <- function(pixels, dates, timeser, IndexToDate) {
  library(zoo)
  library(bfast)
  tspx <- timeser(pixels, dates)
  breaks <- bfastlite(tspx, response ~ trend + harmon, order = 3, breaks = "BIC")
  
   #if two breaks are recorded keep the break with the largest magnitude
  if (length(breaks$breakpoints$breakpoints)>1)
      breaks$breakpoints$breakpoints<-which.max(breaks$breakpoints$breakpoints)
  
  # If no break, return NAs
  if (is.na(breaks$breakpoints$breakpoints))
    return(c(NA,NA))
  
  # Get break with highest magnitude
  mags <- magnitude(breaks$breakpoints)
  maxMag <- which.max(mags$Mag[,"RMSD"])
  
  return(c(IndexToDate(maxMag, tspx, breaks), mags$Mag[maxMag, "RMSD"]))
}

Let’s apply the function above to the whole raster! Note that system.time() is just there to show you how long the command execution taskes. You can also safely skip it if you don’t want to know. The timeser and IndexToDate arguments are also only necessary for parallel execution, you can remove them if you only want to run the processing on one core. Run on a single core:

# This will take a while: BFAST Lite is not as fast as BFAST Monitor.
system.time({
bflR <- app(VI_m, bflRaster, dates=as.Date(names(VI_m), "%Y-%m-%d"), timeser=timeser, IndexToDate=IndexToDate)
})

##    user  system elapsed 
## 199.634   0.056 199.694

Run on 4 cores is two times faster:

# Best to use all cores! But then you have to include the definition of `timeser()` in `bflRaster`.
system.time({
bflR <- app(VI_m, bflRaster, dates=as.Date(names(VI_m), "%Y-%m-%d"), timeser=timeser, IndexToDate=IndexToDate, cores=4)
})

##    user  system elapsed 
##   2.232   0.004  64.462

Plot the results (optionally compare with BFAST Monitor results above):

names(bflR) <- c('time of break', 'magnitude of change')
plot(bflR)

More information

More information can be found on the BFAST website and in the BFAST papers mentioned on the website.

The following section gives extra information about the concept of seasonality monitoring using harmonic analysis. There are no questions from this section, it’s for your interest only,

library(bfast)
## a demo ndvi time series:
ndvi <- ts(rowSums(simts$time.series))
tsp(ndvi) <- tsp(simts$time.series)
## input variable for the sinus and cosinus functions
f <- 23
w <- 1/f
tl <- 1:length(ndvi)
## 3th order harmonic model
co <- cos(2 * pi * tl * w)
si <- sin(2 * pi * tl * w)
co2 <- cos(2 * pi * tl * w * 2)
si2 <- sin(2 * pi * tl * w * 2)
co3 <- cos(2 * pi * tl * w * 3)
si3 <- sin(2 * pi * tl * w * 3)
# fit the seasonal model using linear regression
fitm<- lm(ndvi~co+si+co2+si2+co3+si3) 
predm <- fitted(fitm) ## predict based on the modelfit
plot(co, type = "l", ylab = "cos and sin")
lines(si, type = "l", lty = 2)

#create time series bfast on the 3th order harmonic function
predm <- ts(as.numeric(predm), start=c(2000,4), frequency=23) 
plot(ndvi, lwd = 3, col = "grey", ylab = "NDVI")
lines(predm, type = "l", col = "red") # fitted