如何用terra包对分组栅格栈批量计算中位数与分位数？

解决方案 #

要实现按模型分组计算栅格的中位数、第5和第95百分位数，关键是在terra::app函数中正确传递分组信息并对每个单元格的数值进行分组统计。以下是可行代码：

library(terra)

# 创建测试数据
a <- rast(ncol = 10, nrow = 10, vals=rep(5,100), names=1)
b <- rast(ncol = 10, nrow = 10, vals=rep(10,100), names=1)
c <- rast(ncol = 10, nrow = 10, vals=rep(5,100), names=2)
d <- rast(ncol = 10, nrow = 10, vals=rep(10,100), names=2)
z <- c(a, b, c, d)

# 基于图层名称定义分组
groups <- as.integer(names(z))
groups_unique <- unique(groups)

# 单组统计函数（含NA值处理）
compute_stats <- function(x) {
  c(
    median(x, na.rm = TRUE),
    quantile(x, 0.05, na.rm = TRUE),
    quantile(x, 0.95, na.rm = TRUE)
  )
}

# 单元格级分组统计函数
grouped_stats <- function(x) {
  # 按分组拆分当前单元格的数值
  grouped_data <- split(x, groups)
  # 对每个分组计算统计量
  stats_results <- lapply(grouped_data, compute_stats)
  # 扁平化结果向量
  flat_results <- unlist(stats_results)
  # 为结果分配有意义的名称
  names(flat_results) <- paste0(
    rep(c("median_", "q5_", "q95_"), length(groups_unique)),
    rep(groups_unique, each = 3)
  )
  flat_results
}

# 执行计算
g <- app(z, fun = grouped_stats)

# 查看结果
g

关键说明 #

1. 分组定义：提前基于栅格图层名称创建分组向量，确保在app函数中能正确拆分每个单元格的数值。
2. NA值处理：在统计函数中加入na.rm = TRUE，避免缺失值干扰计算。
3. 结果命名：通过拼接前缀和分组ID，生成符合预期的图层名称（如median_1、q5_2）。

执行后将得到包含6层的栅格栈，完全匹配你期望的输出格式。

内容的提问来源于stack exchange，提问作者TheRealJimShady