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 <- https://github.com/frandadamo/frandadamo/blob/main/lai.tif

plots <- https://github.com/frandadamo/frandadamo/blob/main/plots.rds

Maybe the issue is that the crs of the sf is WKT while the crs of 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.