5

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)
josliber
  • 43,891
  • 12
  • 98
  • 133
Marine
  • 521
  • 1
  • 4
  • 23

2 Answers2

5

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
rsoren
  • 4,036
  • 3
  • 26
  • 37
  • Thank you very much rsoren. In particular, I have a problem to calculate road length (kilometers) because the function `rasterToContour()` gives road edges with edges on left side and right side of the road. I would prefer to have a line passing through a road (i.e., passing through the center of road pixels). – Marine Jan 29 '16 at 17:35
  • Could you not just divide by 2, to calculate the length of one edge? – rsoren Jan 29 '16 at 23:11
  • how will you distinguish a road inside the buffer from a road in the edge of the buffer? in one case you have to divide it by 2, in the other, you can't. – Carlos Alberto Jan 29 '16 at 23:23
4

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.

Community
  • 1
  • 1
Carlos Alberto
  • 598
  • 3
  • 9