Restricting plot to shapefile boundaries in r

3
votes

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

1
rgpsplotshapefile
what packages are you using? maptools and spatstat? – GSee
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
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
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

1 Answers

4
votes

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)