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)

You can also omit `fun = ` and get exactly the same result. Seems like the documentation is correct - only used if `x` consists of points. What about a workaround like `p <- as.points(lns)` and use these as input to `rasterize()`? Might work in this case, considering the resolution of your raster and the complexity of the geometry used. — dimfalk, Jul 28 '22 at 18:58
Eventually in combination with an increase in node density along the lines using `rgeos::gInterpolate()` on a `sp` object? — dimfalk, Jul 28 '22 at 19:47
My real world data is far more complex (44 000 lines, raster with 700 000 cells). I just used this simple example for repex purposes. — chamaoskurumi, Jul 29 '22 at 12:18
Complex (in terms of nodes per grid cell) would've been good for `as.points()` but I assume you meant quantities only. However, `rgeos::gInterpolate()` might still work here when neglecting performance issues. — dimfalk, Jul 29 '22 at 12:41
It worked, performance was not an issue. I used `sf::st_line_sample` instead of `rgeos::gInterpolate`. Thanks @falk-env! Oh and you were right, by complex by meant quantities only. — chamaoskurumi, Aug 02 '22 at 14:51

score 1 · Accepted Answer · answered Aug 02 '22 at 14:50

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")

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

Problem

Example

1 Answers1

Example solution