2

在此处输入图像描述 在此处输入图像描述

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

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的方法也非常受欢迎:)

4

1 回答 1

2

请参阅下面的解决方案，使用terra. 代码terra::extract用来创建两个对应list的s：

  • 像素值
  • 周围的缓冲区值

随后，使用 , 对值进行成对处理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)

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

u = r[[1]]
values(u) = result
plot(u)

输出

于 2021-01-20T08:15:47.300 回答