开发一个自定义 Rcpp 函数，与 terra::focalCpp 一起使用来计算移动窗口内值的中位数

Question

我正在尝试复制 R 中中位数的使用，其中包括 na.rm=TRUE 作为 Rcpp 代码。我发现这个非常有用的链接，其中包括我在代码中使用 na.rm=TRUE 实现 Rccp 中位数所需的确切代码：https://github.com/RcppCore/Rcpp/issues/424。此链接中的median_dbl函数直接在R中工作得很好，但在我用于focalCpp的函数中实现该函数对我来说不起作用。

这是一些示例数据

nr <- nc <- 50
r <- rast(ncols=nc, nrows=nr, ext= c(0, nc, 0, nr))
values(r) <- rep(c(rep(NA, 10), rep(1, 10), seq(1,8), rep(-999,12), rep(NA,5), rep(1,15),seq(1,8), rep(NA,12), seq(1,20)), 50)
 

这是我的移动窗口代码

# moving window matrix 
mat = matrix(1,15,15) # create matrix of 1's that envelopes the extent of the buffer
gr = expand.grid(1:nrow(mat), 1:nrow(mat)) # df of all pairwise values based on row/col index
center = 8 # centroid index of the square grid
gr$dist = sqrt((gr$Var1-center)^2 + (gr$Var2-center)^2) # euclidean distance calucation 
threshold = 200/30 # 200m threshold is converted into number of pixels from center
gr$inside = ifelse(gr$dist < threshold, 1, NA) # if distance is less than threshold, grid value is one, otherwise NA
w = matrix(gr$inside, 15,15)  # Using gr$inside, place indexed values into matrix of original dimensions

我想使用 Rcpp 复制 R 函数

  testFunction = function(x) {
    q=x # make a copy of window matrix
    test.ind = which(q==-999) # find dummy indicies
    q[test.ind] = NA
    median(q, na.rm=T)
  }

这是我的 Rcpp 代码，使用链接中的中值函数https://github.com/RcppCore/Rcpp/issues/424

我认为在 Rcpp 代码中包含 NA 并运行 focusCpp 函数存在问题。运行此函数时出现此错误

Error: [focalCpp] test failed 

。当查看源代码时，这没有帮助。如果我有一个 if 语句来忽略所有 nas

if (!std::isnan(x[j]))

，则该函数可以工作，但这不会创建我想要的输出。

在移动窗口内的单元格计算时使用 NA 似乎存在一些问题

// [[Rcpp::export]]
Rcpp::NumericVector testFunction (Rcpp::NumericVector x, size_t ni, size_t nw) {
  Rcpp::NumericVector out(ni);
    // loop over cells
    size_t start = 0;
    for (size_t i=0; i<ni; i++) {
        size_t end = start + nw;
        // compute something for a window
        Rcpp::NumericVector v;
        // loop over the values of a window
        for (size_t j=start; j<end; j++) {
            double q = x[j];
            if (q == -999) {q = NA_REAL;}
            v[j] = q;
        }
        Rcpp::NumericVector v2 = na_omit(v);
        int n = v2.length();
        if(n == 0){
            out[i] = 0;
        }else{
          std::sort(v2.begin(), v2.end());
          if (n % 2 == 0) {
            out[i] = (v2[(n / 2) - 1] + v2[(n / 2)]) / 2;
          } else {
            out[i] = v2[(n / 2)];
          }
        }
        start = end;
        }
        return out;
    }

focalCpp函数运行Rcpp函数

  output<- focalCpp(r, w=w, fun=testFunction, fillvalue = 0)

Answer 1

从简单开始有助于找到解决方案并克服 R 崩溃问题。第二个循环以及从 double 到 NumericVector 似乎是问题所在。我使用范围代替窗口，这对我有用。

// [[Rcpp::export]]
Rcpp::NumericVector fmediantempC (Rcpp::NumericVector x, size_t ni, size_t nw) {
    Rcpp::NumericVector out(ni);
    // loop over cells
    for (size_t i=0; i<ni; i++) {
        size_t start = i*nw;
        size_t end = start+nw-1;
        Rcpp::NumericVector zw = x[Rcpp::Range(start,end)]; //Current window
        Rcpp::NumericVector v2 = Rcpp::na_omit(zw);
        int n = v2.length();
        if(n == 0){
            out[i] = NA_REAL; //if all NA values then marked as NA
        }else{
          std::sort(v2.begin(), v2.end());
          if (n % 2 == 0) {
            out[i] = (v2[(n / 2) - 1] + v2[(n / 2)]) / 2;
          } else {
            out[i] = v2[(n / 2)];
          }
    
        }
    }
    return out;
}