Cliping several rasters with a multi-polygon shapefile
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
The problem!
Imagine you have this situation: you have several global raster files and a shapefile with a few areas (e.g. Natural Parks). You want to generate a raster file using as mask each of the polygons in the shapefile for each of the original rasters. So, if you have 5 global rasters and a shapefile with 10 polygons, the output of this would be 50 rasters (a smaller raster for each polygon, cut from each of the 5 larger rasters).
Well I had this problem, as you might have guessed. I looked for solutions, maybe there are some, but I could not find any! I’m sure there are solutions out there… but let me show you mine!
This is a simple task, but if it is not automated, it’s a grueling one… As such, the natural step, for me, was using R to automate this work.
The solution!
First we need to load the raster package and the shapefile:
library(raster)
polygon_areas <- raster::shapefile("C:/yourshapefile.shp")
This is the code for the function I created, called crop_save:
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"))
}
}
}
The arguments for this function are:
origin_folder – Where the original rasters are saved.
pattern – This is a character string to identify raster files: in the folder were rasters are saved there are, generally, other files. This argument allows the selection of only rasters (e.g. tif files).
destination_folder – Folder where the otput folder will be created.
name_sub_folder – Name of the sub-folder to be created inside the destination folder. Inside this, a folder is created for each of the original rasters where the smaller rasters for each polygon are saved.
crop_areas – Areas to be used in the raster croping (a SpatialPolygonsDataFrame created by importing the shapefile into R).
name_crop_areas – Column of the SpatialPolygonsDataFrame with the unique names or codes for the regions.
An example (not run, you have to try this with your own rasters):
crop_save(origin_folder = "D:/THIS_FOLDER/"
, pattern = ".tif
, destination_folder = "C:/OUTPUT/"
, name_sub_folder = "Cut_rasters"
, crop_areas = polygon_areas
, name_crop_areas = "polygon_areas$Id"
)
I hope this is useful!
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.