如何在Raster对象中当单元格含过多NA时将均值设为NA？

我正在处理Raster格式的环境数据，每一层对应一个观测时间步。目标是计算每个单元格在时间维度上的均值，但数据中存在大量NA，需要保证均值的稳健性——仅当单元格的有效观测数足够时才计算均值，否则将均值设为NA。比如示例场景中，要求当单元格时间序列里的NA数量≥2时，均值设为NA。

示例数据构建代码：

library(terra)
s <- rast(ncol=10, nrow=10, nlyr=30)
set.seed(1)
values(s) <- rnorm(size(s), 10)
s[3] <- NaN # 某个单元格全层设为NA
s[[4]] <- NaN # 某一层全设为NA
s[[1]][1] <- NaN # 单个单元格设为NA
s[[5]][4] <- NaN # 单个单元格设为NA

解决方案 #

核心思路是先统计每个单元格的NA数量，再基于阈值筛选出符合条件的单元格计算均值，不符合条件的直接设为NA。以下提供两种实现方式：

方式一：分步处理

1. 统计每个单元格的NA数量
   用app()函数遍历每个单元格的时间序列，统计NA个数：

na_count <- app(s, function(x) sum(is.na(x)))

2. 计算基础均值（忽略NA）
   先按常规方法计算均值：

mn <- mean(s, na.rm = TRUE)

3. 过滤NA超标的单元格
   根据设定的阈值（示例为NA≥2），将对应单元格的均值替换为NA：

mn[na_count >= 2] <- NA

方式二：自定义函数一步完成

直接在app()中使用自定义函数，同时完成NA统计和均值计算：

# 定义稳健均值函数：NA数量≥2时返回NA，否则返回有效观测的均值
robust_mean <- function(x) {
  na_num <- sum(is.na(x))
  if (na_num >= 2) {
    NA
  } else {
    mean(x, na.rm = TRUE)
  }
}

# 应用到Raster对象
mn <- app(s, robust_mean)

效果验证 #

• 第3个单元格全层都是NA（NA数量=30），均值会被设为NA
• 第4层全为NA，因此所有单元格的NA数量至少为1；那些额外还有NA的单元格（如第1层第1个、第5层第4个）NA数量≥2，均值也会被设为NA
• 仅NA数量<2的单元格会保留计算出的均值

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