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

BUG: inconsistency between prepared and unprepared polygon contains predicate #740

Open
brendan-ward opened this issue Nov 23, 2022 · 5 comments
Labels

Comments

@brendan-ward
Copy link
Contributor

I'm seeing cases where the contains predicate is true when the polygon is not prepared, but false when the polygon is prepared. The polygon is valid, so I don't think that is the issue. In this case, the line is created by clipping a larger line to the polygon. Example WKT files attached.

In contrast, the line is always within the polygon, regardless of whether or not the line is prepared.

I didn't see a way to use prepared predicates with geosop, but here is what I am doing with Shapely 2.0b2 / GEOS 3.11:

>>> import shapely
>>> poly = shapely.from_wkt(open("test_poly.wkt").read())
>>> line = shapely.from_wkt(open("test_line.wkt").read())
>>> shapely.contains(poly, line)  # calls GEOSContains_r
True
>>> shapely.prepare(poly)  # calls GEOSPrepare_r
>>> shapely.contains(poly, line)  # calls GEOSPreparedContains_r with prepared polygon
False
>>> shapely.within(line, poly)  # calls GEOSWithin_r
True
>>> shapely.prepare(line)
>>> shapely.within(line, poly)  # calls GEOSPreparedWithin_r
True

This might be a special case where the endpoint of the line is within the boundary of the polygon (fails ContainsProperly), which perhaps is causing one of the speedups in the prepared Contains predicate to fail when it shouldn't.

test_poly.wkt.zip
test_line.wkt.zip

@dr-jts
Copy link
Contributor

dr-jts commented Nov 23, 2022

In fact the test_line is NOT contained by the test_poly, since an endpoint of the line falls outside the polygon.
image
So the prepared contains result is correct, and the Relate result is wrong.

This is likely due to the issue described in locationtech/jts#396. Because the current Relate algorithm nodes the input before computing the DE-9IM topology, it can happen that due to robustness the noding causes the line endpoint and the computed node to be (incorrectly) identical. (This is not surprising, since the line endpoint is the computed intersection of a longer line with the polygon.)

Prepared Contains doesn't node the input, but uses (more) robust topological predicates. So it computes the correct answer.

Fixing this is going to require a new Relate algorithm, unfortunately.

@dr-jts dr-jts added the Bug label Nov 23, 2022
@brendan-ward
Copy link
Contributor Author

Thanks for investigating! In my case I was relying on (unprepared) contains giving me the right answer, so that I could attribute the line I clipped by the polygon as being (effectively) within the polygon. Good to know that in the future this could change and I may have to calculate the degree of intersection to get to the same result (i.e., >99% within is good enough).

@dr-jts
Copy link
Contributor

dr-jts commented Nov 24, 2022

I may have to calculate the degree of intersection to get to the same result (i.e., >99% within is good enough).

In the (currently experimental) RelateNG algorithm, I'm hoping to support specifying a tolerance distance for evaluating spatial predicates. That should handle this use case nicely.

@dr-jts
Copy link
Contributor

dr-jts commented Nov 24, 2022

Here's a minimal reproducer:

POLYGON ((1454700 -331500, 1455100 -330700, 1455466.6191038645 -331281.94727476506, 1455467.8182005754 -331293.26796732045, 1454700 -331500))
LINESTRING (1455389.376551584 -331255.3803222172, 1455467.2422460222 -331287.83037053316)

Contains = true
Prepared Contains = false

@dr-jts
Copy link
Contributor

dr-jts commented Nov 24, 2022

Great to see geosop being useful. It does support prepared predicates:

   A B  containsPrep > B - test if geometry A contains geometry B, using PreparedGeometry
   A B  containsProperlyPrep > B - test if geometry A properly contains geometry B using PreparedGeometry
   A B  coversPrep > B - test if geometry A covers geometry B using PreparedGeometry
   A B  intersectsPrep > B - test if geometry A intersects B using PreparedGeometry

Here's the command to evaluate the minimal repro above:

geosop -a "POLYGON ((1454700 -331500, 1455100 -330700, 1455466.6191038645 -331281.94727476506, 1455467.8182005754 -331293.26796732045, 1454700 -331500))" -b "LINESTRING (1455389.376551584 -331255.3803222172, 1455467.2422460222 -331287.83037053316)" containsPrep

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

No branches or pull requests

2 participants