r - 测试几何图形是否至少有一个共同的边界点而不管内部点如何(不像 `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中的其他功能打包但没有成功。 gTouches只返回 FALSE可能是因为标准是“几何图形至少有一个共同的边界点，但没有内部点”。基本上，我正在寻找一种函数来测试“几何图形是否至少有一个共同的边界点”，而不管内部如何。

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

数据:您可以下载两张 map here和 here .两者都是 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 

编辑:

然而，这并不完全正确。从上面的 map 中可以看出，沿着学区西南和东南部内部的许多小多边形本应被识别出来。像这样的结果在 中经常发生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)

这仍然会错过沿海地区的几个多边形，但在某些时候，完整的解决方案可能需要获得一对在细节级别上更匹配的 map 。如果缓冲区大小变得更大，此分析将开始产生误报，例如，拾取学区 NW 角的一些真正内部多边形。