BFAST分析MODIS时间序列数据

前面分享了一篇论文，打算复现研究一下，查了查用到的方法，在这里分享Jan Verbeselt写的教程。教程原文请大家点击文末 阅读原文 跳转。

经过测试代码能够运行通过。由于需要下载MODIS数据，运行时间比较长，我是代码跑起来之后出去吃了个晚饭才运行完的。

【文献阅读】长江流域植被变化监测与潜在驱动力1982-2015

基于R的MODIS数据下载和分析

我们使用R语言从Loobos Site下载MODIS数据。在此之前你需要下载安装R和Rstudio。

安装和加载MODIS数据分析所需的R包

在加载包之前，需要先安装R包，使用下面的代码可以把所需的R包都安好。

# 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("strucchange","zoo", "bfast", "raster", "leaflet", "MODISTools")
for (package in neededPackages){pkgTest(package)}

定义了一个timeser函数用于创建时间序列。

# Function to create time series object
# 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)
  }

使用MODISTools包下载MODIS数据

# 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",
                km_lr = 2,
                km_ab = 2,
                site_name = "testsite",
                internal = TRUE,
                progress = FALSE)

# 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",
                km_lr = 2,
                km_ab = 2,
                site_name = "testsite",
                internal = TRUE,
                progress = FALSE)

创建栅格砖（Raster Brick）清洗MODIS数据

使用mt_to_raster函数（MODISTools包中）创建栅格砖。

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

使用像元可靠性信息清洗MODIS NDVI数据

## 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

使用BFAST监测

下面将数据从栅格中取出，作为一个向量，使用timeser函数生成一个时间序列。

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

使用bfastmonitor函数，利用trend + harmon模型，order 3作为谐波（季节模型）

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

接下来就是栅格断点计算：

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

# here we define the function that we will apply across the brick using the calc 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)) 
}

# calc function 
bfmR <- calc(VI_m, bfmRaster)
names(bfmR) <- c('time of break', 'magnitude of change')
plot(bfmR) # resulting time and magnitude of change

使用click函数选择像元，并对像元进行BFAST分析：

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

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

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

参考文献

1. https://verbe039.github.io/BFASTforAEO/
2. https://vip.arizona.edu/documents/MODIS/MODIS_VI_UsersGuide_June_2015_C6.pdf
3. https://bfast2.github.io/