libgeos
BUG: Intersects predicate error for two valid polygons
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
The issue at this stage is not relevant to libgeos itself as you are using a third-party Python wrapper.
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
JTS shows that the first input is invalid, so the error may be in isValid rather than intersects.
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))
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.
Looks like it was 26af9f377 (January 2021) ... and the JTS commit that changed this behavior (isValid false -> true) is https://github.com/locationtech/jts/commit/b520430f425a41961a2da0f09731dcdfbffeb937
Attached XML test case failing in both GEOS and JTS.
If the hole is removed from the larger polygon, the test passes.
The failure is likely due to the known robustness issue in the spatial predicate algorithm (originally reported as https://github.com/locationtech/jts/issues/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.