4

the gTouches function in the rgeos package tests whether "geometries have at least one boundary point in common, but no interior points". I am looking for a way to test whether "geometries have at least one boundary point in common" without the criteria related to interior points.

Here is the basic setup: I have two shapefiles that are mostly embedded in each other. I want to find the polygons in the file with the smaller areas that are at the border of the larger areas. Here is a graph to describe what I am trying to do:

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

enter image description here

The figure shows census blocks in Staten Island, NY. The green highlighting in one of the larger areas illustrates the blocks I want to identify. Only those that share or cross a border of the larger areas (thick lines). Not the blocks that are in the middle of the larger areas. I tried to do this with with gTouches(map2,map1, byid=TRUE) and other function in the rgeos package but without success. gTouches only returns FALSE probably because the criteria is that "geometries have at least one boundary point in common, but no interior points". Basically, I am looking for a function that tests whether "geometries have at least one boundary point in common" regardless of the interior.

A follow-up question is whether I can get the length of the mutual border?

Data: You can download the two maps here and here. Both are rds files so you can open them like this:

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

1 回答 1

5

您可以使用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

在此处输入图像描述

编辑 ：

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

在此处输入图像描述

这仍然错过了沿海岸的几个多边形，但在某些时候，一个完整的解决方案可能必须涉及获得一对在细节水平上匹配得更好的地图。如果缓冲区大小变得更大，这种分析将开始产生误报，例如，在学区西北角的一些真正的内部多边形。

于 2013-11-14T18:42:55.970 回答