[prev in list] [next in list] [prev in thread] [next in thread] 

List:       r-sig-geo
Subject:    [R-sig-Geo] Calculating road (i.e. linear feature) density using spatstat::density.psp()
From:       Adrian Baddeley <adrian.baddeley () uwa ! edu ! au>
Date:       2015-03-19 0:42:35
Message-ID: CF5661163F77A44781208D9AC4FDEA72014010A2989A () IS-WIN-376 ! staffad ! uwa ! edu ! au
[Download RAW message or body]

Matt Strimas-Mackey <strimas@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
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


[prev in list] [next in list] [prev in thread] [next in thread]