raster - 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 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:

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.

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 (or SpatialLines object with conversion) to a pixel image (or RasterLayer 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 the dimyx parameter, which sets the dimensions (number of rows and columns).

require(sp)
require(raster)
require(maptools)
require(spatstat)

# 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))


# Convert SpatialLines to psp object using maptools library
pspSl <- as.psp(sl)
# Pixellate with resolution of 0.5, i.e. 2x2 pixels
px <- pixellate(pspSl, eps=0.5)
# This can be converted to raster as desired
rLength <- raster(px)
# Values:
round(as.matrix(rLength),3)


     [,1] [,2]
[1,]  0.5  0.0
[2,]  1.0  0.5

The results are exactly as desired.

Adrian also answered my questions about spatstat::density.psp(). He explains that this function:

computes the convolution of the Gaussian kernel with the lines. Intuitively, this means that density.psp 'smears' the lines into two-dimensional space. So density(L) is like a blurred version of pixellate(L). In fact density(L) is very similar to blur(pixellate(L)) where blur is another spatstat function that blurs an image. [The parameter] sigma is the bandwidth of the Gaussian kernel. The value of density.psp(L) at a given pixel u, is something like the total amount of line length in a circle of radius sigma around the pixel u, except that it's really a weighted average of such contributions from different circle radii. Units are length^(-1), i.e. line length per unit area.

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 in pixellate(). I guess I'll have to leave that for the experts.