1

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?

  • any reason that you don't want to use ogr2ogr as outlined in the link that you included in your post? – dmci Jan 12 '16 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 Jan 12 '16 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 Jan 12 '16 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 Jan 12 '16 at 19:54
7

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 :

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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