我试图通过使用focalCpp来加速我正在使用terra::focal进行的一些光栅处理。
下面是一些包含1和NA的示例数据,用于复制实际数据集
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)
字符串
这是我试图在Rcpp中复制的与focal一起使用的原始函数
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
}
型
这是原始聚焦函数,窗口近似为200m缓冲区
# 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")
型
这是我在Rcpp中复制fxnpercent函数的尝试
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;
}'
)
型
经过大量的故障排除,以获得正确的语法在Rcpp,我试图运行这个函数与focalCpp和我有一个错误。
percent = focalCpp(r, w=w, fun=fxnpercent, na.policy="all")
#my current error
Error: [focalCpp] test failed
型
我想我需要在Rcpp函数的窗口循环中做一些计算,但是我在理解如何设置它以正确工作方面遇到了麻烦。
1条答案
按热度按时间k10s72fa1#
NaN检查不是通过is_nan()完成的,而是使用std::isnan()
字符串