Clipping raster using shapefile in R, but keeping the geometry of the shapefile

Question

I am using {raster} to clip (or crop) a raster based on an irregular shapefile (the Amazon biome) but the output always has a rectangular extent. However, I need the output in the exact same geometry of the shapefile. Any tips? Cheers.

library(raster)
library(rgdal)

myshp <- readOGR("Amazon.shp", layer="Amazon")
e <- extent(myshp)
myraster <- raster("Temperature.tif")
myraster.crop <- crop(myraster, e, snap="out", filename="myoutput.tif")

Sterrett

Answer 1

One option is to use raster::mask()

library(maptools)  ## For wrld_simpl
library(raster)

## Example SpatialPolygonsDataFrame
data(wrld_simpl)
SPDF <- subset(wrld_simpl, NAME=="Brazil")

## Example RasterLayer
r <- raster(nrow=1e3, ncol=1e3, crs=proj4string(SPDF))
r[] <- 1:length(r)

## crop and mask
r2 <- crop(r, extent(SPDF))
r3 <- mask(r2, SPDF)

## Check that it worked
plot(r3)
plot(SPDF, add=TRUE, lwd=2)

Visakhapatnam

You need to mask first and then crop it to remove NA's not the other way around. – Chandrachandragupta

this is an incorrect solution... plot(r3, colNA="magenta") clearly shows that the raster has not been cropped but simply "masked" – Mameluke

@Mameluke I believe you are mistaken. Here is a nice explanation of the role of NAs in R rasters, submitted in response to an OP who had a similar concern to yours – Leger

Answer 2

Package terra:: made it simpler, you can ::crop and mask in the same step.

# Load packages
library(maptools) # For geometry
library(terra) # Perform the crop and mask

###--- Preparing polygon and raster ---###

# Example SpatialPolygonsDataFrame
data(wrld_simpl)
polygon <- subset(wrld_simpl, NAME=="Luxembourg")
plot(polygon) # have a look

# Convert from SpatialPolygonsDataFrame to SpatVector (terra package format)
# And create a smaller polygon with buffer (negative to be "inside")
polygon_bf <- buffer(vect(polygon), width= -100000)
plot(polygon_bf, add= T) # have a look on both

# Create a SpatRaster from a file
f <- system.file("ex/elev.tif", package="terra")
r_lux <- rast(f)
plot(r_lux) # have a look on SpatRaster (terra package format)

# See the steps with plot
plot(polygon, add= T)
plot(polygon_bf, add= T)

Click to see intermediate steps (all files)

########################################################
### Crop and mask by any polygon                     ###
raster_cp <- crop(r_lux, polygon_bf, mask= T)
# Note: if mask= F, the crop will be by extent (box) ###
########################################################

### Check the results
plot(raster_cp)
plot(polygon_bf, lwd=1, add=T)

Click to see the final output

In addition, (with raster package) in case you want to perform with a simple geometry (i.g. box), the coordinates of the extent can be place directly:

e <- as(extent(c(xmin= -16, xmax= -7.25, ymin= 4, ymax= 12.75)), 'SpatialPolygons')
crs(e) <- "+proj=longlat +datum=WGS84 +no_defs"
r <- crop(my_raster, e)

Intro