3

I have analysed a dataset of GPS points using density.ppp to generate a sort of heatmap of intensity of the points, as shown below:

enter image description here

However, I would like the image to be limited to the borders of the shapefile, similar to below:

enter image description here

The first image is called as 

x <- readShapePoly("dk.shp")
xlim<-c(min(912),max(920))
ylim<-c(min(8023),max(8030))
a<-ppp(cases@coords[,1], cases@coords[,2], xlim, ylim, unitname=c("km"))
plot(density.ppp(a, 0.1), col=COLORS)
plot(x, add=T, border="white")

where cases@coords are the GPS coordinates of each point of interest, and x is a shapefile which provides the outline for the geographical unit.

The second image is called using this code:

plot(x, axes=T, col=COLORS, border="White")

Does anyone know how this might be done? Perhaps it's not possible with plot() and I will need another package.

As an aside, the next step I plan to do will be to overlay this image over a map imported from GoogleEarth. I'm not yet sure how to do that either, but will post the answer if and when I work it out

many thanks

Jonny
  • 2,703
  • 2
  • 27
  • 35
  • what packages are you using? maptools and spatstat? – GSee Sep 07 '12 at 15:21
  • Well, the first `plot()` uses the `xlim` and `ylim` you provided above. If you want to manually adjust them use `plot(...,ylim=c(0,100),xlim=c(0,100)` You will find the xlim and ylim values from the shapefile calling `x@bbox` – jakob-r Sep 07 '12 at 16:27
  • If you have the borders of your shapes well defined, you should be able to overplot with `polygon` , defining your polygon between the shape borders and the plot borders, and filling that area with white or other background color. But I bet some `ggplot` or `maptools` expert will chime in with the correct answer :-) – Carl Witthoft Sep 07 '12 at 16:37
  • Thanks for your comments - the answer below solved the problem for me, and the packages used are listed at the start of Greg's answer. – Jonny Sep 08 '12 at 15:28

1 Answers1

4

The result of density.ppp has a matrix (v) that contains the information, if the points outside of the polygon of interst are changed to NA before it is plotted then they will not plot. Here is an example of doing this:

library(maptools)
library(sp)
library(spatstat)

xx <- readShapePoly(system.file("shapes/sids.shp", package="maptools")[1],
      IDvar="FIPSNO", proj4string=CRS("+proj=longlat +ellps=clrk66"))

x <- rnorm(25, -80, 2)
y <- rnorm(25, 35, 1 )

tmp <- density( ppp(x,y, xrange=range(x), yrange=range(y)) )
plot(tmp)
plot(xx, add=TRUE)
points(x,y)

tmp2 <- SpatialPoints( expand.grid( tmp$yrow, tmp$xcol )[,2:1],
    proj4string=CRS(proj4string(xx)) )

tmp3 <- over( tmp2, xx )

tmp$v[ is.na( tmp3[[1]] ) ] <- NA

plot(tmp)
plot(xx, add=TRUE)
Greg Snow
  • 48,497
  • 6
  • 83
  • 110