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

huangapple
go评论122阅读模式
英文:

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

问题

以下是您要翻译的内容：

  1. 我正在尝试使用terra::focal来加速一些栅格处理，通过使用focalCpp来实现。
  3. 这里有一些示例数据，其中包含了1和NAs，以复制实际数据集
  5. ```R
  6. nr <- nc <- 50
  7. r <- rast(ncols=nc, nrows=nr, ext= c(0, nc, 0, nr))
  8. 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)
  9.

这是我尝试在Rcpp中复制的原始用于focal的函数：

  1. fxnpercent = function(x) {
  2.   q=x # 复制x
  3.   q[q!=1] = NA # q变为1
  4.   length(q[!is.na(q)])/length(x[!is.na(x)]) * 100 # 获取1的百分比
  5. }

这是原始focal函数，窗口通过200米缓冲区来近似：

  1. # 移动窗口矩阵
  2. mat = matrix(1,15,15) # 创建一个包围缓冲区范围的1矩阵
  3. gr = expand.grid(1:nrow(mat), 1:nrow(mat)) # 基于行/列索引的所有成对值的数据框
  4. center = 8 # 方形网格的中心索引
  5. gr$dist = sqrt((gr$Var1-center)^2 + (gr$Var2-center)^2) # 欧几里德距离计算
  6. threshold = 200/30 # 200米阈值转换为像素数（从中心计算）
  7. gr$inside = ifelse(gr$dist < threshold, 1, NA) # 如果距离小于阈值，则网格值为1，否则为NA
  8. w = matrix(gr$inside, 15,15)  # 使用gr$inside，将索引值放入原始维度的矩阵中
  10. # 输出移动窗口的百分比
  11. percent = terra::focal(x=r, w=w, fun=fxnpercent, na.policy="all")

这是我尝试在Rcpp中复制fxnpercent函数的代码：

  1. cppFunction( 
  2.     'NumericVector fxnpercent(NumericVector x, size_t ni, size_t nw) {
  3. 		NumericVector out(ni);
  5. 		// 循环处理单元格
  6. 		size_t start = 0;
  7. 		for (size_t i=0; i<ni; i++) {
  8. 			size_t end = start + nw;
  9. 			// 计算窗口内的值
  10. 			double v = 0;
  11. 			// 循环处理窗口内的值
  12. 			for (size_t j=start; j<end; j++) {
  13. 				if (x[j] != 1) {
  14. 				  v += std::nan("");
  15. 				}
  16. 				else {
  17. 				  v += x[j];
  18. 				}
  19. 				
  20. 			}
  21. 			
  22. 			NumericVector x1 = x[!is_na(x)];
  23. 			NumericVector v1 = v;
  24. 			NumericVector v2 = v1[!is_nan(x)];
  25. 			
  26. 			size_t v2size = v2.size();
  27. 			size_t x1size = x1.size();
  28. 			
  29. 			out[i] = (v2size / x1size) * 100;
  30. 			
  31. 			start = end;
  32. 		}
  33. 		return out;
  34. 	}'
  35.   )

在对Rcpp的语法进行了大量故障排除以确保正确性之后，我尝试使用focalCpp运行此函数，但出现错误。

  2. percent = focalCpp(r, w=w, fun=fxnpercent, na.policy="all")
  4. # 我当前的错误
  6. Error: [focalCpp] 测试失败

我认为我需要在Rcpp函数的窗口循环内进行一些计算，但我在如何设置它以正常工作方面遇到了问题。

英文:

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

  1. nr &lt;- nc &lt;- 50
  2. r &lt;- rast(ncols=nc, nrows=nr, ext= c(0, nc, 0, nr))
  3. values(r) &lt;- 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)
  4.

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

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

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

  1. # moving window matrix 
  2. mat = matrix(1,15,15) # create matrix of 1&#39;s that envelopes the extent of the buffer
  3. gr = expand.grid(1:nrow(mat), 1:nrow(mat)) # df of all pairwise values based on row/col index
  4. center = 8 # centroid index of the square grid
  5. gr$dist = sqrt((gr$Var1-center)^2 + (gr$Var2-center)^2) # euclidean distance calucation 
  6. threshold = 200/30 # 200m threshold is converted into number of pixels from center
  7. gr$inside = ifelse(gr$dist &lt; threshold, 1, NA) # if distance is less than threshold, grid value is one, otherwise NA
  8. w = matrix(gr$inside, 15,15)  # Using gr$inside, place indexed values into matrix of original dimensions
  10. #output percent from moving window
  11. percent = terra::focal(x=r, w=w, fun=fxnpercent, na.policy=&quot;all&quot;)

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

  1. cppFunction( 
  2.     &#39;NumericVector fxnpercent(NumericVector x, size_t ni, size_t nw) {
  3. 		NumericVector out(ni);
  5. 		// loop over cells
  6. 		size_t start = 0;
  7. 		for (size_t i=0; i&lt;ni; i++) {
  8. 			size_t end = start + nw;
  9. 			// compute something for a window
  10. 			double v = 0;
  11. 			// loop over the values of a window
  12. 			for (size_t j=start; j&lt;end; j++) {
  13. 				if (x[j] != 1) {
  14. 				  v += std::nan(&quot;&quot;);
  15. 				}
  16. 				else {
  17. 				  v += x[j];
  18. 				}
  19. 				
  20. 			}
  21. 			
  22. 			NumericVector x1 = x[!is_na(x)];
  23. 			NumericVector v1 = v;
  24. 			NumericVector v2 = v1[!is_nan(x)];
  25. 			
  26. 			size_t v2size = v2.size();
  27. 			size_t x1size = x1.size();
  28. 			
  29. 			out[i] = (v2size / x1size) * 100;
  30. 			
  31. 			start = end;
  32. 		}
  33. 		return out;
  34. 	}&#39;
  35.   )

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

  2. percent = focalCpp(r, w=w, fun=fxnpercent, na.policy=&quot;all&quot;)
  4. #my current error
  6. 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.

答案1

得分: 1

NaN检查不是通过is_nan()来完成的，而是使用std::isnan()

  1. #include <Rcpp.h>
  2. using namespace Rcpp;
  4. // [[Rcpp::export]]
  5. NumericVector fxnpercent(NumericVector x, size_t ni, size_t nw) {
  6.   NumericVector out(ni);
  7.   
  8.   // 循环遍历单元格
  9.   size_t start = 0;
  10.   for (size_t i = 0; i < ni; i++) {
  11.     size_t end = start + nw;
  12.     // 计算窗口内的某些内容
  13.     double v = 0;
  14.     double totalCount = 0;
  15.     double countOfOnes = 0;
  17.     // 循环遍历窗口内的值
  18.     for (size_t j = start; j < end; j++) {
  19.       if (!std::isnan(x[j])) {
  20.         totalCount++;
  21.         if (x[j] == 1) {
  22.           countOfOnes++;
  23.         }
  24.       }
  25.     }
  27.     if (totalCount > 0) {
  28.       out[i] = (countOfOnes / totalCount) * 100;
  29.     } else {
  30.       out[i] = std::nan(""); // 如果全部都是NaN，则为NaN
  31.     }
  32.     
  33.     start = end;
  34.   }
  35.   return out;
  36. }
英文:

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

  1. #include&lt;Rcpp.h&gt;
  2. using namespace Rcpp;
  4. // [[Rcpp::export]]
  5. NumericVector fxnpercent(NumericVector x, size_t ni, size_t nw) {
  6.   NumericVector out(ni);
  7.   
  8.   // loop over cells
  9.   size_t start = 0;
  10.   for (size_t i=0; i&lt;ni; i++) {
  11.     size_t end = start + nw;
  12.     // compute something for a window
  13.     double v = 0;
  14.     double totalCount = 0;
  15.     double countOfOnes = 0;
  17.     // loop over the values of a window
  18.     for (size_t j=start; j&lt;end; j++) {
  19.       if (!std::isnan(x[j])) {
  20.         totalCount++;
  21.         if (x[j] == 1) {
  22.           countOfOnes++;
  23.         }
  24.       }
  25.     }
  27.     if (totalCount &gt; 0) {
  28.       out[i] = (countOfOnes / totalCount) * 100;
  29.     } else {
  30.       out[i] = std::nan(&quot;&quot;); // it is NaN if all are NaN
  31.     }
  32.     
  33.     start = end;
  34.   }
  35.   return out;
  36. }

huangapple
  • 本文由 发表于 2023年7月28日 04:26:44
  • 转载请务必保留本文链接：https://go.coder-hub.com/76783194.html
  • c++
  • r
  • raster
  • rcpp
  • terra
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定