I'm getting strange/incorrect results from RGeo polygon intersection functions (Ruby 2.3.0, RGeo 0.5.3)
EXAMPLE 1:
I have two polygons which I believe share a boundary but do not share any internal space (i.e. they touch but do not 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)
When we check the intersection between them, it returns a line, as expected for geometries sharing only a boundary:
poly_1.intersection poly_2
=> #<RGeo::Geos::CAPILineStringImpl:0x3fc0249af168 "LINESTRING (-8224757.546680832 4960523.476563589, -8225598.074380083 4961210.51680879)">
However, when running the following checks, we get the opposite of what's expected:
poly_1.overlaps? poly_2
=> true
poly_1.touches? poly_2
=> false
EXAMPLE 2:
We take two legitimately overlapping polygons:
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)
And calculate the differences in both directions:
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))">
Now we check the remainder polygons against their subtractors:
diff_b.touches? poly_3
=> true
diff_b.overlaps? poly_3
=> false
This is good.
diff_a.touches? poly_4
=> false
diff_a.overlaps? poly_4
=> true
This is IMPOSSIBLE- There's no way the remainder of a polygon subtracted from another should be able to overlap that polygon!
And why is it only happening in one direction?
EXAMPLE 3:
The weirdness continues. Let's now get the intersection of poly_3 and 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))">
Now since this is what should have been subtracted from poly_3 to give us diff_a, diff_a should therefore intersect with intersection_a in exactly the same way as that it intersects with poly_4 (the subtractor).
Except it doesn't:
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))">
Worse, neither of those two intersection results make sense. It should be a single, 2-segmented line, which neither of those are.
Short answer
Unfortunately, it seems that you cannot expect reliable and correct output from
touches?
andoverlaps?
when using Float coordinates.It doesn't depend on RGeo or GEOS versions (or for that matter, JTS, the project on which GEOS is based).
If you need information about the position of two polygons, you could use
Geometry#distance
and check that it is smaller than a given epsilon.Geometry#intersection
is a bit more reliable thantouches?
andoverlaps?
, but it isn't guaranteed to work with every example either.Is it a precision problem?
Theory
touches?
is by definition very sensitive : polygons have to share a point or a line on their border, but shouldn't have any common interior points.touches?
is false.touches?
is false.Sensitivity analysis
How sensitive is it exactly?
Let's consider two polygons, one just above the other :
We could move the lower border of the upper polygon by an epsilon, and see which influence it has :
It outputs :
Those are well behaved polygons, with
A difference of
1E-15
of a degree (about 1 Ångström) is enough to completely change the results, and bothoverlaps?
andtouches?
switch betweentrue
andfalse
.The intersection is either empty, a line or a polygon, but at least the results seem consistent between
intersection
,overlaps?
andtouches?
.In your case, having more complex polygons with tilted sides makes the problem harder. When calculating intersection of two segments, it's hard to maintain 1-Ångström accuracy!
Is it a problem with RGeo?
RGeo uses GEOS, which is a C++ port of JTS (Java Topology Suite).
In order to check that the problem isn't specific to RGeo or GEOS, I calculated Example 1 with JTS 1.14 :
It outputs :
The intersection is the exact same as in your example,
intersects?
andtouches?
output the same wrong results as RGeo.Why are results inconsistent?
Why do
intersection
andtouches?
return inconsistend results?DE-9IM
touches?
,intersects?
,overlaps?
and other predicates are derived from the Dimensionally Extended nine-Intersection Model (DE-9IM). It's a matrix describing the dimensions of intersection between interiors, boundaries and exteriors of geometries.This matrix is calculated on line 72 of
src/operation/relate/RelateComputer.cpp
in GEOS :The algorithm seems to require exact noding, which isn't the case in any of your examples.
All the tests that I could find for JTS used integers coordinates, even a test called "complex polygons touching" :
There's not a single test case where predicates are calculated for floating-point coordinates.
Note that even though
wkt_3
andwkt_4
have rounded coordinates in your example, calculating their differences create polygons with inexact coordinates :x1
ofdiff_a
is-8243077.837796713
.Geometry#intersection
Geometry#intersection
is calculated on line 670 ofsrc/operation/overlay/OverlayOp.cpp
in GEOS :Comments in the code seems to indicate that exact noding isn't required, and there are multiple if statements to check if the results appear correct.
This method seems to be more robust than
RelateComputer::computeIM
.Using GeometryPrecisionReducer?
With
GeometryPrecisionReducer
and aPrecisionModel
, it becomes possible to specify a grid of allowable points for all Geometries.GeometryPrecisionReducer
is implemented in GEOS, but isn't available in RGeo. Tests in JTS with your first example have shown that it doesn't solve your problem anyway : inexact coordinates are snapped to the closest point on the grid. Every point is moved a bit to the north/south/east or west, which changes the slopes of every side.It also changes the boundaries and their intersections. Depending on the
PrecisionModel
, the intersection for your first example could be empty, a line or a polygon.