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

[Adrian Baddeley]

Matt Strimas-Mackey <strimas at zoology.ubc.ca> writes:

> Note: I previously posted this question to Stack Exchange, 
> but haven't receive a response.

Questions about the spatstat package can best be addressed directly to the authors.

> 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.

spatstat::pixellate.psp will do this. 
Use the "..." arguments to specify the raster dimensions.
                L <- as.psp(YourData)
                Z <- pixellate(L, dimyx=512)
The resulting image values give the total length of fragments of lines in each pixel. 
To get values which are lengths per unit area, divide by the pixel area:
                Z <- Z/with(Z, xstep * ystep)

> In a previous post on this listserv, the spatstat::density.psp() function 
> has been recommended for this task. 

Uh, no, that function does something else. 

> 1. How can I interpret the output of the density.psp function? 


As the help file states, density.psp 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.

> What are the units?

Units are length^(-1), i.e. line length per unit area

> 2. How can I choose the sigma parameter in density.psp so the results align
> with the more direct, intuitive approach above?

They do not agree unless sigma is very small.

'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.

> 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.

Imagine replacing each line by a fine grid of points, each point having a 'weight'
equal to the length of line that it replaced. Then apply density.ppp to these points with weights.
The result is the density.psp output.

Adrian Baddeley
spatstat author

Prof Adrian Baddeley FAA
Curtin University