Terra中基于另一栅格重分类大型RasterStack的计数异常问题

需求说明 #

• 手里有月度温度栅格栈hist_temp_stack，要完成两个操作：
  1. 计算每个栅格单元的90分位数T90
  2. 统计每个栅格在所有图层中达到或超过该分位数的次数，得到频次汇总栅格，用于对比过去与未来的阈值超频次变化

初始代码踩坑 #

一开始尝试的代码：

hist_extreme_temp <- app(hist_temp_stack, fun = function(x) ifelse(x >= T90[], 1, 0))
hist_extreme_temp_count <- sum(hist_extreme_temp)

出现两个问题：

• 弹出警告：longer object length is not a multiple of shorter object length
• 汇总栅格最大值达62977，远超过图层总数492，结果完全不合理

后来用T90 <- rep(T90, nlyr(hist_temp_stack))尝试修复警告，直接触发内存不足错误（提示需8TB内存）

问题根源 #

• 使用app()时，每个栅格单元的所有图层值会以向量形式传入函数，但T90是单图层栅格，直接比对时R会强制循环T90的长度匹配向量，导致逻辑混乱，统计结果出错
• 用rep()复制T90会生成与整个栅格栈总像素数相同的超大对象，直接耗尽内存

正确解决方案 #

换用calc()函数，直接在每个栅格单元的时间序列内完成分位数计算与频次统计，无需额外复制分位数栅格：

hist_extreme_temp_count <- calc(hist_temp_stack, fun = function(x) {
  # 计算当前栅格单元的90分位数（忽略NA值）
  t90 <- quantile(x, probs = 0.9, na.rm = TRUE)
  # 统计该单元所有图层中≥分位数的次数
  sum(x >= t90, na.rm = TRUE)
})

方案优势 #

• 无需复制分位数栅格，内存占用极低，不会出现内存不足问题
• 每个栅格单元独立计算，逻辑匹配准确，无长度不匹配警告
• 一步完成两个操作，代码简洁高效

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