Convert line shapefile to raster, value=total length of lines within cell

The below is modified from Jeffrey Evans' solution. This solution 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

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:

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

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)
    return(roads_utm_crp_length)
  } else {
    return(0)
  }
})

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

library(RColorBrewer)
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))