[R-sig-Geo] Calculating road (i.e. linear feature) density using spatstat::density.psp()

[Matt Strimas-Mackey]

Note: I previously posted this question to Stack Exchange, but haven't
receive a response. I will update the SE question with any answers I get on
this listserv:
http://gis.stackexchange.com/questions/138861/calculating-road-density-in-r-using-kernel-density

I have a large (~70MB) shapefile of roads and want to convert this to a
raster with road density/length in each cell. Ideally I'd like to do this
in R.

My initial approach was to directly calculate the lengths of line segments
in each cell. 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, road lengths in each cell
     [,1] [,2]
[1,]  0.5  0.0
[2,]  1.0  0.5

--------------

Looks good, but not scaleable! In a previous post on this listserv, 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:
1. How can I interpret the output of the density.psp function? What are the
units?
2. How can I choose the sigma parameter in density.psp so the results align
with the more direct, intuitive approach above?
3. 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.

Thanks!

	[[alternative HTML version deleted]]