Problem intersecting shapefiles with OGR in Python

Question

I'm trying to move from ArcPy for geoprocessing. After searching some posts, I found this solution using the GDAL/OGR Python bindings:

Intersect Shapefiles using Shapely

My code is the following:

import os
act_dir = os.path.dirname(os.path.abspath(__file__))
print act_dir
from osgeo import ogr
ogr.UseExceptions()
ogr_ds = ogr.Open(os.path.join(act_dir,r'data'), True) 
SQL = """\
 SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, A.*, B.*
 FROM Central_America_24h A, municipios_simp_WGS84_dptos_n B
 WHERE ST_Intersects(A.geometry, B.geometry);
"""
layer = ogr_ds.ExecuteSQL(SQL, dialect='SQLITE')
# copy result back to datasource as a new shapefile
layer2 = ogr_ds.CopyLayer(layer, 'mylayer')
#save, close
layer = layer2 = ogr_ds = None

I have a folder called 'data', which contains the file 'Central_America_24h.shp' (point shapefile) and 'municipios_simp_WGS84_dptos_n.shp' (polygon shapefile). I want to intersect the points with the polygons and assign them their attributes. It seems the problem is when I try to save the result and I get the following error:

Traceback (most recent call last):
 File "D:\pyAnaconda\puntos_de_calor_test\puntos_calor_tes2t.py", line 39, in <module>
 layer2 = ogr_ds.CopyLayer(layer, 'mylayer')
 File "C:\Anaconda2\lib\site-packages\osgeo\ogr.py", line 815, in CopyLayer
 return _ogr.DataSource_CopyLayer(self, *args, **kwargs)
ValueError: Received a NULL pointer.

I have checked the documentation, but I don't know what I'm doing wrong:

OGRLayer * OGRDataSource::CopyLayer

Does anyone know what is causing the problem?

• dmci: any reason that you don't want to use ogr2ogr as outlined in the link that you included in your post?
• Ernesto561: Thanks, I am trying to avoid using an os.system.call and don´t want to install OGR. I have already installed the GDAL/OGR python bindings and would like to use them for this purpose.
• user30184: I wonder how you can use GDAL python bindings if you do not have GDAL installed. Not advising to use os.system.call, though.
• Ernesto561: Sorry, I don't know much about GDAL. Yes, I have GDAL installed, and found the ogr2ogr executable in the site_packages folder of my Python installation. But I would like to use the GDAL python bindings if possible.

Answer 1

I do not understand why beginners try to start with the GDAL/OGR Python bindings (not very "Pythonic" and difficult) when other easier alternatives are available.

With your script, you need to know osgeo.ogr and the SQL dialect of SQLite. The solution proposed by Mike T is powerful but not "classic" and performs only the intersection of shapefiles.

What you are trying to do is a Spatial Join (point in Polygon) and not a simple intersection and you can use :

• pyshp (shapefile) (pure Python module) as in Point in Polygon or Point in Polygon 2: Walking the line
• Fiona (other Python binding of OGR with Python dictionaries) and Shapely as in More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (1) and many other examples in GIS SE
• GeoPandas as in More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (2)
• If you want to use OGR, consult the Python GDAL/OGR Cookbook and many examples in GIS SE

• Chang: Shapefile iteself has its limit. So relying on GDAL APIs seems to be a long term strategy.