如果栅格单元格包含一个点，则将栅格单元格的值更改为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上执行此操作。

- B. Davis

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

- Josh O'Brien

---

很棒。非常有帮助。你在Stack中为什么索引了层[[1]]？测试一下，似乎你可以一次对整个堆栈执行操作。例如：ii <- unique(extract(Stack, Pts, cellnumbers=TRUE)[,"cells"]) Stack[ii] <- NA - B. Davis

@B.Davis 没有什么特别的原因，实际上 unique() 的调用也是不必要的。我已经编辑了帖子以进行相应的简化。 - Josh O'Brien