在R中限制绘图范围为shapefile边界

Question

我使用density.ppp分析了一个GPS点数据集，生成了一种点密度的热图，如下所示：

然而，我希望图片被限制在形状文件的边界内，类似于下面这张图片：

第一张图片被称为

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")

其中cases@coords是每个感兴趣点的GPS坐标，x是提供地理单元轮廓线的shapefile文件。

使用以下代码调用第二张图像：

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

有人知道这个怎么做吗？也许使用plot()不可能，我需要另一个包。

顺便说一下，我计划下一步要将这个图像叠加在从GoogleEarth导入的地图上。我还不确定如何做到这一点，但如果我弄清楚了，我会发布答案。

非常感谢

- Jonny

---

你正在使用哪些程序包？是maptools和spatstat吗？ - GSee

首先，第一个 plot() 使用了你上面提供的 xlim 和 ylim。如果你想手动调整它们，可以使用 plot(...,ylim=c(0,100),xlim=c(0,100)。你可以通过调用 x@bbox 来获取 shapefile 的 xlim 和 ylim 值。 - jakob-r

如果您的形状边界定义得很好，那么您应该能够使用“polygon”进行重叠绘制，在形状边界和绘图边界之间定义多边形，并用白色或其他背景颜色填充该区域。但我敢打赌一些“ggplot”或“maptools”专家会发表正确的答案 :-) - Carl Witthoft

感谢您的评论 - 下面的答案解决了我的问题，并且使用的软件包在Greg的答案开头列出。 - Jonny

Answer 1

density.ppp 的结果包含一个矩阵(v)，它包含信息，如果在绘制图形之前将兴趣多边形外的点更改为NA，则它们将不会被绘制。以下是一个示例：

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