如果单元格包含一个点，则将栅格单元格值更改为 NA

Question

我在单个 stack 以及一个包含点坐标的 SpatialPoints 对象中包含多个栅格图层。如果单元格包含一个点，我正在尝试将栅格图层的值更改为 NA。

我在下面提供了一个可重现的示例。然而，为了把它放在上下文中，我的真实数据由 10 层“栖息地”（海拔、树冠覆盖等）组成，它们都包含在一个堆栈中。 rasters 覆盖的区域大约为 34 x 26 公里。此外，我还有大约 10,000 个来自动物的 GPS 位置。

因此，对于可重现的示例。

library(raster)
library(sp)

制作三个rasters并将它们组合成一个stack

set.seed(123)
r1 <- raster(nrows=10, ncols=10)
r1 <- setValues(r1, sample(c(1:50), 100, replace = T))

r2 <- raster(nrows=10, ncols=10)
r2 <- setValues(r2, sample(c(40:50), 100, replace = T))

r3 <- raster(nrows=10, ncols=10)
r3 <- setValues(r3, sample(c(50:55), 100, replace = T))

Stack <- stack(r1, r2, r3)
nlayers(Stack)

制作SpatialPoints

Pts <- SpatialPoints(data.frame(x = sample(-175:175, 50, replace = T),
                                y = sample(-75:75, 50, replace = T)))

绘制其中一个栅格和点

plot(Stack[[2]])
plot(Pts, add = T)

现在，是否可以将每个包含至少一个点的单元格的值更改为 NA？在stack 上执行此操作会很棒。

Answer 1

我会使用extract()，它可以被要求返回每个点所在的栅格像元的像元编号：

ii <- extract(Stack, Pts, cellnumbers=TRUE)[,"cells"]
Stack[ii] <- NA

## Check any one of the layers to see that this worked:
plot(Stack[[2]])
plot(Pts, add=TRUE)

或者，可以使用rasterize()，但对于非常大的栅格，它可能会（？）速度较慢：

ii <- !is.na(rasterize(Pts, Stack))
Stack[ii] <- NA