时间序列的栅格数据Sen+MK分析（R语言）源自徐洋——NDVI时间序列分析之Sen+MK分析全过程梳理

一、代码运行文件夹内为年数据或月数据，
得三个数据：

MKtest（用于检验，不用于制图）

slope(趋势)

Z（显著度）

#install.packages("sp")
#install.packages("raster")
#install.packages("rgdal")
#install.packages("trend")
#install.packages("terra")
#dir.create("D:/rfiles")
#setwd("D:/rfiles")
#getwd()
 
 
#pkgbuild::find_rtools(debug = TRUE)
 
 
library(sp)
library(raster)
library(rgdal)
library(trend)
library(terra)
setwd("J:/PTJPL ET8218/驱动因子/EVI/year") #输入一个文件夹内的单波段TIFF数据，在这里是历年的NDVI年最大值
fl <- list.files(pattern='*tif$')
firs<-raster(fl[1])
for (i in 2:37) {
  r <- raster(fl[i])
  firs <- stack(firs, r)
}
fun <- function(y){
  if(length(na.omit(y)) < 37) return(c(NA, NA, NA))   #删除数据不连续含有NA的像元
  MK_estimate <- sens.slope(ts(na.omit(y), start = 1982, end = 2018, frequency = 1), conf.level = 0.95) #Sen斜率估计
  slope <- MK_estimate$estimate
  MK_test <- MK_estimate$p.value
  Zs <- MK_estimate$statistic
  return(c(Zs, slope, MK_test))
}
 
e <- calc(firs, fun)     #栅格计算
e_Zs <- subset(e,1)      #提取Z值
e_slope <- subset(e,2)   #提取sen斜率
e_MKtest <- subset(e,3)  #提取p值
 
plot(e_Zs)
plot(e_slope)
plot(e_MKtest)
 
writeRaster(e_Zs, "J:/PTJPL ET8218/驱动因子/senEVI/EVI1982_2018_Zs.tif", format="GTiff", overwrite=T)
writeRaster(e_slope, "J:/PTJPL ET8218/驱动因子/senEVI/EVI1982_2018_slope.tif", format="GTiff", overwrite=T)
writeRaster(e_MKtest, "J:/PTJPL ET8218/驱动因子/senEVI/EVI1982_2018_MKtest.tif", format="GTiff", overwrite=T)

二、充分类 变化趋势划分

结合SNDVI和Z统计量划分NDVI变化趋势：

• slope

  • -0.0005~0.0005稳定区域

  • 大于或等于0.0005植被改善区域

  • 小于-0.0005为植被退化区域

• Z统计量

  Slope划分

  • 置信水平0.05

  • Z绝对值大于1.96显著

  • Z绝对值小于等于1.96不显著

• Slope被划分为三级：

  • SNDVI≤−0.0005 植被退化

  • −0.0005≥SNDVI≥0.0005 植被生长稳定

  • SNDVI≥0.0005 植被改善

• 使用重分类（Reclassify）对slope进行划分

  • 由于slope.tif文件研究区范围外的值非空，所以在这里先裁剪了一下

  • 裁剪所用矢量和栅格数据坐标系需要一致，否则范围容易出错

  • 统一使用了WGS84地理坐标系作为空间参考

  • 使用Model builder构建地理处理流

Slope划分过程

• 重分类结果：

  • -1退化

  • 0稳定

  • 1改善

Z值划分

• 对Z值进行重分类，确定显著性

  • |Zs|≤1.96 未通过95%置信度检验，不显著

  • |Zs|≥1.96 通过95%置信度检验，显著

Z值重分类

• 重分类结果：

  • 1不显著

  • 2显著

三、变化趋势计算

使用栅格计算器将Slope和Z值计算结果相乘，最后得到趋势变化划分

• -2严重退化

• -1轻微退化

• 0稳定不变

• 1轻微改善

• 2明显改善