3

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

使用 drawPoly(sp=TRUE) 命令，我绘制了 18 个不同的多边形。该函数似乎不允许将 18 个多边形作为一个 sp 对象，因此我将它们全部单独绘制。然后我给多边形一个与栅格匹配的 proj4string，并将它们设置为一个列表。我通过 lapply 函数运行列表，将它们转换为栅格（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})

无论出于何种原因，该死的 testacool 是一个空白图，我的目标是将它作为一个栅格，其中除了堆栈中的 NA s 之外的所有值都等于 1。

有小费吗？

谢谢。

4

4 回答 4

3

这样做sum(s)会起作用，因为堆栈中只有一个值的任何网格单元的sum()返回。NANA

要查看它是否有效，请比较以下生成的数字：

plot(s)
plot(sum(s))
于 2013-04-29T05:10:58.673 回答
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 社区摇滚。问题在一小时内解决（两次）。谢谢。

于 2013-04-29T05:14:42.460 回答
1

你走在正确的轨道上。该[运算符是为栅格和栅格堆栈定义的，因此您可以只使用一行：

r4[ any(is.na(s) ) ] <- NA
plot(r4)

如果你想使用calc，你可以像这样使用它：

r4 <- calc( s, function(x){ ( ! any( is.na(x) ) ) } )
r4[is.na(r4)] <- NA
plot(r4)
于 2013-04-29T06:17:44.943 回答
1

我不熟悉您使用的包，但是查看代码中的最后一行，我认为问题可能出在此处：

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

它看起来calc不会做太多事情。它正在设置由sr4索引的值并将其设置为。 NAxNA

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

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

虽然，尚不清楚这是否是你所追求的。

于 2013-04-29T04:16:54.773 回答