[GIS] GDAL/OGR Clip Shapefile

gdalogrpython

I cannot find any solid info on this and its a pretty rudimentry GIS process.

I have 2 shapefiles.
One shapefiles is a single large polygon.
The other is a larger multipolygon.

I need to clip the multipolygon to the extents of the smaller polygon. None of the multipolygons are completely outside of the clipping polygon.

I can do this quite easily in Arc/QGIS/MsPaint/Whatever but this has to be done in python and using OGR.

Code Snippet:

import os, sys, csv
from osgeo import ogr

## Input
inDataSource = driver.Open('OAs_Temp.shp', 0)
inLayer = inDataSource.GetLayer()

## Clip
inClipSource = driver.Open('_convexhull.shp', 0)
inClipLayer = inClipSource.GetLayer()

## Clipped Shapefile... Maybe??? 
outDataSource = driver.CreateDataSource('FINAL.shp')
outLayer = outDataSource.CreateLayer('FINAL.shp', geom_type=ogr.wkbMultiPolygon)


featureDefn = outLayer.GetLayerDefn()
outFeature = ogr.Feature(featureDefn)
ogr.Layer.Clip(inLayer, inClipLayer, outLayer)
outFeature.SetGeometry(outLayer)
outLayer.CreateFeature(outFeature)

Issue is with outFeature.SetGeometry(outLayer)

Traceback:
TypeError: in method 'Feature_SetGeometry', argument 2 of type 'OGRGeometryShadow *'

The problem is that I just can't figure out how OGR is doing things and what it wants. I am finding that cookbook thing pretty lacklustre and incomplete and the documentation is not very clear on how any of this works. It is some of the most convoluted open source libraries I have ever worked with and I regret taking this project on.

This is my last option. As always.

Best Answer

I've tested your code, and here is a working version. I use Python 3.6.9 on Ubuntu 18.04, gdal.VersionInfo() '2020300'.

from osgeo import ogr

## Input
driverName = "ESRI Shapefile"
driver = ogr.GetDriverByName(driverName)
inDataSource = driver.Open('your_input.shp', 0)
inLayer = inDataSource.GetLayer()

print(inLayer.GetFeatureCount())
## Clip
inClipSource = driver.Open('your_clip.shp', 0)
inClipLayer = inClipSource.GetLayer()
print(inClipLayer.GetFeatureCount())

## Clipped Shapefile... Maybe??? 
outDataSource = driver.CreateDataSource('your_output.shp')
outLayer = outDataSource.CreateLayer('FINAL', geom_type=ogr.wkbMultiPolygon)

ogr.Layer.Clip(inLayer, inClipLayer, outLayer)
print(outLayer.GetFeatureCount())
inDataSource.Destroy()
inClipSource.Destroy()
outDataSource.Destroy()

Layer.Clip adds features (geometry and attributes) to the output layer, nothing to do with them.

Related Question