BFAST分析MODIS时间序列数据
前面分享了一篇论文，打算复现研究一下，查了查用到的方法，在这里分享Jan Verbeselt写的教程。教程原文请大家点击文末 阅读原文 跳转。
经过测试代码能够运行通过。由于需要下载MODIS数据，运行时间比较长，我是代码跑起来之后出去吃了个晚饭才运行完的。
基于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)
参考文献
https://verbe039.github.io/BFASTforAEO/ R语言安装部署基础 https://vip.arizona.edu/documents/MODIS/MODIS_VI_UsersGuide_June_2015_C6.pdf https://bfast2.github.io/
