Kate Zinszer
2009-Aug-18 21:37 UTC
[R] kernel density estimation for univariate data using splancs
Hi,
I previously received help in extract data from a shapefile and now my question
is about kernel density estimation. My objective is to have 3 kernel density
plots; 2 for the each set of cases and the 3rd is the difference in kernel
densities between the 2 sets of cases. Previously, I used the density function
from the stats package, which worked but I wanted finer control of the
bandwidth. Therefore, I was advised to switch to the package splancs.
I've been trying to copy the example provided in the help file but cannot
produce the same results for my data. I've pasted my code below and I am
not sure how to define the grd1 (cellsize and cells.dim). Here is my code:
coords <- S at polygons[[1]]@Polygons[[1]]@coords
coord <- as.points(coords)
cases1.xy <- as.points(mergedis$jit.x.x, mergedis$pre.fsa.shp.Y2)
grd1 <- GridTopology(cellcentre.offset=c(50281627.2, 580082), cellsize=c(0.2,
0.2), cells.dim=c(75,100))
k1000 <- spkernel2d(cases1.xy, coord, h0=1000, grd1)
k5000 <- spkernel2d(cases1.xy, coord, h0=5000, grd1)
if (.sp_lt_0.9()) {
df <- AttributeList(list(k1000=k1000, k5000=k5000))
} else {
df <- data.frame(k1000=k1000, k5000=k5000)
}
kernels <- SpatialGridDataFrame(grd1, data=df)
spplot(kernels, checkEmptyRC=FALSE, col.regions=terrain.colors(16), cuts=15)
Using this code, I get NAs for all values of k1000 and k5000. If this did work,
can I then subtract the spplots from one another to have the difference between
the two sets of cases?
Any advice would be greatly appreciated!
Kate
. You need to check the @hole
slot of the ... at Polygons[[1]] object if you care whether it's a hole
or an island!
Barry
