将single raster设置为NA，其中光栅堆栈的值为NA

Question

我有两个30m x 30m光栅文件，我想从中采样点。在采样之前，我想从图像中去除云层。我转向R和Hijman的光栅包来完成这项任务。

我使用绘图Poly（sp=TRUE）命令，绘制了18个不同的多边形。该函数似乎不允许18个多边形作为一个sp对象，所以我将它们全部分开绘制。然后，我给多边形一个与光栅匹配的投影4string，并将它们设置为列表。我通过lApplication函数运行列表以将它们转换为光栅（Hijman软件包中的光栅化函数），多边形区域设置为NA，图像的其余部分设置为1。

我的最终目标是一个光栅层，18个区域设置为NA。我尝试堆叠光栅化多边形列表，并将其子集设置为在相同区域中设置一个新的光栅到NA。我的可重复代码在下面。

library(raster)
r1 <- raster(nrow=50, ncol = 50)
r1[] <- 1
r1[4:10,] <- NA
r2 <- raster(nrow=50, ncol = 50)
r2[] <- 1
r2[9:15,] <- NA
r3 <- raster(nrow=50, ncol = 50)
r3[] <- 1
r3[24:39,] <- NA

r4 <- raster(nrow=50, ncol = 50)
r4[] <- 1

s <- stack(r1, r2, r3)

test.a.cool <- calc(s, function(x){r4[is.na(x)==1] <- NA})

无论出于什么原因，该死的测试池是一个空白图，我的目标是将其作为一个光栅，除了堆栈中的NAs，s等于1之外的所有值。

有什么提示吗？

谢谢。

Answer 1

我不熟悉您使用的软件包，但是看看代码中的最后一行，我认为问题可能在这里：

 function(x){r4[is.na(x)==1] <- NA})

看起来calc对此没有太大帮助。它正在设置由x的NAs索引的r4的值并将其设置为NA。

然后呢？如果有的话，也许：

 function(x){r4[is.na(x)==1] <- NA; return(r4) })

虽然，还不清楚这是否是你想要的。

Answer 2

我也在R-Sig-Geo论坛上发布了这个问题，并收到了软件包作者的回复。两种最简单的解决方案：

使用sp包将我的多边形绑定为一个，然后光栅化多边形。

p <- rbind(p1, p2, p3...etc., makeUniqueIDs = TRUE)

r4 <- raster(nrow=50, ncol = 50)
r4[] <- 1
mask <- rasterize(p, r4)
mask[mask %in% 1:18] <- 1
#The above code produces a single raster file with 
#my polygons as unique values, ready for masking.

第二个简单的解决方案，正如Josh O'Brien所指出的：

m <- sum(s)
test <- mask(r4, m)

R社区触礁。问题在一小时内解决了（两次）。谢谢

Answer 3

执行求和将起作用，因为对于堆栈中有一个值为NA的任何网格单元，求和将返回NA。

要查看其是否有效，请比较以下数据：

plot(s)
plot(sum(s))