Search code examples
rterra

Problem at cropping raster data using spatVector

I’d have a very basic question to ask, please. I have to crop a number of MODIS LAI raster using a SpatVector. Using the code below, the cropped raster and the SpatVector do not align so I guess that the cropping is not executed properly. Any ideas? Thanks a lot.

library(terra)
library(geodata)

#get spatVector
es <- geodata::gadm(country="ESP", level=1, path=".")
catal <- es[es$NAME_1 == "Cataluña", ]

#get MODIS data
modis <- terra::rast("see link below")
modis <- modis*0.1
# terra::plot(modis)

modis_crs <- terra::crs(modis) #get crs of MODIS data
catal_2 <- terra::project(catal, modis_crs) #give the spatVector the same crs as MODIS
# terra::crs(catal_2)==terra::crs(modis) #check if the two crs are the same
modis_cr <- terra::crop(modis, catal_2) #crop (here without mask)
terra::plot(modis_cr) #plot cropped MODIS data
terra::plot(catal_2, add=TRUE, lwd=2) #plot spatVector with the MODIS crs

A MODIS LAI raster to run this example is available here <- link removed

This is what I get, the two files do not align (ignore the colour of the raster):

enter image description here

I tried to change the crs of MODIS data rather than the crs of the spatVector, I tried with different options of extent, align, or snap. I tried to play a bit with the sf library and read quite a few similar questions.

Solution

  • It seems like the crs definition of your raster is not correct. Overriding this by the proj4string below results in expected positioning:

    library(terra)
#> terra 1.7.71
library(geodata)

# get data
es <- geodata::gadm(country="ESP", level=1, path=".")
catal <- es[es$NAME_1 == "Cataluña", ] 

# read modis
modis <- terra::rast("cropmodis.tif")

# replace modis crs
terra::crs(modis) <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs"

# reproject vector
modis_crs <- terra::crs(modis)
catal_2 <- terra::project(catal, modis_crs)

# crop
modis_cr <- terra::crop(modis, catal_2)

# inspect
terra::plot(modis_cr)
terra::plot(catal_2, add = TRUE)

    Note:

    The crs definition comes from a random MODIS dataset downloaded from EarthExplorer:

    terra::rast("MOD11A2.A2024089.h18v04.061.2024099144843.hdf") |> terra::crs(proj = TRUE)
#> [1] "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs"