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?

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