Erasing features using Python GDAL/OGR?

Erasing features using Python GDAL/OGR?

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_FEATURESBuffer_5km.shp' paths = r'PATH_TO_SMALL_ROUND_FEATURESBuffer_2km.shp' pathout = r'PATH_TO_OUTPUT_FOLDERclip.shp' big = small = polbig = polsmall = 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, '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!

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 callobject.difference(other)
  • save the result into a new variable

For an exampel see

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.

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

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.

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.