I'm not sure if this is the right place to ask this, but here we go.
I have a python script that uses matplotlib's basemap and another part that uses shapely to do an intersection of 2 polygons. If basemap is imported before shapely and I run the intersection I get this exception:
intersect_poly = grid_poly.intersection(data_poly) File "/sw/lib/python2.7/site-packages/shapely/geometry/base.py", line 334, in intersection return geom_factory(self.impl['intersection'](self, other)) File "/sw/lib/python2.7/site-packages/shapely/topology.py", line 53, in __call__ "This operation produced a null geometry. Reason: unknown") shapely.geos.TopologicalError: This operation produced a null geometry. Reason: unknown
If I import shapely first everything works fine. I would assume this is because of some "funkiness" in the way they are accessing the GEOS library. I've checked that in both situations the same library file is loaded in shapely ("print shapely.geos._lgeos").
Does anyone have an idea as to why this is happening and if there is a right way of doing this? Does this happen for anyone else? In the mean time I can just make sure to import shapely first (not sure if that affects basemap yet). Otherwise maybe I'll skim through the basemap source and ask this question on one of the mailing lists. Any suggestions on which mailing list I should post this to? Thanks.
I'm using OSX(10.7) with a fink install that has "libgeos3.3.3-shlibs", "libgeos3.3.1-shlibs", "libgeos3.3.1", "libgeos3.3.0-shlibs", "libgeos3.3.0", and "shapely-py27 (1.2.16-1)" installed.
And here's a simple test script that reproduces the problem (flip the imports and it works):
from mpl_toolkits import basemap from shapely import geometry g_ring = [(-88.462425, 26.992203), (-57.847187, 26.992203), (-57.847187, 17.599869), (-88.462425, 17.599869), (-88.462425, 26.992203)] grid_g_ring = [(-123.044, 59.844000000000001), (-49.384999999999998, 57.289000000000001), (-65.090999999999994, 14.335000000000001), (-113.133, 16.369), (-123.044, 59.844000000000001)] data_poly = geometry.Polygon(g_ring) grid_poly = geometry.Polygon(grid_g_ring) print grid_poly.intersection(data_poly).area