前面分享了一篇論文，打算復現研究一下，查了查用到的方法，在這裡分享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 maskplot(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/https://vip.arizona.edu/documents/MODIS/MODIS_VI_UsersGuide_June_2015_C6.pdfhttps://bfast2.github.io/