I'm trying to move from ArcPy for geoprocessing. After searching some posts, I found this solution using the GDAL/OGR Python bindings:
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:
Does anyone know what is causing the problem?