I've been using spatstat to analyse point patterns in R. Here's an example to put this in context:

library(maptools) # to convert sp class to ppp

Some point data:

x <- rnorm(20)
y <- rnorm(20)
p <- SpatialPoints(coords = matrix(c(x, y), ncol = 2))

The challenge is to find the kernel density at each point, not in a continuous field aggregated as cells in a raster image, as calculated by density from spatstat:

p <- as.ppp(p) d <- density.ppp(p, sigma = 0.3) plot(d)

enter image description here

From this one can identify clusters and other useful things, but it does not provide imediate insight into the point that has the highest dot density. So this is a two-part question:

  1. How to extract the coordinates of the pixel with the highest density in d?
  2. Is there a way to calculate the kernel density only for the 20 points in p?
  • p=as.ppp(p) fails (R 3.1.1, latest spatstat) "Can't interpret X as a point pattern"... – Spacedman Sep 25 '14 at 9:42
  • Ah yes @Spacedman, forgot to add maptools dependency. Should work now. – RobinLovelace Sep 25 '14 at 9:54

I'd convert the output to a raster object. Then:

x <- runif(20)
y <- runif(20)
p <- SpatialPoints(coords = matrix(c(x, y), ncol = 2))

Then compute densities:

pp = ppp(x,y) # all points in a (0,1) default window
d <- density.ppp(pp, sigma = 0.1)
dp <- density.ppp(pp, sigma = 0.1, at="points")

That's Q2 answered! For Q1 I turn to the raster package:

dr = raster(d)
xyFromCell(dr, which.max(dr))
             x          y
[1,] 0.1523438 0.00390625

Note this is on slightly different data than you because I did it with data on a (0,1) square. Now I've got maptools your max point comes out here:

> dr = raster(d)
> xyFromCell(dr, which.max(dr))
           x         y
[1,] 1.33514 0.3392474
  • There is an extra = in pp = = ppp(x,y). – user32309 Sep 25 '14 at 10:01
  • 1
    Fixed. Must always do a final run through. Thanks! – Spacedman Sep 25 '14 at 10:02
  • This is awesome! How could anyone have guessed about the at argument in density.ppp? Just shows it's important to read the documentation as it's in there. That will save lots of computational time. Your answer will save lots of brain time, thanks! – RobinLovelace Sep 25 '14 at 10:11
  • I didn't guess - I read the help for density.ppp! – Spacedman Sep 25 '14 at 10:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.