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_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?

  • any reason that you don't want to use ogr2ogr as outlined in the link that you included in your post?
    – dmci
    Commented Jan 12, 2016 at 18:51
  • 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.
    – Ernesto561
    Commented Jan 12, 2016 at 19:25
  • I wonder how you can use GDAL python bindings if you do not have GDAL installed. Not advising to use os.system.call, though.
    – user30184
    Commented Jan 12, 2016 at 19:49
  • 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.
    – Ernesto561
    Commented Jan 12, 2016 at 19:54

1 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 :

  • Shapefile iteself has its limit. So relying on GDAL APIs seems to be a long term strategy.
    – Chang
    Commented Dec 12, 2020 at 19:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.