I have a large (~70MB) shapefile of roads and want to convert this to a raster with road density in each cell.
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.
Best Answer
I posted this question on the R-sig-Geo listserv and received a helpful answer from Adrian Baddeley, one of the spatstats authors. I will post my interpretation of his response here for posterity.
Adrian notes that the function
spatstat::pixellate.psp()
is a better match to my task. This function converts a line segment pattern (orSpatialLines
object with conversion) to a pixel image (orRasterLayer
with conversion), where the value in each cell is the length of the line segments passing through that cell. Exactly what I'm looking for!The resolution of the resulting image can be defined with the
eps
parameter or thedimyx
parameter, which sets the dimensions (number of rows and columns).The results are exactly as desired.
Adrian also answered my questions about
spatstat::density.psp()
. He explains that this function:It remains somewhat unclear to me when the Gaussian kernel approach of
density.psp()
would be preferred over the more intuitive approach of directly calculating line lengths inpixellate()
. I guess I'll have to leave that for the experts.