Erasing features using Python GDAL/OGR?

Question

I have a big round polygon feature (5km diameter) and a small round polygon feature (3km diameter) that I would like to clip with Python GDAL/OGR to obtain the "donut". Once this is done, I will apply the same process for 25.000 features, that is why I need to do this with GDAL/OGR and not editing manually in QGIS/ArcGIS. It seems that the Clip operation works at the level of layers. Therefore, what I did was creating a layer for my 5km circle, a layer for my 3km circle and an output layer required by the prototype of the function.

Apparently the operation is finishing correctly, because I get a .shp file with 6KB inside, but when I open this shapefile either with QGIS or ArcGIS it is completely blank. I thought it was the projection, but the .prj file says it is set to EPSG:28992 (Amersfoort_RD_New), which is what I need.

The code I wrote is the following:

outDriver = ogr.GetDriverByName("ESRI Shapefile")    
outDataSource = outDriver.CreateDataSource(pathout)
outLayer = outDataSource.CreateLayer("Big", srs = outSpatialRef, geom_type=ogr.wkbPolygon)


outDriver2 = ogr.GetDriverByName("ESRI Shapefile")    
outDataSource2 = outDriver2.CreateDataSource(pathout)
outLayer2 = outDataSource2.CreateLayer("Small", srs = outSpatialRef, geom_type=ogr.wkbPolygon)


outDriver3 = ogr.GetDriverByName("ESRI Shapefile")    
outDataSource3 = outDriver3.CreateDataSource(pathout)
outLayer3 = outDataSource3.CreateLayer("Clip", srs = outSpatialRef, geom_type=ogr.wkbPolygon)

outLayer.CreateFeature(mybigfeature)
outLayer2.CreateFeature(mysmallfeature)
something = outLayer.Clip(outLayer2, outLayer3)
print something # returns zero

So, I do not understand how the operation seems to finish correctly, writes 6KB in my disk and then I see nothing on QGIS. I have attached an image, showing the features, but not the clipped polygon.

I guess there is something missing: memory flush?

insert a geometry for the output (how if there is no return of Clip)?

close the files?

Hope someone can spot the problem!

---

Following the advice of Benjamin, I installed Shapely and Fiona libraries to solve this problem. Shapely is used "for manipulation and analysis of planar geometric objects" and Fiona does the reading and writing of shapefiles in a Pythonic and very neat way.

So after playing around with these libraries, I came up with the following code:

from shapely.geometry import shape, mapping
import fiona


pathb = r'PATH_TO_BIG_ROUND_FEATURES\Buffer_5km.shp'
paths = r'PATH_TO_SMALL_ROUND_FEATURES\Buffer_2km.shp'
pathout = r'PATH_TO_OUTPUT_FOLDER\clip.shp'

big = fiona.open(pathb)
small = fiona.open(paths)

polbig = big.next()
polsmall = small.next()

a = shape(polbig['geometry'])
b = shape(polsmall['geometry'])
c = a.difference(b)

source_schema = { 'geometry': 'Polygon', 'properties': {'id': 'int'}}
source_crs = {'no_defs': True, 'ellps': 'WGS84', 'datum': 'WGS84', 'proj': 'longlat'}
source_driver = "ESRI Shapefile"

with fiona.open(pathout, 'w', 'ESRI Shapefile', schema=source_schema, crs=source_crs) as out:
    out.write({
        'geometry': mapping(c),
        'properties': {'id': 123},
    })

And this produces the following output:

Which is exactly what I need!

Answer 1

Since you are using python you could also have a look at Shapely. It is a library for manipulation and analysis of geometric objects in the Cartesian plane.

For your problem you would proceed like

• Load both geometries into Shapely objects
• on the obeject you want to substract the other from call object.difference(other)
• save the result into a new variable

For an exampel see http://toblerity.org/shapely/manual.html#object.difference

I personally would recommend you using a library which has a python interface, especially since when you use an IDE for Development you get stuff like autocomplete and type checking.

Answer 2

I would utilize the GDAL SQLite dialect. Here comes an example:

The two polygons in the image are as WKT:

POLYGON (( -113 39, -113 42, -110 42, -110 39, -113 39 ))

POLYGON (( -112 40, -112 41, -111 41, -111 40, -112 40 ))

You can make a hole into polygon 1 with polygon 2 by using Spatialite function ST_Difference https://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html

Demo with ogrinfo:

ogrinfo -ro -dialect sqlite -sql "select ST_Difference(ST_GeomFromText('POLYGON (( -113 39, -113 42, -110 42, -110 39, -113 39 ))'),ST_GeomFromText('POLYGON (( -112 40, -112 41, -111 41, -111 40, -112 40 ))'))" multipolygon.json

INFO: Open of `multipolygon.json'
      using driver `GeoJSON' successful.

Layer name: SELECT
Geometry: Unknown (any)
Feature Count: 1
Extent: (-113.000000, 39.000000) - (-110.000000, 42.000000)
Layer SRS WKT:
(unknown)
Geometry Column = ST_Difference(ST_GeomFromText('POLYGON (( -113 39, -113 42, -1
10 42, -110 39, -113 39 ))'),ST_GeomFromText('POLYGON (( -112 40, -112 41, -111
41, -111 40, -112 40 ))'))
OGRFeature(SELECT):0
  POLYGON ((-113 39,-113 42,-110 42,-110 39,-113 39),(-112 40,-111 40,-111 41,-1
12 41,-112 40))

The result is

POLYGON ((-113 39,-113 42,-110 42,-110 39,-113 39),(-112 40,-111 40,-111 41,-112 41,-112 40))

And it looks like this

What you need to do is to build the SQL with Python and execute it.

Answer 3

As I mentioned in my comment you could split out the inner and larger polygons into two separate layers and use either Difference tool (QGIS) or Erase tool (ArcGIS) to run in a single process. If you want, these two tools could be configured to run in a python script too.