Convert a vector grid created with sf::st_make_grid into a raster stack

Question

I have a vector grid created with sf::st_make_grid() with several columns of values. I would like to convert it to a raster stack, with each column of values becoming a layer of the raster stack.

I provide an example below which is close to what I want to achieve. I have 2 questions problems :

1. is it a good way to approach that problem? Isn't it a more straightforward way to convert what is basically a vectorial grid into pixels for a raster?
2. when I plot the vector grid and the raster on top of each other, the grid cells are not aligned with the pixels of the raster (see last line using mapview.) I would expect that at least their centroid would be aligned. So, I guess I am doing something wrong.

library(sf)
library(terra)

my_cellsize <- 50000

# Load the 'nc' shapefile from the sf package
nc <- st_read(system.file("shape/nc.shp", package = "sf")) |> 
    st_transform(crs = 3359) # 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))
my_grid$value2 <- rnorm(nrow(my_grid))
my_grid$value3 <- 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")
value2_raster <- rasterize(vect(my_grid), template, field = "value2")
value3_raster <- rasterize(vect(my_grid), template, field = "value3")
stacked_raster <- c(value1_raster, value2_raster, value3_raster)


# interactivbe map showing that the 2 set of grid cells do not overlap correctly
library(mapview)
mapview(stars::st_as_stars(value1_raster)) + mapview(my_grid)

Here is a screenshot of the leaflet map showing the misalignment between the vector grid and raster pixels

I also tried to start from the xy coordinates of the centroids of the cells. But if I understand well, this would work only if I had values on a square grid. Here the cells that are outside the area of interest are considered as missing values and trigger an error message (missing values not allowed)

data.frame(
    my_grid |> st_centroid() |> st_coordinates(),
    my_grid[,c("value1", "value2", "value3"), drop = TRUE]
    ) |> 
    rast(type = "xyz", extent= ext(my_grid), crs = "EPSG:3359")

I also tried with {stars} with no more success :

template = st_as_stars(st_bbox(my_grid), dx = my_cellsize, dy= my_cellsize, values = NA_real_)
value1_raster = st_rasterize(my_grid, template)
mapview(stars::st_as_stars(value1_raster)) + mapview(my_grid)

Answer 1

There is no misalignment in the data. You can see that with

plot(value1_raster)
lines(my_grid)

You see misalignment because you use "mapview"

mapview((value1_raster)) + mapview(my_grid)

This is because mapview first transforms the data to the coordinate reference system of its base maps. After transformation, the polygon data is no longer a regular grid. But the raster data is transformed to a new regular grid, so the two sources are now misaligned. To avoid that you could transform nc to the crs that mapview uses. For example:

library(sf)
library(terra)
library(mapview)

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(value1_raster) + mapview(as.lines(vect(my_grid)))

Also, if your goal is to make a raster, then it seems overly complicated and inefficient to first make polygons and rasterize these. You could do something along these lines instead:

library(terra)
nc <- vect(system.file("shape/nc.shp", package = "sf")) |> 
    project(y = "EPSG:3359") 

r <- rast(nc, res=50000)
r$value1 <- rnorm(ncell(r))
r$value2 <- rnorm(ncell(r))
r$value3 <- rnorm(ncell(r))