GeometryCollection#union with slightly overlapping input polygons excludes some polygons in its result

[grimsa]

(Not sure if this is a bug or just me using JTS in an unsupported way)

The need: union ~5-20 polygons that may be (slightly) overlapping

My current implementation uses new GeometryCollection(polygons).union() and it has worked reliably so far.

However, I found one case where the union does not fail, yet the produced polygon is missing large areas.

Input:

Union result (red) over input (blue):

Versions used:

• JTS 1.18.2 (with OverlayNG enabled)
• Java 16

Test case:

    @Test
    void union_slightlyOverlappingPolygons_areasShouldNotBeLost() throws ParseException {
        var gf = new GeometryFactory();
        Polygon initialPolygon = (Polygon) new WKTReader().read(
                "POLYGON ((-17.055671609725188 0.0838585273938056, -17.219533690879903 8.580102931174729, -12.78759751645274 4.315883851167196, -12.707594043193543 0.1677170531469111, -17.055671609725188 0.0838585273938056))"
        );
        Function<double[][], Polygon> createPolygon = points -> gf.createPolygon(
                IntStream.range(0, points[0].length)
                        .mapToObj(index -> new CoordinateXY(points[0][index], points[1][index]))
                        .toArray(CoordinateXY[]::new)
        );
        var p1 = createPolygon.apply(new double[][]{
                {-16.15907384118054, -16.28767148661218, -13.719459554848422, -13.639456293556348, -16.15907384118054},
                {1.0157206673651495, 7.6835050482212575, 5.212481574525698, 1.0643148816523527, 1.0157206673651495}
        });
        var p2 = createPolygon.apply(new double[][]{
                {-12.725226307666448, -13.639456293556348, -13.719459768757051, -12.78759751645274, -12.725226307666448},
                {1.0819470369308437, 1.0643148816523527, 5.212481780339414, 4.315883851167196, 1.0819470369308437}
        });
        var p3 = createPolygon.apply(new double[][]{
                {-13.621824029083443, -16.14144162383529, -16.15907384118054, -13.639456293556348, -13.621824029083443},
                {0.15008489786842003, 0.10149068267229658, 1.0157206673651493, 1.0643148816523527, 0.15008489786842003}
        });
        var p4 = createPolygon.apply(new double[][]{
                {-17.219533690879903, -16.28767148661218, -16.15907384118054, -17.073303827070436, -17.219533690879903},
                {8.580102931174729, 7.683505048221257, 1.0157206673651495, 0.9980885120866584, 8.580102931174729}
        });
        var p5 = createPolygon.apply(new double[][]{
                {-12.707594043193543, -13.621824029083443, -13.639456293556348, -12.725226307666448, -12.707594043193543},
                {0.1677170531469111, 0.15008489786842005, 1.0643148816523527, 1.0819470369308437, 0.1677170531469111}
        });
        var p6 = createPolygon.apply(new double[][]{
                {-17.055671609725188, -17.073303827070436, -16.15907384118054, -16.14144162383529, -17.055671609725188},
                {0.0838585273938056, 0.9980885120866584, 1.0157206673651493, 0.1014906826722966, 0.0838585273938056}
        });
        Polygon[] polygons = {
                p1, p2, p3, p4, p5, p6
        };

        var gcUnion = gf.createGeometryCollection(polygons).union();
        var mpUnion = gf.createMultiPolygon(polygons).union();
        var indUnion = Arrays.stream(polygons)
                .map(Geometry.class::cast)
                .reduce(Geometry::union)
                .orElseThrow();

        assertAll(
                () -> assertEquals(initialPolygon.getArea(), gcUnion.getArea(), 0.0001), // Fail
                () -> assertEquals(initialPolygon.getArea(), mpUnion.getArea(), 0.0001), // Fail
                () -> assertEquals(initialPolygon.getArea(), indUnion.getArea(), 0.0001) // Pass
        );
    }

Some notes:

• Both GeometryCollection#union and MultiPolygon#union produce the same incorrect result (displayed above)
• Unioning polygons one by one produces the correct result, but I assume it to be more expensive than doing a single union operation on the whole set of polygons (if that is supported)
• Removing jts.overlay=ng flag makes the test pass - all 3 ways produce the same correct result

So I guess my questions are:

1. Should GeometryCollection/MultiPolygon produce a correct union for this slightly overlapping set of polygons? (i.e. is it supported)
2. If not, what would be the preferred way to union ~5-20 simple polygons - would it be doing it one by one (like in the test), or is there some better way?

[dr-jts]

Thanks for the detailed issue! I don't see anything wrong with the way you are using JTS (and it's interesting seeing JTS being used with modern functional Java code).

This does look like some sort of bug. I can run your test code and get the failure. The odd thing is that if I load the input GeometryCollection into the TestBuilder and run it it works correctly!

GEOMETRYCOLLECTION (POLYGON ((-16.15907384118054 1.0157206673651495, -16.28767148661218 7.6835050482212575, -13.719459554848422 5.212481574525698, -13.639456293556348 1.0643148816523527, -16.15907384118054 1.0157206673651495)), POLYGON ((-12.725226307666448 1.0819470369308437, -13.639456293556348 1.0643148816523527, -13.719459768757051 5.212481780339414, -12.78759751645274 4.315883851167196, -12.725226307666448 1.0819470369308437)), POLYGON ((-13.621824029083443 0.15008489786842, -16.14144162383529 0.1014906826722966, -16.15907384118054 1.0157206673651493, -13.639456293556348 1.0643148816523527, -13.621824029083443 0.15008489786842)), POLYGON ((-17.219533690879903 8.580102931174729, -16.28767148661218 7.683505048221257, -16.15907384118054 1.0157206673651495, -17.073303827070436 0.9980885120866584, -17.219533690879903 8.580102931174729)), POLYGON ((-12.707594043193543 0.1677170531469111, -13.621824029083443 0.1500848978684201, -13.639456293556348 1.0643148816523527, -12.725226307666448 1.0819470369308437, -12.707594043193543 0.1677170531469111)), POLYGON ((-17.055671609725188 0.0838585273938056, -17.073303827070436 0.9980885120866584, -16.15907384118054 1.0157206673651493, -16.14144162383529 0.1014906826722966, -17.055671609725188 0.0838585273938056)))

I'll keep digging and try and find out what's going on. It's possible it's something to do with numeric precision that only shows up when the data is loaded as Java double-precision literals - although would be surprising, given the extensive robustness heuristics in OverlayNG.

[dr-jts]

I think the problem is related to the fact that you are specifying 17 digits of precision in the input ordinates. For instance, this one. 0.15008489786842003. If I truncate this to 16 digits of precision (0.15008489786842) then the test case works.

I still think the various levels of robust geometric code in OverlayNG should handle this - but obviously it doesn't! Perhaps it's not surprising - 17 digits is right at the limit of double-precision floating-point representation.

I'll keep digging, but in the meantime you can round off your input numbers a bit.

[dr-jts]

The same problem occurs when computing the union of p3 and p5, so that simplifies testing.

[dr-jts]

A test case is here: https://github.com/dr-jts/jts/tree/fix-union-highprec

[grimsa]

Ran into a similar issue of union of polygons with very high precision coordinates producing an empty result. Works fine when serialized to WKT and loaded into test builder. Did not check against latest JTS snapshot.

Test Case:

    @Test
    void geometryCollectionUnion() {
        var gf = new GeometryFactory();
        var p1 = gf.createPolygon(new Coordinate[]{
                new Coordinate(-0.49524584515584785, 2.254321259525603),
                new Coordinate(-0.4518674095819571, 1.6850769474821623),
                new Coordinate(0.11737696012368323, 1.7284553797131785),
                new Coordinate(-3.208506829071115, -2.146129227235706),
                new Coordinate(-3.821129633016156, -1.6202633492159437),
                new Coordinate(-0.49524584515584785, 2.254321259525603)
        });

        var p2 = gf.createPolygon(new Coordinate[]{
                new Coordinate(-0.49524584515584785, 2.254321259525603, 6.986031855278883),
                new Coordinate(-0.4518674095819571, 1.6850769474821623, 6.986031855334046),
                new Coordinate(0.11737696012368318, 1.7284553797131785, 6.986031855278883),
                new Coordinate(-0.49524584515584785, 2.254321259525603, 6.986031855278883)
        });

        // Passes
        assertFalse(p2.union(p1).isEmpty(), "p2 union p1");
        // Passes
        assertFalse(p1.union(p2).isEmpty(), "p1 union p2");
        // Fails in JTS 1.18.2 with OverlayNG
        assertFalse(gf.createGeometryCollection(new Polygon[]{p1, p2}).union().isEmpty(), "geometry collection union");
    }