interference shapely and OGR python 2.7 _gdal.so on Mac

[eric66]

Hello,

Running this program gives wrong results: False,True,True,True :

from shapely.geometry import Polygon
import ogr

a=ogr.CreateGeometryFromWkt("POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0))")
b1=ogr.CreateGeometryFromWkt("POLYGON ((4 3, 4 5, 6 5, 6 3, 4 3))")
b2=ogr.CreateGeometryFromWkt("POLYGON ((4 3, 4 5, 6 5, 6 4, 4 3))")

print b1.Intersects(a)
print a.Intersects(b1)

print b2.Intersects(a)
print a.Intersects(b2)

If I import ogr first, no problem.

After debbugging, it seems to come from _gdal.so (osgeo/init.py) which interferes with shapely.

Eric.

[Even Rouault]

It is a bit odd that shapely and the OGR python bindings could interact... I can't try on Mac, but at least I can say I don't reproduce on Linux :

$ python
Python 2.6.6 (r266:84292, Oct 29 2012, 23:09:05) 
[GCC 4.4.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> from shapely.geometry import Polygon
>>> import ogr
/home/even/gdal/svn/1.10/gdal/swig/python/build/lib.linux-x86_64-2.6/osgeo/gdal.py:80: DeprecationWarning: ogr.py was placed in a namespace, it is now available as osgeo.ogr
  DeprecationWarning)
>>> 
>>> a=ogr.CreateGeometryFromWkt("POLYGON ((0 0, 0 10, 10 10, 10 0, 0 0))")
>>> b1=ogr.CreateGeometryFromWkt("POLYGON ((4 3, 4 5, 6 5, 6 3, 4 3))")
>>> b2=ogr.CreateGeometryFromWkt("POLYGON ((4 3, 4 5, 6 5, 6 4, 4 3))")
>>> 
>>> print b1.Intersects(a)
True
>>> print a.Intersects(b1)
True
>>> 
>>> print b2.Intersects(a)
True
>>> print a.Intersects(b2)
True

[eric66]

It has the same result if I import ogr before shapely and use shapely, and an is_valid call crashes python.

$ python
Python 2.7.2 (default, Oct 11 2012, 20:14:37) 
[GCC 4.2.1 Compatible Apple Clang 4.0 (tags/Apple/clang-418.0.60)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from ogr import OFTInteger
>>> from shapely.geometry import Polygon
>>> a=Polygon([(0,0),(0,10),(10,10),(10,0)])
>>> b1=Polygon([(4,3),(4,5),(6,5),(6,3)])
>>> b2=Polygon([(4,3),(4,5),(6,5),(6,4)])
>>> print b1.intersects(a)
False
>>> print a.intersects(b1)
True
>>> print b2.intersects(a)
True
>>> print a.intersects(b2)
True
>>> print a.is_valid
Assertion failed: (dynamic_cast<DirectedEdgeStar*>(ees)), function linkResultDirectedEdges, file PlanarGraph.cpp, line 235.
Abort trap: 6
... crash

When I debug, I found that replacing my "from ogr import OFTInteger" with the instructions launched in ogr/init.py cause the problem :

import imp
fp, pathname, description = imp.find_module('_gdal', ["/Library/Frameworks/GDAL.framework/Versions/1.9/Python/2.7/site-packages/osgeo"])
imp.load_module('_gdal', fp, pathname, description)

Eric.

[Even Rouault]

Hum, that could be a packaging issue between the GDAL.framework and Shapely if they don't use the same libgeos binary. The crash in "(dynamic_cast<DirectedEdgeStar*>(ees)), function linkResultDirectedEdges, file PlanarGraph.cpp, line 235" is in GEOS ( ​http://svn.osgeo.org/geos/trunk/src/planargraph/PlanarGraph.cpp )...

You should likely report the issue to the packagers of GDAL and/or Shapely.

[Even Rouault]

Not clearly a GDAL problem. Closing