问题描述
我有这个星星对象(也可以格式化为光栅):
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.31525467
## [3,] 1 101 0.94054608
## [4,] 1 101 0.37209238
## [5,] 1 101 -0.38388234
## [6,] 1 101 -0.03120593
现在,可以使用 mapply
顺序处理提取的值。
首先,我们将对象转换为 list
s:
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)