有没有更快的方法来使用多个多边形来掩盖多个栅格？

Posted 2023-02-25

【中文标题】有没有更快的方法来使用多个多边形来掩盖多个栅格？

【英文标题】：Is there a faster way to mask multiple rasters using multiple polygons?

【发布时间】：2021-05-04 00:18:02

【问题描述】：

基本上我有 12 个多光谱图像，我想使用 2 个多边形（小水体）对它们进行遮罩。这 2 个多边形在一个 shapefile 中，但如果可以使过程更容易，我可以将它们分解。在一些优秀用户的帮助下，我在一个多边形上使用 12 个图像对这一切进行了测试，效果很好，但我最终需要对多个多边形执行此操作，所以我想调整我的代码。

使用单个多边形裁剪所有栅格的循环：

#The single polygon 
mask <- st_read(here::here("data", "mask.shp") %>%
  st_as_sf()

#Creates list of input files and their paths
crop_in <- list.files(here::here("data", "s2_rasters"), pattern="tif$", full.names=TRUE)

#Creates list of output files and their directory.
crop_out <- gsub(here::here("data", "s2_rasters"), here::here("data", "s2_cropped"), crop_in)

for (i in seq_along(crop_in)) 
  b <- brick(crop_in[i])
  crop(b, mask, filename = crop_out[i]) 

就像我说的那样，这很好用，但我想遮罩而不是裁剪。此外，我需要使用多个多边形进行遮罩。

我的工作循环做同样的事情，但多 (2) 个多边形：

masks_2 <- st_read(here::here("data", "multiple_masks.shp")) %>% 
  st_as_sf()

for (i in seq_along(crop_in)) 
  
  b <- brick(crop_in[i])
  mask(b, masks_2, filename = crop_out[i], overwrite = TRUE)

这花了大约 2 个小时（这让我很怀疑），我认为它在途中的某个地方丢失了多边形 ID。当我尝试绘制结果时，情节是空的。我的最终输出应该是 24 个栅格堆栈，每个多边形 12 个。我需要做进一步的图像分析，所以我需要保留名字。我希望这是有道理的，谢谢！

【参考方案1】：

这是一个使用terra 的最小、独立、可重现的示例，因为它比raster 快得多（确保您使用的是当前版本）

12 层的栅格数据集

library(terra)
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
r <- rep(r, 12) * 1:12
names(r) <- paste0("band", 1:12)

两个“湖”

v <- vect(system.file("ex/lux.shp", package="terra"))
v <- v[c(1,12)] 

解决方案：

x <- mask(r, v)   

并且总是在运行循环之前尝试单个案例。

所以如果你有 12 个文件，你可以做类似的事情

inf <- list.files("data/s2_rasters", pattern="tif$", full.names=TRUE)
outf <- gsub(".tif$", "_masked.tif", inf)

for (for i in 1:length(inf)) 
    r <- rast(inf[i])
    m <- mask(r, v, filename=outf[i]) 

这样做可能会快一点（只对多边形进行一次光栅化）

msk <- rast(inf[1])
msk <- rasterize(v, msk)
for (for i in 1:length(inf)) 
    r <- rast(inf[i])
    m <- mask(r, msk, filename=outf[i]) 

如果可行，也可以制作一个对象/文件。

rr <- rast(inf)
mm <- mask(rr, v)