Calculate line density in a buffer from a raster

Question

Is there a way to calculate the density of the road (km / km²) in the buffer around the spatial lines? Roads are represented by pixels (1 pixel = 625 m²) in the raster. So I started converting road pixels into polylines using the rasterToContour (package raster) function. Then I think about calculating the total length of the lines in the buffer (in km) and the buffer area (in km²).

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

Here is an example of reproducibility:

 ## 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 are looking for the length of the road (kilometers) divided by the buffer area (square kilometers), aren't you? To calculate the length of the road, you can use your method with rasterToContour() , although this does not reproduce with the example you provided.

But to calculate the buffer area in square kilometers, you can do: n <- length(extract(r, line2_buff)) to get n number of pixels in the buffer equal to 625 m ^ 2. The required conversion factor is 1 km ^ 2 = 1 000 000 m 2. All together, the buffer area in km ^ 2 is determined as follows:

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

Answer 2

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

The problem is that converting a raster to a roadmap (strings) is not as simple as converting it to a polygon (which you get with a ToContour raster). And converting to a polygon will not give you the result you are looking for (as explained earlier). Thus, you will have to do this manually or spend additional time coding (searching for a “skeletonizing raster”, for example, Define a linear function on a raster map and return the object to a linear shape using R ).

I think the usual approach is to work based on the area (km2 / km2), and you can do this quite easily in raster format. If your resolution is good enough, you can do later (Road Area) / (average road width) / (buffer zone) to try to get closer to the km / km2 value.