R – Convert Line Shapefile to Raster with Total Length of Lines Within Cell

linerrasterrasterization

I have a line shapefile representing a road network. I wish to rasterize this data, with the resulting values in the raster showing the total length of the lines that fall within the raster cell.

The data is in the British National Grid projection so the units will be meters.

Ideally I would like to perform this operation using R, and i'm guessing that the rasterize function from the raster package would play a part in achieving this, I just can't work out what the function applied should be.

Best Answer

You can now do this with terra::rasterizeGeom

Example data:

library(terra)
v <- vect(system.file("ex/lux.shp", package="terra"))
r <- rast(v)

Solution:

x <- rasterizeGeom(v, r, "length")    
plot(x); lines(v)

The below is a solution with the older spatial packages. It is based on Jeffrey Evans' answer, but this one is much faster as it does not use rasterize

library(raster)
library(rgdal)
library(rgeos)

roads <- shapefile("TZA_roads.shp")
roads <- spTransform(roads, CRS("+proj=utm +zone=37 +south +datum=WGS84"))
rs <- raster(extent(roads), crs=projection(roads))
rs[] <- 1:ncell(rs)

# Intersect lines with raster "polygons" and add length to new lines segments
rsp <- rasterToPolygons(rs)
rp <- intersect(roads, rsp)
rp$length <- gLength(rp, byid=TRUE) / 1000
x <- tapply(rp$length, rp$layer, sum)
r <- raster(rs)
r[as.integer(names(x))] <- x
Related Question