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?

Improve this question
4
  • any reason that you don't want to use ogr2ogr as outlined in the link that you included in your post? 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. 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. 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. Commented Jan 12, 2016 at 19:54

1 Answer 1

Reset to default
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 :

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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.