I'm trying to replicate the use of median within R that includes na.rm=TRUE as Rcpp code. I found this really useful link that includes the exact code I need for implementing Rccp median with na.rm=TRUE within my code here: https://github.com/RcppCore/Rcpp/issues/424. The median_dbl function from this link works great directly within R but implementing the function within a function I'm using for focalCpp has not worked for me.
Here's some example data
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)
Here is my moving window code
# 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
R function I want to replicate using Rcpp
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)
}
Here is my Rcpp code with using the median function from the link https://github.com/RcppCore/Rcpp/issues/424
I think there is an issue with having NAs within Rcpp code and running the focalCpp function. When running this function I have this error Error: [focalCpp] test failed . Which is not helpful, when looking at the source code. The function works if I have an if statement to ignore all of the nas if (!std::isnan(x[j])) but this does not create my desired output.
There seems to be some issue with having NAs while calculating by cells within the moving window
// [[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 function to run Rcpp function
output<- focalCpp(r, w=w, fun=testFunction, fillvalue = 0)
Starting simple helped to get to the solution and get past R crashing. The second loop and going from a double to a NumericVector seemed to be the issue. I'm using range instead for the window and that worked for me.
// [[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;
}