Small difference in polygon coords results in point intersection instead of polygon intersection

[igheorghita]

Expected behavior and actual behavior.

I have two polygons which are extremely similar to each other and I intersect them both with a third polygon. I expect the intersection to be a large polygon, but in one case I end up with a point.

Steps to reproduce the problem.

from shapely import wkt
from shapely.ops import orient

zone = wkt.loads('POLYGON ((79.57006468414316 227.16715306136757, 43.44048377114246 -7.28462320799008, 127.74283932445364 -13.587090929970145, 154.99848808655608 219.60419326741248, 79.57006468414316 227.16715306136757))')

geom_0 = wkt.loads('POLYGON ((127.74283932445364 -13.587090929970145, 43.44048377114247 -7.28462320799008, 54.329847263553376 63.37854389007508, 54.61818554266568 63.33040101339976, 135.62637080265046 53.862076963892896, 127.74283932445364 -13.587090929970145))')
geom_0.intersection(zone).wkt # outputs 'POINT (127.74283932445364 -13.587090929970145)'

# yet orienting geom_0 first, gives a polygon
orient(geom_0).intersection(zone).wkt # 'POLYGON ((127.74283932445364 -13.587090929970145, 43.44048377114247 -7.28462320799008, 54.329847263553376 63.37854389007508, 54.61818554266568 63.33040101339976, 135.62637080265046 53.862076963892896, 127.74283932445364 -13.587090929970145))'

# this geometry is almost exactly the same as geom_0
geom_1 = wkt.loads('POLYGON ((127.74283932445363 -13.587091019377114, 43.44048377114246 -7.284623297862709, 54.32984726350768 63.37854379976558, 54.61818554261504 63.33040092309082, 135.62637080259915 53.862076873508386, 127.74283932445363 -13.587091019377114))')
geom_1.intersection(zone).wkt # outputs 'POLYGON ((54.32984726350768 63.37854379976558, 54.61818554261504 63.33040092309082, 135.6263707922279 53.862076874720586, 127.74283932445364 -13.587090929970145, 43.440483784834306 -7.284623209013686, 54.32984726350768 63.37854379976558))'

# another observation
geom_0.intersection(zone) # is a point
zone.intersection(geom_0) # is a point

# but
zone.buffer(0).intersection(geom_0) # is a polygon
geom_0.buffer(0).intersection(zone) # is a point

Operating system

Debian 5.19.11

Shapely version and provenance

Shapely 1.8.2 installed with pip

[mwtoews]

Thanks for the report. This is a familiar bug with GEOS. Another workaround is to use normalize() to get expected polygons:

geom_0.normalize().intersection(zone.normalize())
zone.normalize().intersection(geom_0.normalize())

These two results are now identical.

[igheorghita]

Thanks for the quick response. I hadn't heard about normalize().

I also ran into some other confusing behavior - getting ValueError: No Shapely geometry can be created from null value when trying to use unary_union. Small changes in the polygon coords make the difference between getting that error and getting a valid polygon.

buffer_0 = wkt.loads('POLYGON ((-8.060879443418989 -2.6638274935170196, -8.194659017192153 -33.71089267415941, -11.242630721688515 -33.69775918301853, -11.10885114791535 -2.6506940023761434, -8.060879443418989 -2.6638274935170196))')
buffer_1 = wkt.loads('POLYGON ((-84.63155622050482 -2.3338899473058037, -84.61842272936394 0.7140817571905578, -8.047745952278113 0.3841442109793418, -8.06087944341899 -2.6638274935170196, -84.63155622050482 -2.3338899473058037))')
buffer_2 = wkt.loads('POLYGON ((-84.63155622050482 -2.3338899473058037, -84.61842272936394 0.7140817571905576, -11.095717656774472 0.3972777021202178, -11.108851147915349 -2.6506940023761434, -84.63155622050482 -2.3338899473058037))')
buffers = [buffer_0, buffer_1, buffer_2]

new_buffer_0 = wkt.loads('POLYGON ((-8.060879443418989 -2.66382749351702, -8.194659017192153 -33.71089267415941, -11.242630721688515 -33.69775918301853, -11.10885114791535 -2.650694002376143, -8.060879443418989 -2.66382749351702))')
new_buffer_1 = wkt.loads('POLYGON ((-84.63155622050482 -2.333889947305804, -84.61842272936394 0.714081757190558, -8.047745952278113 0.384144210979342, -8.06087944341899 -2.66382749351702, -84.63155622050482 -2.333889947305804))')
new_buffer_2 = wkt.loads('POLYGON ((-84.63155622050482 -2.333889947305804, -84.61842272936394 0.714081757190558, -11.095717656774472 0.397277702120218, -11.108851147915349 -2.650694002376143, -84.63155622050482 -2.333889947305804))')
new_buffers = [new_buffer_0, new_buffer_1, new_buffer_2]

# unary_union(buffers) results in ValueError: No Shapely geometry can be created from null value
unary_union(new_buffers).wkt # 'POLYGON ((-8.194659017192153 -33.71089267415941, -11.242630721688515 -33.69775918301853, -11.10885114791535 -2.650694002376144, -84.63155622050482 -2.333889947305804, -84.61842272936394 0.714081757190558, -8.047745952278113 0.384144210979342, -8.06087944341899 -2.66382749351702, -8.060879443418989 -2.66382749351702, -8.194659017192153 -33.71089267415941))'

all([b.is_valid for b in buffers]) # True
all([b.is_valid for b in new_buffers]) # True

# unary_union([b.buffer(0) for b in buffers]) still gives the value error
# same with unary_union([b.normalize() for b in buffers])
unary_union([b.buffer(1e-6) for b in buffers]) # works
unary_union([orient(b) for b in buffers]) # also works

Any ideas on what might be able to consistently fix something like this in practice? And is this related to the same GEOS issue?

[gy-mate]

This is a familiar bug with GEOS.

@mwtoews Could you please share the link of the upstream bug? Thank you in advance!