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: Intersects predicate error for two valid polygons #766

Open
chuchulk opened this issue Dec 8, 2022 · 9 comments
Open

BUG: Intersects predicate error for two valid polygons #766

chuchulk opened this issue Dec 8, 2022 · 9 comments
Labels

Comments

@chuchulk
Copy link

chuchulk commented Dec 8, 2022

from shapely import wkt
poly1 = wkt.loads('POLYGON ((26639.240191093646 6039.3615818717535, 26639.240191093646 5889.361620883223,28000.000095100608 5889.362081553552, 28000.000095100608 6039.361620882992, 28700.000190214026039.361620882992, 28700.00019021402 5889.361822800367, 29899.538842431968 5889.362160452064,32465.59665091549 5889.362882757903, 32969.2837182586 -1313.697771558439, 31715.832811969216 -1489.87008918589, 31681.039836323587 -1242.3030298361555, 32279.3890331618 -1158.210534269224, 32237.63710287376 -861.1301136466199, 32682.89764107368 -802.0828534499739, 32247.445200905553 5439.292852892075, 31797.06861513178 5439.292852892075, 31797.06861513178 5639.36178850523, 29899.538849750803 5639.361268079038, 26167.69458275995 5639.3602445643955, 26379.03654594742 2617.0293071870683, 26778.062167926924 2644.9318977193907, 26792.01346261031 2445.419086759444, 26193.472956813417 2403.5650586598513, 25939.238114175267 6039.361685403233, 26639.240191093646 6039.3615818717535), (32682.89764107368 -802.0828534499738, 32682.89764107378 -802.0828534499669, 32247.445200905655 5439.292852892082, 32247.445200905553 5439.292852892075, 32682.89764107368 -802.0828534499738))')
poly2 = wkt.loads('POLYGON ((32450.100392347143 5889.362314133216, 32050.1049555691 5891.272957209961, 32100.021071878822 16341.272221116333, 32500.016508656867 16339.361578039587, 32450.100392347143 5889.362314133216))')
poly1.intersects(poly2)

When I perform the above intersects operation,both polygons are valid, but the line of poly1.intersects(poly2) raise TopologyException and tell me the input set may be is invalid.
I am doing with Shapely 2.0b2 / GEOS 3.5

@szymor
Copy link

szymor commented Dec 8, 2022

The issue at this stage is not relevant to libgeos itself as you are using a third-party Python wrapper.

@dbaston
Copy link
Member

dbaston commented Dec 8, 2022

I can confirm this in main using geosop, which shows both inputs are valid but throws a TopologyIntersection during intersection testing.

geosop -a 'POLYGON ((26639.240191093646 6039.3615818717535, 26639.240191093646 5889.361620883223,28000.000095100608 5889.362081553552, 28000.000095100608 6039.361620882992, 28700.00019021402 6039.361620882992, 28700.00019021402 5889.361822800367, 29899.538842431968 5889.362160452064,32465.59665091549 5889.362882757903, 32969.2837182586 -1313.697771558439, 31715.832811969216 -1489.87008918589, 31681.039836323587 -1242.3030298361555, 32279.3890331618 -1158.210534269224, 32237.63710287376 -861.1301136466199, 32682.89764107368 -802.0828534499739, 32247.445200905553 5439.292852892075, 31797.06861513178 5439.292852892075, 31797.06861513178 5639.36178850523, 29899.538849750803 5639.361268079038, 26167.69458275995 5639.3602445643955, 26379.03654594742 2617.0293071870683, 26778.062167926924 2644.9318977193907, 26792.01346261031 2445.419086759444, 26193.472956813417 2403.5650586598513, 25939.238114175267 6039.361685403233, 26639.240191093646 6039.3615818717535), (32682.89764107368 -802.0828534499738, 32682.89764107378 -802.0828534499669, 32247.445200905655 5439.292852892082, 32247.445200905553 5439.292852892075, 32682.89764107368 -802.0828534499738))' -b 'POLYGON ((32450.100392347143 5889.362314133216, 32050.1049555691 5891.272957209961, 32100.021071878822 16341.272221116333, 32500.016508656867 16339.361578039587, 32450.100392347143 5889.362314133216))' intersects    

@dbaston
Copy link
Member

dbaston commented Dec 8, 2022

JTS shows that the first input is invalid, so the error may be in isValid rather than intersects.

@dr-jts
Copy link
Contributor

dr-jts commented Dec 11, 2022

The format of the first WKT is invalid, since there's a missing space between two numbers:

28700.000190214026039.361620882992

If this is fixed, in JTS both polygons test as valid, and the OverlayNG intersection works, with result

POLYGON ((32449.982270224027 5889.362878362695, 32450.100395042435 5889.362878395946, 32450.100392347143 5889.362314133216, 32449.982270224027 5889.362878362695))

@dbaston
Copy link
Member

dbaston commented Dec 11, 2022

image

@dbaston
Copy link
Member

dbaston commented Dec 11, 2022

If I update JTS to the latest commit the polygon shows as valid. Unfortunately I wasn't smart enough to record what commit I was on before updating.

@dbaston
Copy link
Member

dbaston commented Dec 11, 2022

Looks like it was 26af9f377 (January 2021) ... and the JTS commit that changed this behavior (isValid false -> true) is locationtech/jts@b520430

@dbaston
Copy link
Member

dbaston commented Dec 11, 2022

Attached XML test case failing in both GEOS and JTS.

geos-issue-786.xml.txt

If the hole is removed from the larger polygon, the test passes.

@dr-jts dr-jts changed the title BUG: Intersection error with two valid polygons BUG: Intersects error with two valid polygons Dec 12, 2022
@dr-jts dr-jts changed the title BUG: Intersects error with two valid polygons BUG: Intersects predicate error with two valid polygons Dec 12, 2022
@dr-jts dr-jts added the Bug label Dec 12, 2022
@dr-jts dr-jts changed the title BUG: Intersects predicate error with two valid polygons BUG: Intersects predicate error for two valid polygons Dec 12, 2022
@dr-jts
Copy link
Contributor

dr-jts commented Dec 12, 2022

The failure is likely due to the known robustness issue in the spatial predicate algorithm (originally reported as locationtech/jts#396). Because a hole edge is very close to an edge of the polygon, these likely collapse together during noding, and thus create invalid topology during predicate evaluation.

As noted in #740, the hope is that an improved spatial predicate algorithm will avoid this situation.

Note that a detailed evaluation of the polygon with a hole confirms that it is valid. It might be worth determining what caused the valid algorithm to fail originally, however, and add a unit test for this.

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

4 participants