Develop a custom Rcpp function to be used with terra::focalCpp to calculate the percent of a specific value within a moving window

Question

I'm trying to speed up some raster processing I'm doing using terra::focal by using focalCpp.

Here is some example data with 1s and NAs included to replicate an actual dataset

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(1,12), rep(NA,5), rep(1,15),seq(1,8), rep(NA,12), seq(1,20)), 50)
 

This is the original function used with focal that I'm trying to duplicate in Rcpp

fxnpercent = function(x) {
  q=x # make a copy of x
  q[q!=1] = NA # q becomes just 1s
  length(q[!is.na(q)])/length(x[!is.na(x)]) * 100 # gets percentage of 1s
}

This is the original focal function with a window approximated by a 200m buffer

# 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

#output percent from moving window
percent = terra::focal(x=r, w=w, fun=fxnpercent, na.policy="all")

And this is my attempt to duplicate the fxnpercent function in Rcpp

cppFunction( 
    'NumericVector fxnpercent(NumericVector x, size_t ni, size_t nw) {
        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
            double v = 0;
            // loop over the values of a window
            for (size_t j=start; j<end; j++) {
                if (x[j] != 1) {
                  v += std::nan("");
                }
                else {
                  v += x[j];
                }
                
            }
            
            NumericVector x1 = x[!is_na(x)];
            NumericVector v1 = v;
            NumericVector v2 = v1[!is_nan(x)];
            
            size_t v2size = v2.size();
            size_t x1size = x1.size();
            
            out[i] = (v2size / x1size) * 100;
            
            start = end;
        }
        return out;
    }'
  )

After lots of troubleshooting to get the syntax correct in Rcpp, I try to run this function with focalCpp and I have an error.

percent = focalCpp(r, w=w, fun=fxnpercent, na.policy="all")

#my current error

Error: [focalCpp] test failed

I think I need to do some calculations within the window loop in the Rcpp function but I'm having trouble with understanding how to set it up to work correctly.

Answer 1

NaN checking is not done via is_nan() but rather using std::isnan()

#include<Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector fxnpercent(NumericVector x, size_t ni, size_t nw) {
  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
    double v = 0;
    double totalCount = 0;
    double countOfOnes = 0;

    // loop over the values of a window
    for (size_t j=start; j<end; j++) {
      if (!std::isnan(x[j])) {
        totalCount++;
        if (x[j] == 1) {
          countOfOnes++;
        }
      }
    }

    if (totalCount > 0) {
      out[i] = (countOfOnes / totalCount) * 100;
    } else {
      out[i] = std::nan(""); // it is NaN if all are NaN
    }
    
    start = end;
  }
  return out;
}