R – How to Clip Several Rasters with Multipolygon and Save as Separate Rasters Using R

cliprrastershapefile

I'm trying to replicate this [example][1] but I'm not getting separate raster as an output. In my example, I have a multipolygon which has 2 polygons and 2 raster. My output I want to be 4 raster layers and the name of the raster I'd like to be 'id_band' where the 'id' part will come from the polygon's id column and the 'band' will come from the original raster (i.e. the final should be 14985_b6.tif).

Here is the code I'm using:

 library(raster)
polygon_areas <- raster::shapefile("test.shp")
plot(polygon_areas)
crop_save <-
  function(origin_folder,
           pattern,
           destination_folder,
           name_sub_folder,
           crop_areas,
           name_crop_areas) {
    file_list <- list.files(path = origin_folder, pattern)
    #Create folder
    dir.create(paste0(destination_folder, "/", name_sub_folder))
    how_many_areas <- nrow(crop_areas)
    #Create raster stack
    raster_stack <- stack()
    #File paths
    paths1 <- paste0(origin_folder, file_list)
    #Load rasters to stack
    for (i in 1:length(file_list)) {
      raster_stack <- stack(raster_stack, raster(paths1[i]))
    }
    names_list <-  eval(parse(text = name_crop_areas))
    numbers <- 1:length(names_list)
    names_list <-
      paste0(as.character(numbers), "_polygon_", names_list)
    polyRR_list <- list()
    for (x in 1:nrow(crop_areas)) {
      pol1 <- assign(names_list[x], crop_areas[x, ])
      polyRR_list[[x]] <- pol1
    }
    for (j in 1:nlayers(raster_stack)) {
      dir.create(paste0(
        destination_folder,
        "/",
        name_sub_folder,
        "/",
        names(raster_stack)[j]
      ))
      for (k in 1:length(polyRR_list)) {
        a <- crop(raster_stack[[j]], polyRR_list[[k]])
        a <-
          mask(
            a,
            polyRR_list[[k]],
            filename = paste0(
              destination_folder,
              "/",
              name_sub_folder,
              "/",
              names(raster_stack)[j],
              "/",
              "RR",
              polyRR_list[[k]]$Id,
              ".tif"
            )
          )
      }
    }
  }

crop_save(
  origin_folder = "C:/Users/Geography/Desktop/a/" #Where your rasters are
  ,
  pattern = ".tif$"
  ,
  destination_folder = "C:/Users/Geography/Desktop/a/output/" #folder to save the rasters resulting from this operation
  ,
  name_sub_folder = "any_name" # Name of the sub-folder to be created inside the destination folder
  ,
  crop_areas = polygon_areas # The shapefile
  ,
  name_crop_areas = "polygon_areas$id" # The unique Id of the polygons in the shapefile
)

Also, the data I'm using:

list(new("SpatialPolygonsDataFrame", data = structure(list(id = c(14985, 
23552), left = c(2665569.79702, 2670699.79702), top = c(1412205.45671, 
1407615.45671), right = c(2665839.79702, 2670969.79702), bottom = c(1411935.45671, 
1407345.45671), id_sample = c(1L, 1L)), class = "data.frame", row.names = 0:1), 
    polygons = list(new("Polygons", Polygons = list(new("Polygon", 
        labpt = c(80.3367355387047, 16.2821390031233), area = 6.1440311016086e-06, 
        hole = FALSE, ringDir = 1L, coords = structure(c(80.3354759967946, 
        80.3379981575006, 80.3379950690668, 80.3354729314431, 
        80.3354759967946, 16.2833585003835, 16.2833555285503, 
        16.2809195021546, 16.2809224739735, 16.2833585003835), .Dim = c(5L, 
        2L)))), plotOrder = 1L, labpt = c(80.3367355387047, 16.2821390031233
    ), ID = "0", area = 6.1440311016086e-06), new("Polygons", 
        Polygons = list(new("Polygon", labpt = c(80.3523654989901, 
        16.2402853833502), area = 6.1434245505222e-06, hole = FALSE, 
            ringDir = 1L, coords = structure(c(80.3511062267076, 
            80.3536279907441, 80.353624759728, 80.3511030187666, 
            80.3511062267076, 16.2415050209594, 16.2415019107378, 
            16.2390657420334, 16.2390688522397, 16.2415050209594
            ), .Dim = c(5L, 2L)))), plotOrder = 1L, labpt = c(80.3523654989901, 
        16.2402853833502), ID = "1", area = 6.1434245505222e-06)), 
    plotOrder = 1:2, bbox = structure(c(80.3354729314431, 16.2390657420334, 
    80.3536279907441, 16.2833585003835), .Dim = c(2L, 2L), .Dimnames = list(
        c("x", "y"), c("min", "max"))), proj4string = new("CRS", 
        projargs = "+proj=longlat +datum=WGS84 +no_defs")), "b6.tif", 
    "b6.tif.aux.xml", "b6.tif.ovr", "b6.tif.vat.cpg", "b6.tif.vat.dbf", 
    "b6.tif.xml", "b7.tif", "b7.tif.aux.xml", "b7.tif.vat.cpg", 
    "b7.tif.vat.dbf")

I'm using R.4.1.1, RStudio 1.4.1717 and Windows10 x64
[1]: https://www.r-bloggers.com/2021/02/cliping-several-rasters-with-a-multi-polygon-shapefile/

Best Answer

I managed to solve the problem but instead of a (single) multipolygon I used separate shapefiles. Below is the code:

library(raster)
library(sf)

file_list <- list.files("mydir", pattern = "*shp", full.names = TRUE)
shapefile_list <- lapply(file_list, read_sf)
polygons <- do.call(rbind, shapefile_list)


crop_save <- function(origin_folder, pattern, destination_folder, name_sub_folder, crop_areas, name_crop_areas){
  file_list <- list.files(path = origin_folder, pattern)
  #Create folder
  dir.create(paste0(destination_folder,"/",name_sub_folder))
  how_many_areas <- nrow(crop_areas)
  #Create raster stack
  raster_stack <- stack() 
  #File paths
  paths1 <- paste0(origin_folder,file_list)
  #Load rasters to stack
  for(i in 1:length(file_list)){
    raster_stack <- stack(raster_stack, raster(paths1[i]))  
  }  
  names_list <-  eval(parse(text=name_crop_areas))
  numbers <- 1:length(names_list)
  names_list <- paste0(as.character(numbers),"_polygon_", names_list)
  polyRR_list <- list()
  for(x in 1: nrow(crop_areas)){
    pol1 <- assign(names_list[x],crop_areas[x,])
    polyRR_list[[x]] <- pol1
  }
  #dir.create(paste0(destination_folder,"/",name_sub_folder, "/", raster_folder))
  for(j in 1:nlayers(raster_stack)){
    for(k in 1:length(polyRR_list)){
      a<-crop(raster_stack[[j]], polyRR_list[[k]])
      a<-mask(a,polyRR_list[[k]], filename = paste0(destination_folder,"/",name_sub_folder, "/", polyRR_list[[k]]$id,"_",names(raster_stack)[j], ".tif"))
    }
  }
}


crop_save(origin_folder = "mydir/"
          , pattern = "tif$"
          , destination_folder = "mydir" 
          , name_sub_folder = "crop_Rasters"
          , crop_areas = polygons 
          , name_crop_areas = "polygons$id" 
)
Related Question