使用 R 中的简单特征库识别多边形的公共边界答案

【问题标题】：Identifying common borders of polygons using the Simple Features library from R使用 R 中的简单特征库识别多边形的公共边界
【发布时间】：2019-12-18 22:06:24
【问题描述】：

我正在尝试使用 sf 包中的 sf_intersection() 函数来识别两个不同多边形的公共边界。

我在我的数据中尝试了这个 simple approach，它来自一个 shapefile，但它并没有完全按照我的预期工作。

我的数据是来自this repository 的 shapefile 名称“区域”，这是我尝试过的：

library(sf)
library(ggplot2)

zones <- st_read('./Data/zones.shp')
zones$id <- seq(nrow(zones))
borders <- st_intersection(zones, zones)
borders <- borders[borders$id != borders$id.1, ]

ggplot() +
  geom_sf(data = zones, color='red', fill=NA) +
  geom_sf(data = borders, color = 'navy')

最终的情节产生了这个结果：

如果您仔细观察，您会注意到多边形内线的某些部分不是borders 中的线的一部分（它们是红色而不是蓝色）。我不知道为什么会这样。任何提示或建议将不胜感激。谢谢！

【问题讨论】：

正如@David_O 提到的，边界实际上并没有共享，因此不会相交。试试，例如：ggplot(data = zones) +geom_sf() + xlim(c(768824.7, 770092.5)) + ylim(c(2944198.6, 2944940.3)) 你会看到两个线串分叉的部分。
您可以在相交之前尝试st_buffer您的区域对象。
我用来修复多边形格子以修复非重叠边界或重叠条的最好方法是：github.com/tudelft3d/pprepair - 需要一些组装。

标签： r gis sf

【解决方案1】：

这是边界的局部不精确。对于大多数矢量数据格式，共享的 POLYGON 边界在每个邻居中重复。坐标上的细微差别并不需要太多的时间就可以使两个边界的交点不完整。

恐怕这不是解决方案。

这部分显示了问题的类型。视图窗口是：

> par('usr')
[1]  764968.2  765650.8 2945266.2 2945890.9

那条条子只有大约 3 米宽。

编辑：只是为了添加我对使用st_snap 的解决方案的尝试。这似乎在某些地方起到了作用，但并非始终如一。感觉它不像预期的那样工作。另外，只是要注意投影使用美国英尺作为单位，这让我感到困惑。

z1 <- st_geometry(zones[1,])
z2 <- st_geometry(zones[2,])

z1 <- st_cast(z1, 'LINESTRING')
z2 <- st_cast(z2, 'LINESTRING')

z1s <- st_snap(z1, z2, 1000)

border <- st_intersection(z1s, z2)

捕捉容差超出了顶部 - 区域之间的间隙似乎小于 10 英尺 - 但即使有这么大的容差，实际边界也缺少部分。更奇怪的是，结果有一个完全出乎意料的延伸，距离实际交叉口超过 6500 英尺。

【讨论】：

谢谢，这是我看到结果时想到的第一个想法。我会尝试寻找可能的解决方案来解决这个问题，如果我能找到它会发布它。
我玩过st_snap，我的直觉说它应该以给定的容差将一个多边形的顶点捕捉到另一个多边形，但它似乎并不像我预期的那样工作

【解决方案2】：

@David_O 确定了问题 - POLYGON 边界实际上并未触及整个共享边界，因此st_intersection 不会识别它们。

一种解决方法可能是在相交之前st_buffer您的区域对象，尽管这无疑是一种粗略的解决方法：

borders <- st_intersection(st_buffer(zones, 5), st_buffer(zones, 5))

borders <- borders[borders$id != borders$id.1, ]

ggplot() +
  geom_sf(data = zones, color='red', fill="transparent") +
  geom_sf(data = borders, color = 'navy')

【讨论】：

这是一个很好的解决方案。一个问题是它会给你一个围绕边界的 10m 宽的 POLYGON 而不是 LINESTRING。诚然，我怀疑您已经获得了 Polygon 或 GeometryCollection。如果这只是为了显示，可能很适合您，但如果您出于分析目的将边界作为几何图形，这可能是一个问题。
太棒了！事实上，使用st_buffer() 是解决我的问题的方法
嗯，是的 - 我认为它适用于简单的绘图，但不是很准确。我认为有一个st_snap 解决方案。让我调查一下..
我对@987654326@ 的尝试一直在返回 GEOMETRYCOLLECTION，其一端带有一个奇怪的多边形。不幸的是，我对 sf 对象和空间谓词还不够好，无法理解为什么会发生这种情况或如何解决它。