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

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.

Answer

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