R语言计算sen趋势

library(terra)
library(trend)

# 设置输入路径（NDVI数据存放位置）
input_path <- "G:/影响因子驱动分析/NDVIdata/ndvi/"

# 设置输出路径（斜率结果保存位置）
output_path <- "G:/论文/lucc/00-20ndvi/"

# 读取NDVI时间序列数据
fl <- list.files(path = input_path, 
                pattern = ".tif$", 
                full.names = TRUE)

# 按文件名排序确保时间顺序正确（重要！）
fl <- fl[order(as.numeric(gsub("\\D", "", fl)))] # 提取文件名中的数字排序
firs <- rast(fl)

# 定义Sen斜率计算函数
fun_slope <- function(x) {
  if (sum(!is.na(x)) < 21) return(NA)  # 21年数据完整性检查
  ts_data <- na.omit(x)
  if (length(ts_data) < 2) return(NA)  # 至少需要2个有效值
  MK_estimate <- trend::sens.slope(ts_data)
  return(MK_estimate$estimate)
}

# 并行计算(4核)
firs_slope <- app(firs, fun_slope, cores = 4)
names(firs_slope) <- "slope"

# 创建输出目录（如果不存在）
if (!dir.exists(output_path)) {
  dir.create(output_path, recursive = TRUE, showWarnings = FALSE)
}

# 保存结果
writeRaster(firs_slope,
            filename = paste0(output_path, "NDVI_Sen_Slope_2000-2020.tif"),
            overwrite = TRUE,
            datatype = "FLT4S")  # 使用浮点型保存小数

# 可视化结果
plot(firs_slope, 
     main = "NDVI变化趋势 (2000-2020)",
     col = colorRampPalette(c("red", "white", "green"))(100))