将绘图限制为 r 中的 shapefile 边界答案

【问题标题】：Restricting plot to shapefile boundaries in r将绘图限制为 r 中的 shapefile 边界
【发布时间】：2012-09-07 13:53:02
【问题描述】：

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

但是，我希望将图像限制在 shapefile 的边界内，类似于以下：

第一张图片被称为

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

其中 case@coords 是每个兴趣点的 GPS 坐标，x 是提供地理单位轮廓的 shapefile。

使用此代码调用第二个图像：

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

有谁知道如何做到这一点？也许用 plot() 是不可能的，我需要另一个包。

顺便说一句，我计划做的下一步是将这张图片覆盖在从 GoogleEarth 导入的地图上。我也不确定该怎么做，但如果我解决了，我会发布答案

非常感谢

【问题讨论】：

你用的是什么包？地图工具和spatstat？
嗯，第一个 plot() 使用您在上面提供的 xlim 和 ylim。如果您想手动调整它们，请使用plot(...,ylim=c(0,100),xlim=c(0,100) 您将从调用x@bbox的 shapefile 中找到 xlim 和 ylim 值@
如果您的形状边界定义明确，您应该能够使用 polygon 进行重叠绘图，在形状边界和绘图边界之间定义多边形，并用白色或其他填充该区域背景颜色。但我敢打赌，一些ggplot 或maptools 专家会给出正确答案:-)
感谢您的 cmets - 下面的答案为我解决了问题，并且使用的软件包列在 Greg 答案的开头。

标签： r gps plot shapefile

【解决方案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)

【讨论】：