Search code examples
rterra

Change the crs of a raster to match the crs of a simple feature point object

Say I have a MODIS LAI raster and a simple feature (sf) object (points) with two different crs. I need to transform the crs of the raster to match the crs of the sf data (vice versa would also be ok, but I’d prefer to get rid of the MODIS sinusoidal system). Various attempts did not get me to the solution, i.e., having the raster and sf in the same CRS for further visualization, processing, etcetera. Any help on how can I do this, please? Below are a few options I tried:

library(terra)
library(sf)

lai <- terra::rast("lai.tif")
terra::crs(lai) #check crs of lai

plots <- readRDS("plots.rds")
sf::st_crs(plots) #check crs of sf

#Trial 1
plots_tr <- sf::as_Spatial(plots) #transform the sf and get the crs
lai_pr1 <- terra::project(lai, "+proj=utm +zone=31 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs") #project the raster using the crs of the transformed sf
terra::plot(lai_pr1) #plot raster
plot(plots_tr, add=T, pch=19, col="black", cex=0.1) #overlay points

#Trial 2
lai_pr2 <- lai
terra::crs(lai_pr2) <- "+proj=utm +zone=31 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" # replace lai crs
terra::plot(lai_pr2) #plot raster
plot(plots_tr, add=T, pch=19, col="black", cex=0.1) #overlay points

#Trial 3
lai_pr3 <- terra::project(lai, "epsg:25831") #project the raster using the EPSG code of the sf data (i.e., ETRS89 / UTM zone 31N)
terra::plot(lai_pr3) #plot raster
plot(plots_tr, add=T, pch=19, col="black", cex=0.1) #overlay points

Link to the data to run this example:

lai <- link removed

plots <- link removed

Maybe the issue is that the crs of the sf is WKT while the crs of the raster isn't? Also, using the sf without sf::as_Spatial() did not change the outcome. I also had a go at changing the crs of the sf with st_transform(). I read the many questions already posted on this but could not get the solution. Maybe spTransform()?

I hope this is not too silly. Many thanks.

Solution

  • Having fixed the coordinate reference system of your MODIS file, this is pretty straightforward making use of sf::st_transform() or terra::project(), depending on what crs you want to keep. For analysis, I'd continue with reprojecting vector data and keep the raster as-is, but feel free to get rid of the MODIS projection if you like (just don't forget about resampling). Also, absolutely no need to use {sp}.

    library(terra)
#> terra 1.7.71
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

lai <- terra::rast("lai.tif")
plots <- readRDS("plots.rds")

# fix MODIS crs
modis_crs <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs"
terra::crs(lai) <- modis_crs

# option 1: reproject vector data 
plots_reproj <- sf::st_transform(plots, modis_crs)

terra::plot(lai, main = "MODIS Sinusoidal")
plot(plots_reproj, add = TRUE, pch = 19, col="black", cex = 0.1)


    
# option 2: reproject raster data
lai_reproj <- terra::project(lai, "epsg:25831")

terra::plot(lai_reproj, main = "EPSG:25831")
plot(plots, add = TRUE, pch = 19, col="black", cex = 0.1)

    Created on 2024-04-27 with reprex v2.1.0