search

libgeos / geos

Geometry Engine, Open Source
https://libgeos.org
GNU Lesser General Public License v2.1
1.2k stars 358 forks source link

BUG: Intersects predicate error for two valid polygons #766

Open chuchulk opened 1 year ago

chuchulk commented 1 year ago 
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 commented 1 year ago

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

dbaston commented 1 year ago

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 commented 1 year ago

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

dr-jts commented 1 year ago

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 commented 1 year ago

image

dbaston commented 1 year ago

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 commented 1 year ago

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

dbaston commented 1 year ago

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 commented 1 year ago

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.

spawn-guy commented 1 week ago

if i may, i'll add my 2 cents here.

geosop -a "0103000000010000000700000008B0FE7FFFEC10400BE7465F40CF494086BA757B61EF104044A905FFD0CD4940E999AD3361011140BA085A1A66CD49407257708100111140DB11BD956ACE49403BA7ADC89F0E11401F228EFFD9CF49407DF85B7F9EFC104034196DE444D0494008B0FE7FFFEC10400BE7465F40CF4940" -b "0103000000010000000700000061B6409764951040E0CE3EE7E5CF4940A881F769978A10402C106A9402CC4940D8698A58F5B01040D8DA8F31BEC949405221CEB526E21040277D6EF05CCB494009B0FE7FFFEC10400BE7465F40CF494003A368509BC610400E434EF384D1494061B6409764951040E0CE3EE7E5CF4940" intersects

results in

terminate called after throwing an instance of 'geos::util::TopologyException'
  what():  TopologyException: side location conflict at 4.2314434050749767 51.619151982899062. This can occur if the input geometry is invalid.
Aborted (core dumped)

on geosop - GEOS 3.10.6 on Amazon Linux 2023 that i deployed Yesterday and tested geosop Today

but now i reverted back to v3.13.0 (from the exact old rpm binaries) - AND THE PROBLEM HAVE DISAPPEARED!

geosop - GEOS 3.13.0
...
$ geosop -a "0103000000010000000700000008B0FE7FFFEC10400BE7465F40CF494086BA757B61EF104044A905FFD0CD4940E999AD3361011140BA085A1A66CD49407257708100111140DB11BD956ACE49403BA7ADC89F0E11401F228EFFD9CF49407DF85B7F9EFC104034196DE444D0494008B0FE7FFFEC10400BE7465F40CF4940" -b "0103000000010000000700000061B6409764951040E0CE3EE7E5CF4940A881F769978A10402C106A9402CC4940D8698A58F5B01040D8DA8F31BEC949405221CEB526E21040277D6EF05CCB494009B0FE7FFFEC10400BE7465F40CF494003A368509BC610400E434EF384D1494061B6409764951040E0CE3EE7E5CF4940" intersects
true

even more interestingly: the downgraded version(yesterday, via shapely+python) stopped throwing this errors after 2024-10-10 14:23:50,142 UTC! but today the old cli-geosop was failing...

what is happening?!? are these errors time-dependant? RNG dependant? was there a solar flare/eclipse/anything yesterday and now it's gone?

originally reported here: https://github.com/shapely/shapely/issues/2170 and was advised to escalate/"go deeper"(inception.jpg)