在R语言中计算栅格几何均值的正确语法及NA值处理方法

关于(raster1*raster2)^(1/2)的适用性 #

这个代码在两个栅格对应位置均无NA时是有效的，但只要其中一个栅格的对应位置为NA，计算结果就会变成NA（R中NA参与算术运算的结果默认是NA）。如果你的数据不存在NA，或者可以接受NA直接传递到结果中，那这个代码完全可以用；但如果需要忽略NA值计算几何均值，这个方法就不适用了。

处理NA值的方法 #

单/双栅格场景 #

使用raster包的calc函数结合几何均值的对数转换公式（几何均值的对数等于各值对数的均值），通过na.rm=TRUE参数忽略NA值：

library(raster)
# 堆叠两个栅格
r_stack <- stack(raster1, raster2)
# 计算忽略NA的几何均值
geo_mean <- calc(r_stack, fun = function(x) {
  exp(mean(log(x), na.rm = TRUE))
})

注：如果某个位置的所有栅格值都是NA，结果仍会是NA；若仅部分为NA，则会基于非NA值计算几何均值。

批量栅格场景（替代list.files+for循环） #

如果是批量处理用list.files读取的栅格文件，直接将所有栅格堆叠后用calc计算，比for循环更高效简洁：

# 读取指定文件夹下的所有tif栅格
raster_paths <- list.files(path = "你的栅格文件路径", pattern = "\\.tif$", full.names = TRUE)
# 堆叠所有栅格
r_stack <- stack(raster_paths)
# 批量计算忽略NA的几何均值
geo_mean <- calc(r_stack, fun = function(x) {
  exp(mean(log(x), na.rm = TRUE))
})

注意事项 #

如果栅格中存在0或负值，log()运算会报错，此时需要提前处理（比如过滤负值，或给0添加一个极小的正数）——因为几何均值仅适用于正值数据。

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