# Calculating road density in R using kernel density? [closed]

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) * res(rKernDensity)
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.

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.