Calculate line density within a buffer from a raster

Question

Is there a way to calculate road density (km/km²) within a buffer around spatial lines? The roads are represented by pixels (1 pixel = 625 m²) in a raster. So I began to convert the road pixels into polylines by using the function rasterToContour (package raster). Then, I'm thinking of calculating total length of lines within the buffer (in km) and the buffer area (in km²).

r <- raster(rast_path)
x <- rasterToContour(r)

Here is a reproducible example:

## To create raster:
library(raster)
library(rgeos)
r <- raster(ncols=90, nrows=50)
values(r) <- sample(1:10, ncell(r), replace=TRUE)

## Road raster
r[r[] < 10] <- 0
r[r[] >= 10] <- 1
plot(r)

## To create spatial lines 
line1 <- rbind(c(-125,0), c(0,60))
line2 <- rbind(c(0,60), c(40,5))
line3 <- rbind(c(40,5), c(15,-45))
line1_sp <- spLines(line1)
line2_sp <- spLines(line2)
line3_sp <- spLines(line3)

## To create buffer around lines
line2_buff <- gBuffer(line2_sp, width=20)
plot(line2_sp,add=T)
plot(line2_buff,add=T)

Answer 1

You're looking for road length (kilometers) divided by buffer area (square kilometers), right? To calculate road length, you can use your method with rasterToContour(), although that isn't reproducible with the example you provided.

But to calculate area of the buffer in square kilometers, you can do: n <- length(extract(r, line2_buff)) to get the n number of pixels in the buffer, which are 625 m^2 each. The conversion factor you need is 1 km^2 = 1,000,000 m^2. All together, area of the buffer in km^2 is given by:

length(extract(r, line2_buff)) * 625 / 1000000

Answer 2

The function that you are looking to get the intersection between a line and a polygon is gIntersection (see this link http://robinlovelace.net/r/2014/07/29/clipping-with-r.html). However, if you sum the border of a road crossing a buffer you will be counting twice the length of the road (left side + right side of the road).

The problem is that converting a raster to a road map (lines) is not as straightforward as to converting it to a polygon (what you will obtain with rasterToContour). And converting to a polygon will not give you the result you're looking for (as previously explained). So, you will have to to do it manually or invest some extra time in coding (search for "skeletonizing raster", for instance Identify a linear feature on a raster map and return a linear shape object using R).

I think the usual approach is to work on an area basis (km2/km2) and you can do it pretty easily on raster format. If your resolution is good enough, you can do later (Area of the roads)/(average width of a road)/(buffer area) to try to approximate the value to km/km2.