本文介绍了将绘图限制在r中的shapefile边界的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧! 问题描述 29岁程序员,3月因学历无情被辞! 我已经使用density.ppp分析了GPS点的数据集,以生成点强度的热图,如下所示: 然而,我希望图像是仅限于shapefile的边界,类似于以下内容: 第一张图片称为 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使用此代码调用法师: density.ppp 的结果有一个包含信息的矩阵(v),如果外部点在绘制之前,interst的多边形变为 NA ,那么它们不会绘制。这是一个这样做的例子: library(maptools) library(sp) library spatstat) xx IDvar =FIPSNO,proj4string = CRS(+ proj = longlat + ellps = clrk66)) x y tmp plot(tmp) plot(xx, add = TRUE) points(x,y) tmp2 proj4string = CRS(proj4string(xx))) tmp3 tmp $ v [is.na(tmp3 [[1] ])] 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 outmany 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]] ) ] <- NAplot(tmp)plot(xx, add=TRUE) 这篇关于将绘图限制在r中的shapefile边界的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持! 上岸,阿里云! 09-03 10:43