Issue in Union of two intersecting polygons with non-zero areas returns empty

[stalhabukhari]

I encountered few rare cases where the ops.unary_union and Polygon.union methods misbehave, returning an empty polygon.

Expected behavior and actual behavior.

In the scenario shown below, there are two polygons with non-zero areas and calculating their union results in an empty polygon with zero area. Instead, their union should have roughly been the larger polygon (since one is a subset of the other) and hence should have a non-zero area.

from shapely.geometry import Polygon
from shapely.ops import unary_union

check1 = Polygon([[10.38343007464573, 945.7998193082476], [1682.0, 1335.7904894634128], [1682.0, 914.7369958740137],
                  [103.54509436006265, 546.4810826847133], [10.38343007464573, 945.7998193082476]])
check2 = Polygon(
    [[103.54509436006265, 546.4810826847133], [10.383430074645716, 945.7998193082476], [1682.0, 1335.7904894634128],
     [1775.161664285417, 936.4717528398785], [103.54509436006265, 546.4810826847133]])

print(check1.area, check2.area, check1.intersects(check2), unary_union([check1, check2]).area,
      check1.union(check2).area)

Output:

684226.9746010096 703839.9967089986 True 0.0 0.0

Plot:

Operating system

Ubuntu 20.04.3 (focal)

Shapely version and provenance

1.8.1.post1 installed from PyPI using pip

[jorisvandenbossche]

@stalhabukhari thanks for the report!

I can confirm this issue using GEOS 3.10.2 (didn't test master). The union gives an empty polygon, while for example the intersection seems to work correctly:

>>> check1.intersection(check2).wkt
'POLYGON ((1682 1335.7904894634128, 1682 914.7369958740137, 103.54509436006265 546.4810826847133, 10.38343007464573 945.7998193082476, 1682 1335.7904894634128))'

cc @dr-jts

[jorisvandenbossche]

A direct geos reproducer using geosop:

# intersection
$ geosop -a "POLYGON ((10.38343007464573 945.7998193082476, 1682 1335.7904894634128, 1682 914.7369958740137, 103.54509436006265 546.4810826847133, 10.38343007464573 945.7998193082476))" -b "POLYGON ((103.54509436006265 546.4810826847133, 10.383430074645716 945.7998193082476, 1682 1335.7904894634128, 1775.161664285417 936.4717528398785, 103.54509436006265 546.4810826847133))" -f wkt intersection
POLYGON ((1682 1335.7904894634128, 1682 914.7369958740137, 103.54509436006265 546.4810826847133, 10.38343007464573 945.7998193082476, 1682 1335.7904894634128))
# union
$ geosop -a "POLYGON ((10.38343007464573 945.7998193082476, 1682 1335.7904894634128, 1682 914.7369958740137, 103.54509436006265 546.4810826847133, 10.38343007464573 945.7998193082476))" -b "POLYGON ((103.54509436006265 546.4810826847133, 10.383430074645716 945.7998193082476, 1682 1335.7904894634128, 1775.161664285417 936.4717528398785, 103.54509436006265 546.4810826847133))" -f wkt union
POLYGON EMPTY

[stalhabukhari]

@jorisvandenbossche thanks for the prompt response!

I downgraded shapely to 1.7.1 and that seems to work fine for my cases for now.

[mwtoews]

This is fixed in GEOS main with union result POLYGON ((1682 1335.7904894634128, 1682 914.7369958740137, 103.54509436006265 546.4810826847133, 10.38343007464573 945.7998193082476, 1682 1335.7904894634128)). But I see the issue with POLYGON EMPTY from union in GEOS 10.2.

Commands for testing:

geosop -a 'POLYGON ((10.38343007464573 945.7998193082476, 1682 1335.7904894634128, 1682 914.7369958740137, 103.54509436006265 546.4810826847133, 10.38343007464573 945.7998193082476))' -b 'POLYGON ((103.54509436006265 546.4810826847133, 10.383430074645716 945.7998193082476, 1682 1335.7904894634128, 1775.161664285417 936.4717528398785, 103.54509436006265 546.4810826847133))' intersection -f wkt

[jorisvandenbossche]

That's good to hear! Do you know if the fix is also backported to 3.10.x?

[dr-jts]

This was fixed in GEOS main by libgeos/geos@5ea36b1.

It hasn't been backported to the 3.10 series yet.

[mwtoews]

This is a (not obvious) duplicate of #1216 ; it would be nice to see this backported, if feasible.

[dr-jts]

Now ported to GEOS 3.10 series: libgeos/geos@27fed83

[theroggy]

This is a (not obvious) duplicate of #1216 ; it would be nice to see this backported, if feasible.

Because the #1216 fix was shipped a while back in version 3.11, I had a look here as well and with geos 3.11 the test script above gives the folowing output, which looks good I think:

684226.9746010096 703839.9967089986 True 703839.9967089986 703839.9967089986

So, I suppose this issue can be closed as well...