将绘图限制在r中的shapefile边界 [英] Restricting plot to shapefile boundaries in r

问题描述

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

然而，我希望图像是仅限于shapefile的边界，类似于以下内容：

第一张图片称为

 <$ 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，它提供了地理单元的轮廓。

第二i使用此代码调用法师：

  plot（x，axes = T，col = COLORS，border =White）
  

有人知道这可能会怎么做吗？也许这是不可能的情节（），我将需要另一个包。

另外，我打算做的下一步是将这个图像覆盖在地图上从GoogleEarth导入。我还不确定如何做到这一点，但是如果当我解决问题时会发布答案。

非常感谢
<解决方案 density.ppp 的结果有一个包含信息的矩阵（v），如果外部点在绘制之前，interst的多边形变为 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 < -  rorm（25,35,1）
 
 tmp < -  density（ppp（x，y，xrange = range（x），y​​range = 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）
  

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:

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

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

解决方案

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)

这篇关于将绘图限制在r中的shapefile边界的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！