如何基于共同索引对两个Raster Stack执行匹配运算？

方法1：重复参考栈后直接运算（最优解） #

这是最简洁高效的方式，利用terra内置的rep()函数将参考栈Stack2重复至与Stack1相同的层数，随后直接执行元素级除法：

# 将Stack2重复3次，使其层数与Stack1（18层）一致
Stack2_rep <- rep(Stack2, times = 3)
# 执行对应月份的除法：Stack1的每个月层自动匹配Stack2_rep的对应月层
result <- Stack1 / Stack2_rep

terra会自动优化这种向量化运算，无需额外分组或循环，性能最优。

方法2：使用tapp实现分组运算 #

若坚持使用tapp，可以结合lapply遍历月份分组，对每个分组内的Stack1图层匹配Stack2的对应月份层：

idx <- rep(1:6, 3)
# 定义运算函数：分组内的Stack1子栈 ÷ Stack2对应月份层
monthly_div <- function(x, ref, month_idx) {
  x / ref[[month_idx]]
}

# 遍历每个唯一月份索引，处理后合并结果
result_list <- lapply(unique(idx), function(i) {
  # 提取Stack1中所有属于第i个月的图层并运算
  tapp(Stack1, idx, function(x) monthly_div(x, Stack2, i), which = i)
})
result <- do.call(c, result_list)

方法3：显式子集匹配（直观易读） #

如果需要更直观的分组逻辑，可以直接按索引拆分图层后运算：

idx <- rep(1:6, 3)
result_list <- list()

# 遍历Stack2的每个月份层
for(i in seq_len(nlyr(Stack2))) {
  # 提取Stack1中所有对应第i个月的图层
  s1_month_layers <- Stack1[[idx == i]]
  # 执行除法并存入结果列表
  result_list[[i]] <- s1_month_layers / Stack2[[i]]
}

# 将所有结果图层合并为一个栈
result <- do.call(c, result_list)

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