我已经使用 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")
其中 cases@coords 是每个兴趣点的 GPS 坐标,x 是提供地理单元轮廓的 shapefile。
第二张图片使用以下代码调用:
plot(x, axes=T, col=COLORS, border="White")
有谁知道如何做到这一点?也许 plot() 是不可能的,我需要另一个包。
顺便说一句,我计划做的下一步是将此图像覆盖在从 GoogleEarth 导入的 map 上。我也不确定该怎么做,但如果我解决了问题,我会发布答案
非常感谢
最佳答案
density.ppp
的结果有一个包含信息的矩阵 (v)密谋那么他们就不会密谋。这是一个这样做的例子:
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 边界,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/12319504/