一、代码运行文件夹内为年数据或月数据,
得三个数据:
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明显改善