Fedor, It's a classic floating-point precision problem. To my eye, all your coordinates have 2 decimal places precision, and watch what happens when we trim down the WKT:
import shapely.wkt #wkt1 = 'POLYGON ((-660.0000000000000000 1.5500000000000000, -640.0000000000000000 6.1600000000000001, 420.0000000000000000 -5.2000000000000002, 430.0000000000000000 -5.4000000000000004, 440.0000000000000000 -5.2999999999999998, 550.0000000000000000 -6.0999999999999996, 560.0000000000000000 -6.2000000000000002, 570.0000000000000000 -6.4000000000000004, 2000.0000000000000000 -8.8000000000000007, 2000.0000000000000000 -10.0999999999999996, -660.0000000000000000 -10.0999999999999996, -660.0000000000000000 1.5500000000000000))' #wkt2 = 'LINESTRING (-660.0000000000000000 -5.3000000000000007, 2000.0000000000000000 -5.3000000000000007)' wkt1 = 'POLYGON ((-660.00 1.55, -640.00 6.16, 420.00 -5.20, 430.00 -5.40, 440.00 -5.30, 550.00 -6.10, 560.00 -6.20, 570.00 -6.40, 2000.00 -8.80, 2000.00 -10.10, -660.00 -10.10, -660.00 1.55))' wkt2 = 'LINESTRING (-660.00 -5.30, 2000.00 -5.30)' g1 = shapely.wkt.loads(wkt1) g2 = shapely.wkt.loads(wkt2) assert g1.is_valid assert g2.is_valid print g1.intersection(g2) I get: GEOMETRYCOLLECTION (POINT (440.0000000000000000 -5.2999999999999998), LINESTRING (-660.0000000000000000 -5.2999999999999998, 425.0000000000000000 -5.2999999999999998)) The right answer GEOMETRYCOLLECTION (POINT (440.0 -5.3), LINESTRING (-660.0 -5.3, 425.0 -5.3)) is computed, but the WKT writer can't perfectly represent that answer. The less you use WKT as an intermediary format in your code, the better. And if you have to use it, trim your values to their actual precision. Cheers, On Tue, Jul 13, 2010 at 11:16 AM, Emmanuel Lambert <emmanuel.lamb...@intec.ugent.be> wrote: > Hi, > > >From now and then I also experience simimar issues and it is indeed a > numerical error. When using the same polygons in a larger geometry for > example, the problem is often not occurring. > I had some improvement by using Geos 3.2.2, but the problem is still not > completely solved. Hopefully somebody with good knowledge of the code > can look into it... I might look into it myself after summer but > currently don't have time. > > wbr > Emmanuel > > > > can On Tue, 2010-07-13 at 11:10 +0200, Fedor Baart wrote: >> Hi all, >> >> I am running into a strange issue with shapely. >> I can't seem to be able to intersect two geometries. The code I use is: >> >> >> import shapely.wkt >> wkt1 = 'POLYGON ((-660.0000000000000000 1.5500000000000000, >> -640.0000000000000000 6.1600000000000001, 420.0000000000000000 >> -5.2000000000000002, 430.0000000000000000 -5.4000000000000004, >> 440.0000000000000000 -5.2999999999999998, 550.0000000000000000 >> -6.0999999999999996, 560.0000000000000000 -6.2000000000000002, >> 570.0000000000000000 -6.4000000000000004, 2000.0000000000000000 >> -8.8000000000000007, 2000.0000000000000000 -10.0999999999999996, >> -660.0000000000000000 -10.0999999999999996, -660.0000000000000000 >> 1.5500000000000000))' >> wkt2 = 'LINESTRING (-660.0000000000000000 -5.3000000000000007, >> 2000.0000000000000000 -5.3000000000000007)' >> g1 = shapely.wkt.loads(wkt1) >> g2 = shapely.wkt.loads(wkt2) >> >> assert g1.is_valid >> assert g2.is_valid >> >> g1.intersection(g2) >> >> >> I get the following error: >> >> File "<stdin>", line 1, in <module> >> File "/var/folders/V5/V5I4l+zDH4uYZ15MPmQ4ek+++TI/-Tmp-/ >> python-19066a9y.py", line 11, in <module> >> g1.intersection(g2) >> File "/Users/fedorbaart/Library/Python/2.6/site-packages/shapely/ >> geometry/base.py", line 318, in intersection >> return geom_factory(self.impl['intersection'](self, other)) >> File "/Users/fedorbaart/Library/Python/2.6/site-packages/shapely/ >> topology.py", line 55, in __call__ >> "This operation produced a null geometry. Reason: unknown") >> shapely.geos.TopologicalError: This operation produced a null >> geometry. Reason: unknown >> >> I think there might be some floating point error because if I remove >> point with y=-5.2999999999999998, which is very close to >> -5.3000000000000007, it works. >> I created the same script using geos+c++, that works fine. >> >> Can someone suggest a fix for this, or a way around it? >> >> Tested using geos 3.22, gcc 4.4, python 2.6, Shapely 1.2.1 >> >> >> Thanks and kind regards, >> >> Fedor Baart >> >> >> >> _______________________________________________ >> Community mailing list >> Community@lists.gispython.org >> http://lists.gispython.org/mailman/listinfo/community > > > _______________________________________________ > Community mailing list > Community@lists.gispython.org > http://lists.gispython.org/mailman/listinfo/community > -- Sean Gillies Programmer Institute for the Study of the Ancient World New York University _______________________________________________ Community mailing list Community@lists.gispython.org http://lists.gispython.org/mailman/listinfo/community
