1

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.

5
  • shapely may be a good option, see: gis.stackexchange.com/questions/312907/… – Zoltan Apr 14 '20 at 13:19
  • I cannot use shapely as I cannot install additional modules to this system. In any case shapely uses GDAL/OGR so why not do the same? – Ingwe Apr 14 '20 at 13:37
  • Shapely is based on GEOS/JTS and its aim to manipulate planar geometric objects. It extends the functionality of GDAL/ORG. – Zoltan Apr 14 '20 at 15:11
  • Hi Zoltan. Thanks, semantics aside, I would really like to accomplish my original goal. DO you have any experience with GDAL/OGR? – Ingwe Apr 15 '20 at 8:47
  • Maybe this helps: gis.stackexchange.com/q/359867/2948, just make sure use intersection rather than st_intersects (the former returns the intersection, the former is a binary test). – Jose Apr 28 '20 at 23:35
2

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.

1
  • It does indeed run to completion without errors... but not with the expected output. Trying with manufactured data I get everything working exactly as expected as well, yet the shapefiles produced in an earlier part of the code, in particular, "OA_Temp.shp" do not give the desired result. The issue would appear to be that each polygon in this file has the same FID=0. This is a correct answer all the way however and I really appreciate the assistance. – Ingwe Apr 15 '20 at 16:46

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.