Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

A bug caused by numerical accuracy #1871

Open
itmorn opened this issue Aug 19, 2023 · 2 comments
Open

A bug caused by numerical accuracy #1871

itmorn opened this issue Aug 19, 2023 · 2 comments

Comments

@itmorn
Copy link

itmorn commented Aug 19, 2023

from shapely import Polygon

poly1 = Polygon([[446, -207], [315, -75.99999999999999], [557, -76]])
poly2 = Polygon([[184, -207], [315, -76], [446, -207], [315, -338]])
intersection = poly1.intersection(poly2)
print(intersection.wkt)  # ouptput: POLYGON ((315 -75.99999999999999, 557 -76, 446 -207, 315 -75.99999999999999))
print(intersection.area)  # output: 15850.999999999998

poly1 = Polygon([[446, -207], [315, -76], [557, -76]])
poly2 = Polygon([[184, -207], [315, -76], [446, -207], [315, -338]])
intersection = poly1.intersection(poly2)
print(intersection.wkt)  # output: LINESTRING (446 -207, 315 -76)
print(intersection.area)  # output: 0.0
@mwtoews
Copy link
Member

mwtoews commented Aug 20, 2023

Older shapely versions (before OverlayNG) would have always given the second answer, with a LineString intersection output, without area. Newer shapely versions with OverlayNG unfortunately return an unexpected intersection in the first example, where the result geometry does not really intersect poly1. @dr-jts might be able to weigh in to say if this is expected or not:

A: POLYGON ((446 -207, 315 -75.99999999999999, 557 -76, 446 -207))
B: POLYGON ((184 -207, 315 -76, 446 -207, 315 -338, 184 -207))

If you are using Shapley 2+ with a modern GEOS (most are), then you can use the optional grid size to round the overlay coordinates, here using a very small grid size close to machine epsilon:

intersection = poly1.intersection(poly2, 1e-15)
print(intersection.wkt)  # output: LINESTRING (446 -207, 315 -76)

@dr-jts
Copy link

dr-jts commented Aug 20, 2023

This is probably caused by a known issue with OverlayNG, which creates incorrect results in some situations with almost-coincident linework.

An example is: libgeos/geos#942

So far there's no fix for this. The suggested workaround of using rounded coordinates is the best available right now.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

4 participants