RGeo 交叉函数的问题

Question

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

例 1：

我有两个多边形，我相信它们共享一个边界但不共享任何内部空间（即它们接触但不重叠）：

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 应该以与与 poly_4（减法器）相交的方式完全相同的方式与intersection_a 相交。

除了它没有：

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 段线，这两者都不是。

Answer 1

简答

不幸的是，似乎您不能期望从Float 坐标touches?和overlaps?使用 Float 坐标时获得可靠和正确的输出。

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

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

Geometry#intersection比touches?and更可靠一些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埃）足以完全改变的结果，并且两个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?和其它谓词衍生自的二维延伸的九交模型（DE-9IM） 。它是一个矩阵，用于描述几何图形的内部、边界和外部之间的交集维度。

该矩阵src/operation/relate/RelateComputer.cpp在 GEOS 的第 72 行计算：

IntersectionMatrix* RelateComputer::computeIM()

该算法似乎需要精确点头，而在您的任何示例中都不是这种情况。

我能找到的所有 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具有圆坐标，在你的榜样，计算它们的差别造成的多边形与非精确坐标：x1的diff_a是-8243077.837796713。

几何#intersection

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

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

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

这种方法似乎比RelateComputer::computeIM.

使用 GeometryPrecisionReducer？

使用GeometryPrecisionReducer和 a PrecisionModel，可以为所有几何图形指定一个允许点的网格。

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

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