[Matplotlib-users] importing basemap then shapely causes error

[David H. <dh...@gm...>]

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: 
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",  line334,  in  intersection
   return  geom_factory(self.impl['intersection'](self,  other))
File  "/sw/lib/python2.7/site-packages/shapely/topology.py",  line53,  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_toolkitsimport  basemap
from  shapelyimport  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|


Thanks again.  Please CC me in any replies.

-Dave