importing basemap then shapely causes error

I’ve asked this question on GIS stack exchange site, but thought it
would be good to post here too. The SE question is here:
Thanks again. Please CC me in any replies.
-Dave

···

http://gis.stackexchange.com/questions/50394/importing-matplotlib-basemap-and-shapely
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.

    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. The current basemap version in fink is
1.0.2.

    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, 2
6.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