I have a large (~70MB) shapefile of roads and want to convert this to a raster with road density in each cell. Ideally I'd like to do this in R along with GDAL command line tools if necessary.
My initial approach was to directly calculate the lengths of line segments in each cell as per this thread. This produces the desired results, but is quite slow even for shapefiles much smaller than mine. Here's a very simplified example for which the correct cell values are obvious:
require(sp)
require(raster)
require(rgeos)
require(RColorBrewer)
# Create some sample lines
l1 <- Lines(Line(cbind(c(0,1),c(.25,0.25))), ID="a")
l2 <- Lines(Line(cbind(c(0.25,0.25),c(0,1))), ID="b")
sl <- SpatialLines(list(l1,l2))
# Function to calculate lengths of lines in given raster cell
lengthInCell <- function(i, r, l) {
r[i] <- 1
rpoly <- rasterToPolygons(r, na.rm=T)
lc <- crop(l, rpoly)
if (!is.null(lc)) {
return(gLength(lc))
} else {
return(0)
}
}
# Make template
rLength <- raster(extent(sl), res=0.5)
# Calculate lengths
lengths <- sapply(1:ncell(rLength), lengthInCell, rLength, sl)
rLength[] <- lengths
# Plot results
spplot(rLength, scales = list(draw=TRUE), xlab="x", ylab="y",
col.regions=colorRampPalette(brewer.pal(9, "YlOrRd")),
sp.layout=list("sp.lines", sl),
par.settings=list(fontsize=list(text=15)))
round(as.matrix(rLength),3)
#### Results
[,1] [,2]
[1,] 0.5 0.0
[2,] 1.0 0.5
Looks good, but not scaleable! In a couple other questions the spatstat::density.psp()
function has been recommended for this task. This function uses a kernel density approach. I am able to implement it and it seem faster than the above approach, but I'm unclear how to choose the parameters or interpret the results. Here's the above example using density.psp()
:
require(spatstat)
require(maptools)
# Convert SpatialLines to psp object using maptools library
pspSl <- as.psp(sl)
# Kernel density, sigma chosen more or less arbitrarily
d <- density(pspSl, sigma=0.01, eps=0.5)
# Convert to raster
rKernDensity <- raster(d)
# Values:
round(as.matrix(rKernDensity),3)
#### Results
[,1] [,2]
[1,] 0.100 0.0
[2,] 0.201 0.1
I thought it might be the case that the kernel approach calculates density as opposed to length per cell, so I converted:
# Convert from density to length per cell for comparison
rKernLength <- rKernDensity * res(rKernDensity)[1] * res(rKernDensity)[2]
round(as.matrix(rKernLength),3)
#### Results
[,1] [,2]
[1,] 0.025 0.000
[2,] 0.050 0.025
But, in neither case, does the kernel approach come close to aligning with the more direct approach above.
So, my questions are:
- How can I interpret the output of the
density.psp
function? What are the units? - How can I choose the
sigma
parameter indensity.psp
so the results align with the more direct, intuitive approach above? - Bonus: what is the kernel line density actually doing? I have some sense for how these approaches work for points, but don't see how that extends to lines.