在第二层/属性上提取具有缓冲条件的栅格平均值

debugcn 发表于 Dev

哈利勒·泰伯

我有这个星星对象（也可以格式化为栅格）：

stars object with 2 dimensions and 2 attributes
attribute(s):
   LST_mean        elevation    
 Min.   :14.98   Min.   :296.0  
 1st Qu.:16.89   1st Qu.:346.9  
 Median :17.64   Median :389.3  
 Mean   :17.52   Mean   :389.2  
 3rd Qu.:18.18   3rd Qu.:428.3  
 Max.   :20.11   Max.   :521.6  
dimension(s):
  from to  offset    delta                       refsys point values    
x   71 83 4387654  860.241 DHDN / 3-degree Gauss-Kru... FALSE   NULL [x]
y   33 41 5598885 -860.241 DHDN / 3-degree Gauss-Kru... FALSE   NULL [y]

具有2个属性（在栅格情况下为图层）：温度和海拔。使用温度，我想选择落入缓冲区内的像素并返回均值，仅针对每次与所考虑的像素之间的高程差小于90米的像素。

任何想法如何做到这一点？计算落入缓冲区的像素的平均值非常容易，但是我找不到一种对它们设置任何条件的方法。

我将非常感谢您的帮助和建议。使用其他软件包的方法satrs也非常受欢迎:)

迈克尔·多曼

请使用下面的解决方案terra。该代码用于terra::extract创建两个对应list的：

像素值
周围的缓冲区值

随后，使用mapply与您建议的功能相似的功能，使用来成对处理这些值。

这是我第一次使用，terra但似乎terra::extract比快得多raster::extract，因此，即使对于较大的栅格，此解决方案也可能可行。

创建样本数据：

library(sf)
library(terra)

r = rast(ncol = ncol(volcano), nrow = nrow(volcano), xmin = 0, xmax = ncol(volcano), ymin = 0, ymax = nrow(volcano))
values(r) = volcano
s = r
s[] = rnorm(ncell(s))
r = c(r, s)
crs(r) = ""
plot(r)

计算缓冲区：

pnt = as.points(r, values = FALSE)
pol = buffer(pnt, 10)

从点提取栅格值：

x = extract(r, pnt)
head(x)

##      ID lyr.1       lyr.1
## [1,]  1   100 -0.03525223
## [2,]  2   100  0.31525467
## [3,]  3   101  0.94054608
## [4,]  4   101  0.37209238
## [5,]  5   101 -0.38388234
## [6,]  6   101 -0.03120593

从缓冲区提取栅格值：

y = extract(r, pol)
head(y)

##      ID lyr.1       lyr.1
## [1,]  1   100 -0.03525223
## [2,]  1   100  0.31525467
## [3,]  1   101  0.94054608
## [4,]  1   101  0.37209238
## [5,]  1   101 -0.38388234
## [6,]  1   101 -0.03120593

现在，可以使用依次处理提取的值mapply。首先，我们将对象转换为lists：

x = as.data.frame(x)
x = split(x, x$ID)
y = as.data.frame(y)
y = split(y, y$ID)

接下来，我们使用mapply每次考虑当前焦点值x和周围缓冲区值进行必要的计算y：

f = function(x, y) {
    d = abs(x[, 2] - y[, 2])  ## differences
    values = y[, 3]  ## values
    mean(values[d < 5], na.rm = TRUE)  ## Mean of subset
}
result = mapply(f, x, y)

最后，将结果放回到栅格模板中：