CODEWITHSUNDEEP
rloopscroprastermask
seak23
seak23

Reputation: 317

Is there a faster way to mask multiple rasters using multiple polygons?

Basically I have 12 multispectral images, and I want to mask them using 2 polygons (small waterbodies). The 2 polygons are in one shapefile, but I can break them up if it would make the process easier. With the help of some nice users on here, I tested this all out using the 12 images on one polygon and it works just fine, but I'll eventually need to do this for multiple polygons so I want to adapt my code.

The loop to crop all rasters using a single polygon:

#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]) 
}

Like I said this works just fine, but I want to mask instead of crop. Additionally, I need to mask using multiple polygons.

My working loop to do the same thing but for multiple (2) polygons:

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)
}

This took around 2 hours (which makes me suspicious) and I think it lost the polygon id somewhere along the way. When I tried plotting the results the plot was empty. My final output should be 24 rasterstacks, 12 for each polygon. I will need to do further image analysis so I will need to keep the names. I hope this makes sense and thank you!

Upvotes: 0

Views: 1093

Answers (1)

Robert Hijmans
Robert Hijmans

Reputation: 47061

Here is a minimal, self-contained, reproducible example using terra because it is much faster than raster (make sure you are using the current version)

Raster dataset with 12 layers

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)

Two "lakes"

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

Solution:

x <- mask(r, v)

And always try things for a single case before running the loop.

So if you have 12 files, you can do something like

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]) 
}

It might be a little faster to instead do this (only rasterize the polygons once)

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]) 
}

Or make one object/file, if that is practical.

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

Upvotes: 2

Related Questions