Mask and crop a raster to multiple polygons in R

Question

Is it possible to mask and/or crop a raster to multiple polygons, with the output being multiple raster's (i.e., a raster stack)?

Here's a reproducible example:

library(tidyverse)
library(raster)
library(sf)
filename <- system.file("external/test.grd", package="raster")
r <- raster(filename)

#  Create a data frame containing the centroids of a polygon 
(df <- data.frame(
  x = c(180000, 180500),
  y = c(331700, 332200)))

# create a point shapefile from the dataframe, then convert to polygon
polygon <- df %>%
  st_as_sf(coords = c("x", "y"), crs = crs(r)) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POINT") %>% 
  st_buffer(dist=200)

# plot the raster and polygon 
plot(r)
plot(polygon, add=TRUE)

# mask and crop the raster
        m <- r %>% mask(polygon) %>%
          crop(polygon)

The above code returns a single raster. How to return multiple rasters for each polygon?

Answer 1

You can do this with an lapply function for the iterator, in leu of a for loop. The result is a list object containing the feature subset rasters.

library(sf)
library(raster)

r <- raster(system.file("external/test.grd", package="raster"))
p <- as(sf::st_buffer(as(sampleRandom(r, 4, sp=TRUE), "sf"), 200), "Spatial")

( rp <- lapply(1:nrow(p), function(x) 
          raster::mask(crop(r, extent(p[x,])), p[x,]) ) )

If you subset (crop) to the common polygon extent you can create a stack.

( rp <- do.call(stack,lapply(1:nrow(p), function(x) 
          raster::mask(crop(r, extent(p)), p[x,]) )) )
plot(rp)

Answer 2

I would just loop through your data frame and collect the outputs in a new brick. You can't crop the outputs if you want them to be in the same brick/stack because the extents of the outputs needs to be the same.

library(tidyverse)
library(raster)
library(sf)

filename <- system.file("external/test.grd", package="raster")
r <- raster(filename)

df <- data.frame(
  x = c(180000, 180500),
  y = c(331700, 332200))

#Make a copy of the input raster to collect the outputs
output <- brick(r)

for(i in 1:nrow(df)){
  polygon <- df[i,] %>%
    st_as_sf(coords = c("x", "y"), crs = crs(r)) %>%
    summarise(geometry = st_combine(geometry)) %>%
    st_cast("POINT") %>% 
    st_buffer(dist=200)
  
  #You can't crop here because the extents of the outputs need to be the same
  m <- r %>% mask(polygon)
 
  #Add the new raster to the output brick
  output <- addLayer(output, m)
   
}

#Drop the first layer, which is just in the input raster
output <- dropLayer(output, 1)

If you want the output rasters to be cropped to the individual polygons, you could collect the outputs in a list instead of a brick/stack.

I hope I understood your question correctly.