前面的推文说过Sen MK分析

• NDVI分析时间序列Sen MK整个过程分析梳理

上述方法的一个很大缺点是R语言计算过程太慢，我改进了前面的代码，并使用了它terra包重写，打开并行计算，提高计算速度。

• rast读取网格数据
• function中长度和年份需要根据实际情况进行修改
• 输出结果包括三个波段，Z值、slope和p值(详见参考文献推文)

library(terra) library(trend) #输入文件夹中的单波段TIFF这里有历年的数据NDVI年最大值 flnames <- list.files(path = './ChinaYearMean/', pattern = '.tif$') fl <- paste0("./ChinaYearMean/", flnames) firs <- rast(fl)  #Sen MK计算 fun_sen <- function(x){   if(length(na.omit(x)) <34) return(c(NA, NA, NA))   #删除数据不连续NA的像元   MK_estimate <- trend::sens.slope(ts(na.omit(x), start = 1982, end = 2015, frequency = 1), conf.level = 0.95) #Sen斜率估计  #部分代码省略，完整代码请关注微信微信官方账号，搜索同名文章获取 #部分代码省略，完整代码请关注微信微信官方账号，搜索同名文章获取 #部分代码省略，完整代码请关注微信微信官方账号，搜索同名文章获取   return(c(Zs, slope, MK_test)) }  #部分代码省略，完整代码请关注微信微信官方账号，搜索同名文章获取  #显示输出 plot(firs_sen) writeRaster(firs_sen, filename = "./firs_sen.tif", names=firs_sen@ptr[["names"]]) 

1. 对比R语言Raster包和Terra包栅格计算
2. NDVI分析时间序列Sen MK整个过程分析梳理
3. NDVI时间序列逐像元计算Hurst指数