有没有更快的方法来使用多个多边形来屏蔽多个栅格?
Is there a faster way to mask multiple rasters using multiple polygons?
基本上我有 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 个。我需要做进一步的图像分析,所以我需要保留名字。我希望这是有道理的,谢谢!
这是一个使用 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])
}
或者做一个object/file,如果可行的话。
rr <- rast(inf)
mm <- mask(rr, v)
基本上我有 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 个。我需要做进一步的图像分析,所以我需要保留名字。我希望这是有道理的,谢谢!
这是一个使用 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])
}
或者做一个object/file,如果可行的话。
rr <- rast(inf)
mm <- mask(rr, v)