projectRaster projections seem to be wrong (maybe related to PROJ6 ?)

[GillesSanMartin]

It seems that projectRaster fails to reproject properly when provided a proj4 string or an sp CRS object.

Possibly related to the changes involved in PROJ6 and discarded datum in proj4 string ? #78

Example from this SO exchange : https://stackoverflow.com/a/63812649/2417437

raster_3.3-13    sp_1.4-2      sf_0.9-5
GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
library(sp)
library(raster)

# create a raster
r <- raster::raster(system.file("external/test.grd", package="raster"))

# create an sf spatial object
coo <- data.frame(x = c(246000, 247000, 246500), y = c(184000, 186000, 185000))
SF <- st_as_sf(coo, coords = c("x", "y"), crs = 31370)

# create an equivalent sp object
SP <- coo
coordinates(SP) <- ~x+y
proj4string(SP) <- CRS(SRS_string = "EPSG:31370")

Reprojecting the points towards the raster crs seems to work

SF_to_r <- st_transform(SF, crs = raster::crs(r))
raster::extract(r, SF_to_r) # result (correct) : 351.7868 236.4216 309.0073

Reprojecting the raster towards the crs of the points seem to create a new raster with the wrong crs
Extract provides a different answer that seems to be incorrect.

# CRS object
r_to_sp <- raster::projectRaster(r, crs = sp::CRS("+init=epsg:31370"))
raster::extract(r_to_sp, SP) # result (wrong) : 341.6399 222.1028 301.2286

# same result with a slightly different syntax for CRS
r_to_sp <- raster::projectRaster(r, crs = sp::CRS(SRS_string = "EPSG:31370"))
raster::extract(r_to_sp, SP) # result (wrong) : 341.6399 222.1028 301.2286

# proj4 string
r_to_sf <- raster::projectRaster(r, crs = sf::st_crs(31370)$proj4string)
raster::extract(r_to_sf, SF) # result (wrong) : 341.6399 222.1028 301.2286

r_to_sp <- raster::projectRaster(r, crs = sp::CRS("+init=epsg:31370")@projargs)
raster::extract(r_to_sp, SF) # result (wrong) : 348.5775 329.1199 277.2260

With stars we obtain the same answer as when we reproject the points toward the crs of the raster

library(stars)
STARS <- stars::read_stars(system.file("external/test.grd", package="raster"))

# reproject the points into the same crs as the stars raster
SF_to_STARS <- st_transform(SF, crs = st_crs(STARS))
aggregate(STARS, SF_to_STARS, function(x) x[1], as_points = FALSE)$test.grd # result = 351.7868 236.4216 309.0073

# reproject the stars raster into the same crs as the points
STARS_to_SF <- st_transform(STARS, crs = st_crs(SF))
aggregate(STARS_to_SF, SF, function(x) x[1], as_points = FALSE)$test.grd # result = 351.7868 236.4216 309.0073

[rhijmans]

This issue, and your post on SO are interesting and useful, but a bit difficult to respond to, as they cover a number of different things. I will try to respond to some of them. For now: transforming a raster will lead to different values. That is expected, and thus that is a bad approach when extracting values with "vector" data. This does not mean that the values you get are right or wrong. But the fact that they are different does not mean they are wrong. Again, that is expected. Stars is different. It appears that st_transform creates a curvilinear grid, not a raster in the strict (rectangular) sense, .

[GillesSanMartin]

OK thanks for the explanation. I have very little experience with stars (first time I use it in fact).
Because I get different results depending on the approach I usually used in the past with raster sp and sf only, I was trying to test different approaches to understand what was going wrong (hence the additional test with stars).
The initial problem was that I obtain different results when I reproject the points toward the crs of the raster or when I reproject the raster toward the crs of the points (if I use sf but also if i use sp as shown on StackOverflow). I would expect to find exactly the same results whatever the approach used. For several reasons I have the impression that the values I obtain when I reproject the points toward the crs of the raster are the "correct" ones but I'm not totally certain of that.

[rhijmans]

Projecting a raster to a new (rectangular) raster requires estimation of the new cell values based on the overlap of the original cell values. The difference is Illustrated here:

library(raster)
r <- raster::raster(system.file("external/test.grd", package="raster"))
a <- aggregate(r, 10, na.rm=T)
b <- as(a, "SpatialPolygons")
aa <- projectRaster(a, crs = "+init=epsg:31370")
bb <- spTransform(b, "+init=epsg:31370")
plot(aa)
lines(bb)