【问题标题】：Problems with RGeo intersection functionsRGeo 交集函数的问题
【发布时间】：2016-12-27 17:18:42
【问题描述】：

我从 RGeo 多边形相交函数（Ruby 2.3.0、RGeo 0.5.3）得到奇怪/不正确的结果

示例 1：

我有两个多边形，我认为它们共享一个边界但不共享任何内部空间（即它们touch 但不overlap）：

wkt_1 = "POLYGON ((-8226874.27782158 4962626.76394919, -8223358.174520462 4961756.817075645, -8223358.174520462 4960289.557693501, -8224471.369428394 4960289.557693501, -8226874.27782158 4962253.674727506, -8226874.27782158 4962626.76394919))"
wkt_2 = "POLYGON ((-8224757.546680832 4960523.476563589, -8225269.1002275925 4959296.105368667, -8226993.791361805 4959219.668340384, -8226420.900079966 4961883.087589158, -8224757.546680832 4960523.476563589))"
poly_1 = RGeo::Geos.factory(:srid => 3857).parse_wkt(wkt_1)
poly_2 = RGeo::Geos.factory(:srid => 3857).parse_wkt(wkt_2)

当我们检查它们之间的交点时，它会返回一条线，这与仅共享边界的几何图形一样：

poly_1.intersection poly_2
=> #<RGeo::Geos::CAPILineStringImpl:0x3fc0249af168 "LINESTRING (-8224757.546680832 4960523.476563589, -8225598.074380083 4961210.51680879)">

但是，在运行以下检查时，我们得到的结果与预期相反：

poly_1.overlaps? poly_2
=> true
poly_1.touches? poly_2
=> false

示例 2：

我们采用两个合法重叠的多边形：

wkt_3 = "POLYGON ((-8243237.0 4970203.0, -8243237.0 4968735.0, -8242123.0 4968735.0, -8242123.0 4970203.0, -8243237.0 4970203.0))"
wkt_4 = "POLYGON ((-8244765.0 4966076.0, -8244765.0 4964608.0, -8243652.0 4964608.0, -8243652.0 4966076.0, -8242680.0 4969362.0, -8244765.0 4966076.0))"
poly_3 = RGeo::Geos.factory(:srid => 3857).parse_wkt(wkt_3)
poly_4 = RGeo::Geos.factory(:srid => 3857).parse_wkt(wkt_4)

并计算两个方向的差异：

diff_a = poly_3.difference poly_4
=> #<RGeo::Geos::CAPIPolygonImpl:0x3fe3fca26028 "POLYGON ((-8243077.837796713 4968735.0, -8243237.0 4968735.0, -8243237.0 4970203.0, -8242123.0 4970203.0, -8242123.0 4968735.0, -8242865.466828971 4968735.0, -8242680.0 4969362.0, -8243077.837796713 4968735.0))">
diff_b = poly_4.difference poly_3
=> #<RGeo::Geos::CAPIPolygonImpl:0x3fe3fd1dda28 "POLYGON ((-8242865.466828971 4968735.0, -8243652.0 4966076.0, -8243652.0 4964608.0, -8244765.0 4964608.0, -8244765.0 4966076.0, -8243077.837796713 4968735.0, -8242865.466828971 4968735.0))">

现在我们根据减法器检查剩余的多边形：

diff_b.touches? poly_3
=> true
diff_b.overlaps? poly_3
=> false

这很好。

diff_a.touches? poly_4
=> false
diff_a.overlaps? poly_4
=> true

这是不可能的 - 从另一个多边形中减去的多边形的剩余部分不可能与该多边形重叠！

为什么它只发生在一个方向上？

示例 3：

怪事还在继续。现在让我们得到 poly_3 和 poly_4 的交集

intersection_a = poly_3.intersection poly_4
=> #<RGeo::Geos::CAPIPolygonImpl:0x3fd1084ece88 "POLYGON ((-8242865.724520766 4968734.582337743, -8243078.32501591 4968734.582337743, -8242680.062418439 4969362.301390027, -8242865.724520766 4968734.582337743))">

现在既然这是应该从 poly_3 中减去得到 diff_a 的值，因此 diff_a 应该与 intersection_a 相交的方式与它与 poly_4（减法器）相交的方式完全相同。

除非它没有：

diff_a.touches? poly_4
=> false
diff_a.touches? intersection_a
=> true
diff_a.intersection poly_4
=> #<RGeo::Geos::CAPILineStringImpl:0x3fe3f98fb854 "LINESTRING (-8242680.0 4969362.0, -8243077.837796713 4968735.0)">
diff_a.intersection intersection_a
=> #<RGeo::Geos::CAPIMultiLineStringImpl:0x3fe3fca157b4 "MULTILINESTRING ((-8242865.466828971 4968735.0, -8242680.0 4969362.0), (-8242680.0 4969362.0, -8243077.837796713 4968735.0))">

更糟糕的是，这两个交集结果都没有意义。它应该是一条 2 段的线，但两者都不是。

【问题讨论】：

很有趣，但我不知道这可能来自哪里。克隆github.com/rgeo/rgeo 并使用git grep -A5 "overlaps?" 没有显示太多信息。你试过用另一个 srid 吗？
@EricDuminil - 刚刚用 srid 4326 尝试过 - 结果相同。此外，我认为如果交点是一条线，那么无论坐标空间如何，什么是接触而不是重叠的交点？
是的。您可能必须计算交点并手动检查结果是否为空、（多）点、（多）线串或（多）多边形。
@EricDuminil 有点违背了使用库的目的
当然，但是如果库失败...

标签： ruby geometry geospatial rgeo

【解决方案1】：

简答

不幸的是，seems 在使用 Float 坐标时，您无法期望来自 touches? 和 overlaps? 的可靠且正确的输出。

它不依赖于 RGeo 或 GEOS 版本（或者就此而言，JTS，GEOS 所基于的项目）。

如果您需要有关两个多边形位置的信息，可以使用 Geometry#distance 并检查它是否小于给定的 epsilon。

Geometry#intersection 比touches? 和overlaps? 更可靠一点，但也不能保证它适用于每个示例。

是精度问题吗？

理论

touches? 根据定义非常敏感：多边形必须在其边界上共享一个点或一条线，但不应有任何共同的内部点。

如果它们相距太远，则它们没有任何共同点，touches? 为假。
如果距离太近，它们的交点是多边形，touches? 为假。

敏感性分析

它到底有多敏感？

让我们考虑两个多边形，一个在另一个之上：

POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))
POLYGON((0 1, 0 2, 1 2, 1 1, 0 1))

我们可以将上多边形的下边界移动一个 epsilon，看看它有什么影响：

require 'rgeo'

epsilon = 1E-15

deltas = [-epsilon, 0, epsilon]

deltas.each do |delta|
  puts "--------------------------------"
  puts "Delta : #{delta}\n\n"
  simple1 = 'POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))'
  simple2 = "POLYGON((0 #{1+delta}, 0 2, 1 2, 1 #{1+delta}, 0 #{1+delta}))"

  puts "Polygon #1 : #{simple1}\n"
  puts "Polygon #2 : #{simple2}\n\n"

  poly_1 = RGeo::Geos.factory(:srid => 4326).parse_wkt(simple1)
  poly_2 = RGeo::Geos.factory(:srid => 4326).parse_wkt(simple2)

  puts "Intersection : #{poly_1.intersection poly_2}"
  puts

  puts "Distance between polygons :"
  puts format('%.30f°', poly_1.distance(poly_2))
  puts

  puts "Overlaps? : #{poly_1.overlaps? poly_2}"
  puts "Touches? : #{poly_1.touches? poly_2}"
end

它输出：

--------------------------------
Delta : -1.0e-15

Polygon #1 : POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))
Polygon #2 : POLYGON((0 0.999999999999999, 0 2, 1 2, 1 0.999999999999999, 0 0.999999999999999))

Intersection : POLYGON ((0.0 0.999999999999999, 0.0 1.0, 1.0 1.0, 1.0 0.999999999999999, 0.0 0.999999999999999))

Distance between polygons :
0.000000000000000000000000000000°

Overlaps? : true
Touches? : false
--------------------------------
Delta : 0

Polygon #1 : POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))
Polygon #2 : POLYGON((0 1, 0 2, 1 2, 1 1, 0 1))

Intersection : LINESTRING (0.0 1.0, 1.0 1.0)

Distance between polygons :
0.000000000000000000000000000000°

Overlaps? : false
Touches? : true
--------------------------------
Delta : 1.0e-15

Polygon #1 : POLYGON((0 0, 0 1, 1 1, 1 0, 0 0))
Polygon #2 : POLYGON((0 1.000000000000001, 0 2, 1 2, 1 1.000000000000001, 0 1.000000000000001))

Intersection : GEOMETRYCOLLECTION EMPTY

Distance between polygons :
0.000000000000001110223024625157°

Overlaps? : false
Touches? : false

这些是表现良好的多边形，具有

平行于轴的边
大小相同
整个长度的共享边

1E-15 一个度数的差异（大约1 Ångström）足以完全改变结果，overlaps? 和touches? 在true 和false 之间切换。

交叉点要么是空的，要么是一条线，要么是一个多边形，但至少intersection、overlaps? 和touches? 之间的结果似乎是一致的。

在您的情况下，具有倾斜边的更复杂的多边形会使问题变得更加困难。在计算两条线段的交点时，很难保持 1-Ångström 的精度！

RGeo 有问题吗？

RGeo 使用GEOS，这是JTS (Java Topology Suite) 的C++ 端口。

为了检查问题不是特定于 RGeo 或 GEOS，我计算了使用 JTS 1.14 的示例 1：

WKTReader wktReader = new WKTReader();
String wkt1 = "POLYGON ((-8226874.27782158 4962626.76394919, -8223358.174520462 4961756.817075645, -8223358.174520462 4960289.557693501, -8224471.369428394 4960289.557693501, -8226874.27782158 4962253.674727506, -8226874.27782158 4962626.76394919))";
String wkt2 = "POLYGON ((-8224757.546680832 4960523.476563589, -8225269.1002275925 4959296.105368667, -8226993.791361805 4959219.668340384, -8226420.900079966 4961883.087589158, -8224757.546680832 4960523.476563589))";
Polygon poly1 = (Polygon) wktReader.read(wkt1);
Polygon poly2 = (Polygon) wktReader.read(wkt2);
System.out.println("Intersection : " + poly1.intersection(poly2));
System.out.println("Overlaps?    : " + poly1.overlaps(poly2));
System.out.println("Intersects?  : " + poly1.intersects(poly2));
System.out.println("Touches?     : " + poly1.touches(poly2));
showMatrixWith(poly1, poly2);

它输出：

Intersection : LINESTRING (-8224757.546680832 4960523.476563589, -8225598.074380083 4961210.51680879)
Overlaps?    : true
Intersects?  : true
Touches?     : false

212
101
212

交叉点与您的示例完全相同，intersects? 和 touches? 输出与 RGeo 相同的错误结果。

为什么结果不一致？

为什么intersection 和touches? 返回不一致的结果？

DE-9IM

touches?、intersects?、overlaps? 和其他谓词都是从the Dimensionally Extended nine-Intersection Model (DE-9IM) 派生的。它是一个矩阵，描述几何内部、边界和外部之间的交集维度。

这个矩阵是在 GEOS 中 src/operation/relate/RelateComputer.cpp 的第 72 行计算得出的：

IntersectionMatrix* RelateComputer::computeIM()

该算法似乎是require exact noding，在您的任何示例中都不是这种情况。

我能为 JTS 找到的所有测试都使用整数坐标，甚至是一个名为“复杂多边形接触”的测试：

  # line 477 in jts-1.14/testxml/general/TestFunctionAA.xml
  <desc>mAmA - complex polygons touching</desc>
  <a>
    MULTIPOLYGON(
      (
        (100 200, 100 180, 120 180, 120 200, 100 200)),
      (
        (60 240, 60 140, 220 140, 220 160, 160 160, 160 180, 200 180, 200 200, 160 200,
        160 220, 220 220, 220 240, 60 240),
        (80 220, 80 160, 140 160, 140 220, 80 220)),
      (
        (280 220, 240 180, 260 160, 300 200, 280 220)))
  </a>
  <b>
    MULTIPOLYGON(
      (
        (80 220, 80 160, 140 160, 140 220, 80 220),
        (100 200, 100 180, 120 180, 120 200, 100 200)),
      (
        (220 240, 220 220, 160 220, 160 200, 220 200, 220 180, 160 180, 160 160, 220 160,
        220 140, 320 140, 320 240, 220 240),
        (240 220, 240 160, 300 160, 300 220, 240 220)))
  </b>

没有一个测试用例为浮点坐标计算谓词。

请注意，即使在您的示例中 wkt_3 和 wkt_4 具有舍入坐标，计算它们的差异会创建坐标不精确的多边形：diff_a 的 x1 是 -8243077.837796713。

几何#intersection

Geometry#intersection 在 GEOS 中 src/operation/overlay/OverlayOp.cpp 的第 670 行计算：

void OverlayOp::computeOverlay(OverlayOp::OpCode opCode)

代码中的注释似乎表明不需要精确的点头，并且有多个 if 语句来检查结果是否正确。

这个方法似乎比RelateComputer::computeIM更健壮。

使用 GeometryPrecisionReducer？

使用GeometryPrecisionReducer 和PrecisionModel，可以为所有几何体指定允许点的网格。

GeometryPrecisionReducer 在 GEOS 中实现，但在 RGeo 中不可用。在 JTS 中使用您的第一个示例进行的测试表明，它无论如何都不能解决您的问题：不精确的坐标被捕捉到网格上最近的点。每个点都向北/南/东或西移动一点，这会改变每一侧的坡度。

它还改变了边界和它们的交点。根据PrecisionModel，您的第一个示例的交点可能是空的、一条线或一个多边形。

【讨论】：

我的荣幸。我喜欢 Ruby、GIS 和大赏 :)

【解决方案2】：

在第一个示例中，您为两个多边形提供了一些浮点数。是什么让你认为它们接触而不重叠？你怎么知道？在一般情况下，您不能用绝对术语（不使用公差）说两条线段是否完全重叠。您可以提供具有圆角坐标的线条，在这种情况下，可以为其设置参数。

使用浮点数是实数空间离散化的一种形式。当使用不同的算法或方法来计算相同的结果时，这会导致结果不一致。例如，考虑二维空间中三条线的交点，它们应该在同一点相交。现在，每次使用不同的线对独立计算交点 3 次。您每次都可能得到相似但不平等的答案。

我没有使用 RGeo，但是在计算相交几何形状时有没有办法调整公差？尝试降低容差值（使该阈值更小）。这将允许该函数在不合并彼此太近的点的情况下计算几何形状。

【讨论】：

我以为你在做某事，所以我用四舍五入的值重新编写了示例 2-3（请参阅上面的编辑），但不幸的是我得到了相同的结果。