r - 如果单元格包含点，则将栅格单元格值更改为 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