Mask a list of rasters using a list of spatial vectors (shapefiles)

I have a list of raster files and a list of shapefiles ordered by date (one for each of 52 weeks of the year). I need to mask each raster with its corresponding shapefile and I am stuck on the code to make it happen. I've completed other data manipulations by making use of lapply, but since I need to iterate through both the rasters and the shapefiles, I'm stuck.

library(terra)
#> terra 1.7.78
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

# take a raster 
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
plot(r)


# take a shapefile
pathshp <- system.file("ex/lux.shp", package = "terra")
v <- vect(pathshp)
v2 <- subset(v, v$ID_1 ==1) # make a Spatial Vector to use for masking
plot(v2)


# mask the raster using the shapefile
masked_raster <- mask(r, v2, inverse=TRUE)
plot(masked_raster)


# repeat this masking for a list of rasters usign a list of spatial vectors

# make some more rasters
 s <- rast(system.file("ex/elev.tif", package="terra"))   
 x <- rep(s, 4)
 
# create a function to mask rasters
 mask_raster <- function(raster, shapefile) {
   masked_raster <- mask(raster, shapefile, inverse = TRUE)
   return(masked_raster)
 }

# use mappply and the new function to cycle through both rasters and Spatvectors
 masked_rasters <- mapply(mask_raster, x, v, SIMPLIFY = FALSE)
#> Error: unable to find an inherited method for function 'mask' for signature 'x = "SpatRaster", mask = "data.frame"'

I'm not sure whether mapply is the best way to go about solving this problem or if there is better way. When I run the code on a test of my actual data with just one raster and one shapefile, I get the following error as well:

Error: [vect] the variable name(s) in argument `geom` are not in `x`

TIA for your help!

Solution

mapply() suits just fine, but with mapply(mask_raster, x, v, SIMPLIFY = FALSE) , v being a SpatVector, you are calling mask_raster() with n-th layer of x and n-th column of v (a single-column df). You can add some verbosity to mask_raster() to check:

mask_raster <- function(raster, shapefile) {
  message(str(raster))
  message(str(shapefile))
  masked_raster <- mask(raster, shapefile, inverse = TRUE)
  return(masked_raster)
}

masked_rasters <- mapply(mask_raster, x, v , SIMPLIFY = FALSE)
#> S4 class 'SpatRaster' [package "terra"]
#> 
#> 'data.frame':    12 obs. of  1 variable:
#>  $ ID_1: num  1 1 1 1 1 2 2 2 3 3 ...
#> 
#> Error: unable to find an inherited method for function 'mask' for signature 'x = "SpatRaster", mask = "data.frame"'

With this particular dataset you could split v by ID_1 or NAME_1 first to get a list. And adjust number of raster layers to match the lenght of v

library(terra)
#> terra 1.7.78

pathshp <- system.file("ex/lux.shp", package = "terra")

# split by ID_1 attribute
v <- vect(pathshp) |> split("ID_1")
str(v)
#> List of 3
#>  $ Diekirch    :S4 class 'SpatVector' [package "terra"]
#>  $ Grevenmacher:S4 class 'SpatVector' [package "terra"]
#>  $ Luxembourg  :S4 class 'SpatVector' [package "terra"]

s <- rast(system.file("ex/elev.tif", package="terra"))   
# adjust number of layers to 3 to match v lenght
x <- rep(s, 3)
# 3 layers, layer can be accessed as a list items, e.g. x[[1]]
x
#> class       : SpatRaster 
#> dimensions  : 90, 95, 3  (nrow, ncol, nlyr)
#> resolution  : 0.008333333, 0.008333333  (x, y)
#> extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326) 
#> source      : elev.tif 
#> names       : elevation, elevation, elevation 
#> min values  :       141,       141,       141 
#> max values  :       547,       547,       547

mask_raster <- function(raster, shapefile) {
  masked_raster <- mask(raster, shapefile, inverse = TRUE)
  return(masked_raster)
}

masked_rasters <- mapply(mask_raster, x, v, SIMPLIFY = FALSE)

par(mfrow=c(1,length(masked_rasters)))
for (i in seq_along(masked_rasters)){
  plot(masked_rasters[[i]])
}

^{Created on 2024-06-20 with reprex v2.1.0}