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

水野真纪

已于 2024-06-17 14:06:37 修改

阅读量2k

点赞数 12

分类专栏： R 文章标签： r语言 python 开发语言

于 2024-06-16 21:51:48 首次发布

本文链接：https://blog.csdn.net/weixin_71625923/article/details/139727046

版权

R 专栏收录该内容

2 篇文章

订阅专栏

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

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构建地理处理流