Misalignment problem when curvilinear raster grid are ploted as regular ones

[GillesSanMartin]

I came across an unexpected behavior of mapview (?) : I create a vector grid with sf:st_make_grid then I transform it into a terra raster object. When I plot both with mapview they are misaligned.

After my question on stackoverflow, Robert Hijman suggested that it was a problem with the Mercator projection used in mapview that should produce a curvilinear grid but is represented with a regular grid (NB: I paraphrase with my own words and understanding)

Here is a reproducible example :

library(sf)
library(terra)
library(stars)
library(ggplot2)
library(mapview)

my_cellsize <- 50000

# Load the 'nc' shapefile from the sf package
nc <- st_read(system.file("shape/nc.shp", package = "sf")) |> 
    st_transform(crs = 32617) # projection in meters

# Create a grid within the extent of the 'nc' dataset
my_grid <-  st_make_grid(nc, cellsize = my_cellsize)
my_grid <- st_sf(GridID = 1:length(my_grid), my_grid)
my_grid <- my_grid[nc,]

# add some random values
my_grid$value1 <- rnorm(nrow(my_grid))

# convert the vector grid into a raster stack
template = rast(vect(my_grid), res = my_cellsize)
value1_raster <- rasterize(vect(my_grid), template, field = "value1")

Various calls to mapview showing the misalignment

# With a simple call to mapview the vector grid and the raster grid are not aligned
mapview(st_as_stars(value1_raster) ) + mapview(my_grid)

# Transforming to the Mercator projection used by mapview does not help
# But the misalignment is different
mapview(st_as_stars(value1_raster) |> st_transform(crs = "+proj=webmerc +datum=WGS84")) + 
    mapview(my_grid |> st_transform(crs = "+proj=webmerc +datum=WGS84"))

# using st_warp instead of st_transform seems to simply reproduce what is done by mapview in the backgroud
mapview(st_as_stars(value1_raster) |> st_warp(crs = "+proj=webmerc +datum=WGS84")) + 
    mapview(my_grid |> st_transform(crs = "+proj=webmerc +datum=WGS84"))

Various calls to other plotting approaches that seem to work fine :

# The raster and the vector grid are correctly 
# plotted with a simple call to plot on the terra objects : 

plot(value1_raster, reset = TRUE)
lines(my_grid, col = "red")

# Using stars, we can also have a correct alignment of a curvilinear raster grid with the vector
# This is what I would expect to see in mapview

stars::st_as_stars(value1_raster) |> 
    st_transform(crs = "+proj=webmerc +datum=WGS84") |> 
    plot(reset = FALSE)
my_grid |> st_transform(crs = "+proj=webmerc +datum=WGS84") |> st_geometry() |> 
    plot(add = TRUE, reset = FALSE, border = "red", lwd = 2, color = NA)

# Same with ggplot : 
ggplot() +
    geom_stars(data = st_transform(st_as_stars(value1_raster), 
                                   crs = "+proj=webmerc +datum=WGS84"), 
               alpha = 0.8, downsample = c(0, 0, 0)) +
    geom_sf(data = st_transform(my_grid, crs = "+proj=webmerc +datum=WGS84"), 
            color = "red", fill = NA) +
    theme_void()

Again, as pointed out by Robert Hijmans, one way to get a correct plot in mapview would be to use the Mercator projection from the beginning. But that does not seem to be a usable solution in real life where you might need other projections adapted to your local conditions to compute areas, distances, etc.

my_cellsize <- 50000
nc <- st_read(system.file("shape/nc.shp", package = "sf")) |> 
    st_transform(crs = "+proj=webmerc +datum=WGS84") 
my_grid <-  st_make_grid(nc, cellsize = my_cellsize)
my_grid <- st_sf(GridID = 1:length(my_grid), my_grid)
my_grid <- my_grid[nc,]
my_grid$value1 <- rnorm(nrow(my_grid))
template = rast(vect(my_grid), res = my_cellsize)
value1_raster <- rasterize(vect(my_grid), template, field = "value1")

mapview(st_as_stars(value1_raster)) + mapview(st_as_sf(as.lines(vect(my_grid))))