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
------------------------------------------------------------------------------
Free Next-Gen Firewall Hardware Offer
Buy your Sophos next-gen firewall before the end March 2013
and get the hardware for free! Learn more.
http://p.sf.net/sfu/sophos-d2d-feb
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
