Search code examples
rterra

Crop raster with shapefile that has different projections

I tried to terra::crop() a SpatRaster (e.g. from the ForestClim database downloaded from here [6GB download]) with a shapefile, to keep only data for Belgium:

# shapefile
ISO <- countrycode("Belgium", origin = 'country.name', destination = 'iso3c')
shp <- gadm(country=ISO, level=0, 
            path = tempdir(), 
            version="latest")
crs(shp) <- "epsg:4326"

# raster
r <- rast(~/Downloads/ForestClim_01.tif) # download link above
crs(r) <- "epsg:4326"

# crop
cr <- crop(r, shp, mask=TRUE, overwrite=TRUE)
Error: [crop] extents do not overlap

I have tried to understand this error but without success. I don't understand why this doesn't work as it worked perfectly with my CHELSA data.

Solution

  • As noted in the previous answer, both files have difference CRS. In these cases, it is often easier to project an sf/SpatVector object to match the target raster rather than the other way around. I have included some code to help you identify the CRS for your r object. The CRS is a PROJ4 string, which R is moving away from in favour of EPSG codes. Fortunately, there are online resources such as https://epsg.io/ where you can search for equivalent codes (in this case, EPSG:3035). If you need the output in WGS84/EPSG:4326, comment below and I will update my answer.

    library(countrycode)
library(geodata)
library(terra)
library(sf)
library(ggplot2)
library(tidyterra)

# Shapefile converted to sf object (no need to assign CRS as it already has one)
ISO <- countrycode("Belgium", origin = 'country.name', destination = 'iso3c')
shp <- gadm(country=ISO, level=0, 
            path = tempdir(), 
            version="latest") %>%
  st_as_sf()

# Raster
r <- rast("ForestClim_01.tif")

# Return CRS of r
crs(r, proj = TRUE)
# [1] "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"

# Project shp to match CRS of ForestClim_01.tif using EPSG code found @ https://epsg.io/
shp3035 <- st_transform(shp, 3035)

# First try crop and mask combined
cr <- crop(r, shp3035, mask = TRUE, overwrite = TRUE)

# If crop and mask combined does not work, then run them separately
cr <- crop(r, shp3035, overwrite = TRUE)
cr <- mask(cr, shp3035, overwrite = TRUE)

ggplot() +
  geom_spatraster(data = cr) +
  geom_sf(data = shp3035, fill = NA) +
  scale_fill_continuous(na.value = "transparent")

    result