[GIS] Erasing features using Python GDAL/OGR

clipgdalogrpython

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!

enter image description here

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:

enter image description here

Which is exactly what I need!

Best Answer

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.

Related Question