基于terra::focalCpp自定义Rcpp函数计算窗口内特定值占比的问题

问题分析 #

原Rcpp函数存在多处逻辑错误，导致focalCpp调用失败：

• 错误地将整个输入向量的非NA值作为统计基数，而非当前窗口内的非NA值
• 用标量累加的方式完全偏离了原R函数的计数逻辑
• NA/NaN的判断逻辑混乱，导致统计结果完全错误

修正后的Rcpp函数 #

以下是适配terra::focalCpp的正确实现，完全对齐原R函数的功能：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector fxnpercent_cpp(NumericVector x, size_t ni, size_t nw) {
  NumericVector out(ni);
  
  size_t start = 0;
  for (size_t i = 0; i < ni; ++i) {
    size_t end = start + nw;
    int count_total = 0;
    int count_ones = 0;
    
    // 遍历当前窗口的所有元素
    for (size_t j = start; j < end; ++j) {
      // 跳过窗口掩码中的NA（对应不在缓冲区范围内的像素）
      if (is_na(x[j])) {
        continue;
      }
      count_total++;
      // 统计等于1的有效像素
      if (x[j] == 1.0) {
        count_ones++;
      }
    }
    
    // 处理窗口内全为NA的边界情况
    if (count_total == 0) {
      out[i] = NA_REAL;
    } else {
      // 计算占比并转为百分比
      out[i] = (static_cast<double>(count_ones) / count_total) * 100.0;
    }
    
    start = end;
  }
  
  return out;
}

使用说明 #

1. 编译上述函数后，直接替换原函数调用focalCpp：

percent_cpp = terra::focalCpp(r, w=w, fun=fxnpercent_cpp, na.policy="all")

2. 可通过以下代码验证结果与原terra::focal的一致性：

# 原方法生成结果
percent = terra::focal(x=r, w=w, fun=fxnpercent, na.policy="all")
# 对比两者结果
all.equal(values(percent), values(percent_cpp), na.rm=TRUE)

关键修正点 #

• 针对每个窗口单独统计非NA总数和1的数量，而非整个输入向量
• 正确跳过窗口掩码中的NA（即w矩阵中标记为NA的位置）
• 处理窗口内全为NA的边界情况，返回NA避免除零错误
• 使用static_cast进行类型转换，确保浮点计算的准确性

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