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
  • 1
    Basemap and Shapely each use GEOS, but in a different ways. Basemap has a C extension that links GEOS in the usual manner, but Shapely accesses GEOS functions via FFI (libffi and ctypes). I scanned the basemap code and don't see why the packages shouldn't be compatible. It's worth asking on the basemap list. – sgillies Feb 8 '13 at 16:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.