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 > [email protected] > http://lists.gispython.org/mailman/listinfo/community _______________________________________________ Community mailing list [email protected] http://lists.gispython.org/mailman/listinfo/community
