测试几何是否具有至少一个共同的边界点,而不管内部点(不像"gTouches")

Question

包中的gTouches函数rgeos测试"几何是否至少有一个共同的边界点,但没有内部点".我正在寻找一种方法来测试"几何具有至少一个共同的边界点",而没有与内部点相关的标准.

这是基本设置:我有两个大多数相互嵌入的shapefile.我想找到文件中的多边形,其中较小的区域位于较大区域的边界.这是一个描述我想要做的事情的图表:

plot(map2, col=NA, border='black', lwd=0.4)
plot(map1, col=NA, border='#666666', lwd=0.2, add=TRUE)

该图显示了纽约史坦顿岛的人口普查区块.其中一个较大区域的绿色突出显示了我想要识别的块.只有那些与较大区域(粗线)共享或跨越边界的人.不是位于较大区域中间的区块.我尝试gTouches(map2,map1, byid=TRUE)在rgeos包中使用with 和其他功能但没有成功.gTouches只返回FALSE可能是因为标准是"几何具有至少一个共同的边界点,但没有内部点".基本上,我正在寻找一种功能,测试"几何形状是否至少有一个共同的边界点",无论内部是什么.

后续问题是我是否可以获得相互边界的长度？

数据:您可以在此处和此处下载两张地图.两者都是rds文件,所以你可以像这样打开它们:

library('rgdal')
library('rgeos')
library('sp')
map1 = readRDS('map1.rds')
map2 = readRDS('map2.rds')

Answer 1

您可以使用组合gIntersects()(查找与学区任何部分相交的所有小多边形)和gContainsProperly()(查找完全包含在学区内并且不与学区边界相交的所有小多边形).然后简单地组合两个结果逻辑矩阵来识别您所追求的多边形.

## Identify polygons that intersect but aren't fully contained within the
## school district whose polygon is given by SD = map2[13,]
SD <- map2[13,]
ii <- gIntersects(SD, map1, byid=TRUE) &
      !gContainsProperly(SD, map1, byid=TRUE)
ii <- apply(ii, 1, any)  ## Handy construct if both layers contain >1 polygon

## Plot that area, to show that results are  correct
plot(SD, col=NA, border='black')                  ## Establish plotted area
plot(map1, col=NA, border='#666666', lwd=0.2, add=TRUE)
plot(map1[ii,], col="lightgreen", add=TRUE)
plot(SD, col=NA, border='black', lwd=2, add=TRUE) ## Put SD boundary on top 

编辑:

然而,那不太对劲.从上面的地图可以看出,学区的SW和SE内部的许多小多边形本应该被识别出来.像这样的结果在rgeos操作中经常发生,并且由于这对层(或GEOS引擎的中间表示)的微小错误注册而产生.

解决方案是gBuffer()在执行拓扑查询之前使用少量缓冲其中一个层.在这里,坐标以米为单位,一些试验和错误表明,20米的缓冲区足以解决问题:

## Expand every polygon in map1 by a 20-meter wide buffer 
map1_buff <- gBuffer(map1, byid=TRUE, width=20)

## and then use the buffered version of map1 in the topological queries
ii <- gIntersects(SD, map1_buff, byid=TRUE) &
      !gContainsProperly(SD, map1_buff, byid=TRUE)
ii <- apply(ii, 1, any)  ## Handy construct if both layers contain >1 polygon

## Plot that area, to show that results are  correct
plot(SD, col=NA, border='black')                  ## Establish plotted area
plot(map1, col=NA, border='#666666', lwd=0.2, add=TRUE)
plot(map1[ii,], col="lightgreen", add=TRUE)
plot(SD, col=NA, border='black', lwd=2, add=TRUE)

这仍然错过沿海岸的几个多边形,但在某些时候,一个完整的解决方案可能需要让一对地图在其细节水平上更好地匹配.如果缓冲区大小变得更大,则该分析将开始产生误报,例如,拾取学区西北角的一些真正的内部多边形.