Rasterize lines in R: How to assign max value to a cell?

Question

Problem

I want to rasterize spatial lines and assign the maximum value of all lines that touch/intersect a raster cell to that cell.

I work in terra with lines as a SpatVect object, I would hence prefer a terra::rasterize solution. I would also be happy with a solution using stars::st_rasterize (see also this question).

The documentation of terra::rasterize seems to suggest it does not support lines properly so far - using the fun argument seems to be limited to points vector data only. I tried with terra::rasterize nevertheless, see the example below.

I'm aware of raster::rasterize, but it seems a bit outdated since it's still sp object based. Or is it the only way to do it atm?

Example

Here you can see that neither the max nor the mean function seems to work properly when rasterizing via terra::rasterize(... fun = "max"/"mean"):

library("terra")

### Example data ###
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
lns <- as.lines(v)
r <- rast(v, res=.2)

### Rasterize via terra::rasterize ###
x_max <- rasterize(lns, r, fun="max", field = "POP")
x_mean <- rasterize(lns, r, fun="mean", field = "POP")

### Plot results ###
fontsize <- 0.7
par(mfrow=c(1,3))
plot(lns, y = "POP", main = "Lines with original values")
text(lns, "POP", cex = fontsize)
plot(x_max, main = "Rasterized via fun 'max'")
text(x_max, cex = fontsize)
plot(lns, add = T)
plot(x_mean, main = "Rasterized via fun 'mean'")
text(x_mean, cex = fontsize)
plot(lns, add = T)

Answer 1

I have found a somewhat hacky solution. As proposed in the comments, I turned the lines into points by sampling along them.

I used sf::st_line_sample() for this instead of rgeos::gInterpolate() - it was cleaner and easier. Then terra::rasterize() can handle the points correctly and applies the max fun as expected.

Confirmed by the plot below.

Example solution

library("terra")
library("dplyr")
library("sf")
library("units")

### Example data ###
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
lns <- as.lines(v)
r <- rast(v, res=.2)

### Turn SpatVector lines into sf object
lns_sf <- lns %>% 
  st_as_sf() %>% 
  st_transform(2169) # reprojection needed for st_line_sample

### Sample points along all lines every kilometer
pts_geometries <- lns_sf %>% 
  st_line_sample(density = units::set_units(1, 1/km))

### Add attributes and make MULTIPOINTs simple POINTS
pts_sf <- st_sf(lns_sf,
                geometry = pts_geometries) %>% 
  st_cast("POINT") %>% 
  st_transform(crs(lns))

### Go back to terra: Turn sf into SpatVector object
pts <- pts_sf %>% 
  vect()

### Now rasterization works and "max" function is applied corretly
x_max <- rasterize(pts, r, fun="max", field = "POP")

fontsize <- 0.7
par(mfrow=c(1,3))
plot(lns, y = "POP", main = "Lines with original values")
text(lns, "POP", cex = fontsize)
plot(x_max, main = "Rasterized with fun 'max' using points generated from lines")
text(x_max, cex = fontsize)
plot(lns, add = T)
plot(pts, main = "Points used for rasterization")