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.

  • Maybe vignette('over', package = 'sp') might help.
    – user32309
    Commented Oct 28, 2014 at 10:38

6 Answers 6


You can now do this with terra::rasterizeGeom

Example data:

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


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


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
  • Seems the most efficient and elegant method of those listed. Also, hadn't seen raster::intersect() before, I like that it combines the attributes of the intersected features, unlike rgeos::gIntersection().
    – Matt SM
    Commented Mar 19, 2015 at 16:53
  • +1 always nice seeing more efficient solutions! Commented Mar 2, 2016 at 18:19
  • @RobertH, I tried using this solution for another problem, where I want to do the same as asked in this thread, but with a very large shapefile of roads (for an entire continent). It seems to have worked, but when I try to do the figure as done by @ fdetsch, I do not get a contiguous grid but just a few colored squares in the picture (see in tinypic.com/r/20hu87k/9).
    – Doon_Bogan
    Commented May 23, 2016 at 11:42
  • And by most efficient... with my sample dataset: run time 0.6sec for this solution vs. 8.25sec for the most upvotes solution. Commented Dec 21, 2018 at 19:10
  • 1
    For my use case it takes to much memory.
    – Funkwecker
    Commented Jun 26, 2020 at 6:23

Following a recent question, you may want to make use of the functionalities offered by the rgeos package to solve your problem. For reasons of reproducibility, I downloaded a shapefile of Tanzanian roads from DIVA-GIS and put it in my current working directory. For the upcoming tasks, you will need three packages:

  • rgdal for general spatial data handling
  • raster for rasterization of the shapefile data
  • rgeos to check intersection of roads with raster template and calculate road lengths

Consequently, your first lines of could should look like this:


After that, you need to import the shapefile data. Note that DIVA-GIS shapefiles are distributed in EPSG:4326, so I will project the shapefile to EPSG:21037 (UTM 37S) to deal with meters rather than degrees.

roads <- readOGR(dsn = ".", layer = "TZA_roads")
roads_utm <- spTransform(roads, CRS("+init=epsg:21037"))

For subsequent rasterization, you will need a raster template that covers the spatial extent of your shapefile. The raster template consists of 10 rows and 10 columns by default, thus avoiding too extensive computation times.

roads_utm_rst <- raster(extent(roads_utm), crs = projection(roads_utm))

Now that the template is set up, loop through all cells of the raster (that currently consists of NA values only). By assigning a value of '1' to the current cell and subsequently executing rasterToPolygons, the resulting shapefile 'tmp_shp' automatically holds the extent of the currently processed pixel. gIntersects detects whether this extent overlaps with roads. If not, the function will return a value of '0'. Otherwise, the road shapefile is cropped by the current cell and the total length of 'SpatialLines' within that cell is being calculated using gLength.

lengths <- sapply(1:ncell(roads_utm_rst), function(i) {
  tmp_rst <- roads_utm_rst
  tmp_rst[i] <- 1
  tmp_shp <- rasterToPolygons(tmp_rst)

  if (gIntersects(roads_utm, tmp_shp)) {
    roads_utm_crp <- crop(roads_utm, tmp_shp)
    roads_utm_crp_length <- gLength(roads_utm_crp)
  } else {

Finally, you can insert the calculated lengths (which are converted to kilometers) into the raster template and visually verify your results.

roads_utm_rst[] <- lengths / 1000

spplot(roads_utm_rst, scales = list(draw = TRUE), xlab = "x", ylab = "y", 
       col.regions = colorRampPalette(brewer.pal(9, "YlOrRd")), 
       sp.layout = list("sp.lines", roads_utm), 
       par.settings = list(fontsize = list(text = 15)), at = seq(0, 1800, 200))


  • This is awesome! I changed sapply() to pbsapply() and used the cluster argument cl = detectCores()-1. Now I'm able to run this example in parallel! Commented Nov 17, 2018 at 0:18

Here's yet another approach. It deviates from those already given by using the spatstat package. As far as I can tell, this package has its own version of spatial objects (e.g. im vs. raster objects), but the maptools package allows conversion back and forth between spatstat objects and standard spatial objects.

This approach is taken from this R-sig-Geo post.


# Load data and transform to UTM
roads <- shapefile('data/TZA_roads.shp')
roadsUTM <- spTransform(roads, CRS("+init=epsg:21037"))

# Need to convert to a line segment pattern object with maptools
roadsPSP <- as.psp(as(roadsUTM, 'SpatialLines'))

# Calculate lengths per cell
roadLengthIM <- pixellate.psp(roadsUTM, dimyx=10)

# Convert pixel image to raster in km
roadLength <- raster(dtanz / 1000, crs=projection(roadsUTM))

# Plot
spplot(rtanz, scales = list(draw=TRUE), xlab="x", ylab="y", 
       col.regions=colorRampPalette(brewer.pal(9, "YlOrRd")), 
       sp.layout=list("sp.lines", roadsUTM), 
       par.settings=list(fontsize=list(text=15)), at=seq(0, 1800, 200))


The slowest bit is converting the roads from SpatialLines to a Line Segment Pattern (i.e. spatstat::psp). Once that's done the actual length calculation parts is quite quick, even for much higher resolutions. For example, on my old 2009 MacBook:

system.time(pixellate(tanzpsp, dimyx=10))
#    user  system elapsed 
#   0.007   0.001   0.007

system.time(pixellate(tanzpsp, dimyx=1000))
#    user  system elapsed 
#   0.146   0.032   0.178 
  • Hmmmm...I sure wish those axes weren't in scientific notation. Anyone have any idea how to fix that?
    – Matt SM
    Commented Mar 19, 2015 at 17:31
  • You can modify the global R settings and turn off notation using: options(scipen=999)) however, I do not know if lattice will honor the global environment settings. Commented Mar 2, 2016 at 18:25

You do not need a for loop. Just intersect everything at once and then add line lengths to the new line segments using the "SpatialLinesLengths" function in sp. Then, using the raster package rasterize function with the fun=sum argument you can create a raster with the sum of the line length(s) intersecting each cell. Using the above answer and associated data here is code that will generate the same results.


roads <- readOGR(getwd(), "TZA_roads")
  roads <- spTransform(roads, CRS("+init=epsg:21037"))
    rrst <- raster(extent(roads), crs=projection(roads))

# Intersect lines with raster "polygons" and add length to new lines segments
rrst.poly <- rasterToPolygons(rrst)
  rp <- gIntersection(roads, rrst.poly, byid=TRUE)
    rp <- SpatialLinesDataFrame(rp, data.frame(row.names=sapply(slot(rp, "lines"), 
                               function(x) slot(x, "ID")), ID=1:length(rp), 
                               length=SpatialLinesLengths(rp)/1000) ) 

# Rasterize using sum of intersected lines                            
rd.rst <- rasterize(rp, rrst, field="length", fun="sum")

# Plot results
spplot(rd.rst, scales = list(draw=TRUE), xlab="x", ylab="y", 
       col.regions=colorRampPalette(brewer.pal(9, "YlOrRd")), 
       sp.layout=list("sp.lines", rp), 
       par.settings=list(fontsize=list(text=15)), at=seq(0, 1800, 200))
  • First time I see SpatialLinesLengths. Guess it is never too late to learn, thank you (: rasterize takes quite long, though (7 times longer than the upper approach on my machine).
    – fdetsch
    Commented Nov 5, 2014 at 8:27
  • I did notice that rasterize was being slow. However, for large problems a for loop is going to really slow things down and I think that your will see a much more optimized solution with rasterize. Besides the developer of the raster package is striving to make every release more optimized and faster. Commented Nov 5, 2014 at 15:27
  • One potential issue I found with this technique is that the rasterize() function includes all lines that touch a given cell. This results in line segment lengths being counted twice in some cases: once in the cell they're supposed to and once in the adjacent cell that the line endpoint just touches.
    – Matt SM
    Commented Mar 14, 2015 at 19:09
  • Yes, but "rp" is the object that is being rasterized which it the intersection of polygons and points. Commented Mar 14, 2015 at 21:22

Let me present you the package vein with several functions to work spatial lines and imports sf and data.table

netsf <- st_as_sf(net) #Convert Spatial to sf
netsf <- st_transform(netsf, 31983) # Project data
netsf$length_m  <- st_length(netsf)
netsf <- netsf[, "length_m"]
g <- make_grid(netsf, width = 1000) #Creat grid of 1000m spacing with columns id for each feature
# Number of lon points: 12
# Number of lat points: 11

gnet <- emis_grid(netsf, g)

sf_to_raster <- function(x, column, ncol, nrow){
  x <- sf::as_Spatial(x)
  r <- raster::raster(ncol = ncol, nrow = nrow)
  raster::extent(r) <- raster::extent(x)
  r <- raster::rasterize(x, r, column)

rr <- sf_to_raster(gnet, "length_m", 12, 11)
spplot(rr, sp.layout = list("sp.lines", as_Spatial(netsf)),
       col.regions = cpt(), scales = list(draw = T))
spplot(rr, sp.layout = list("sp.lines", as_Spatial(netsf)),
       col.regions = cpt(pal = 5176), scales = list(draw = T))
spplot(rr, sp.layout = list("sp.lines", as_Spatial(netsf)),
       col.regions = lucky(), scales = list(draw = T))
# Colour gradient: neota_flor_apple_green, number: 6165

enter image description here

enter image description here

enter image description here


This may sound a bit naive but if it is a road system then select the roads and save to a clip board then look find a tool that allows you to add a buffer to the clipboard set it to the legal width of the road ie 3 meters +/- remember that the buffer is from center line to edge * 2 i for each side so a 3 meter buffer is actually a 6 meter road from side to side.

  • What does road width have to do with road length. This answer does not address the question. Commented Mar 2, 2016 at 18:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.