5

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!

  • 1
    Can you not do it through geoprocessing process? Example, split out your inner and outer polygons into two separate layers and use either the Difference (QGIS), or Erase (ArcGIS) tool? – artwork21 Sep 15 '15 at 16:30
  • I did not think about that, but the question is doing that for the 25.000 features, not only this one. When you say geoprocessing, you mean Python inside either QGIS or ArcGIS? Can you point me to any example? I am totally new on that! Thanks! – Irene Sep 15 '15 at 20:31
6

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.

  • 1
    I didn't know about this library (together with Fiona) and how easy is to handle shapefiles with this! Python GDAL/OGR has a tricky and definitely not very pythonic syntax. I will edit my response to include the new code with this neat library. Thanks! – Irene Sep 16 '15 at 13:57
3

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

enter image description here

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

enter image description here

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

1

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.

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.