7

I have several shapefiles that overlap and I would like to union their geometries to get rid of those overlaps. I know how to do this with OGR and Python for single geometries.

from osgeo import ogr

wkt1 = "POLYGON((-755226.75354172603692859 -1009167.14398036629427224,-755262.54693246085662395 -1016290.02873659669421613,-742770.65356600657105446 -1016111.06178292259573936,-743629.69494364247657359 -1009453.49110624496825039,-755226.75354172603692859 -1009167.14398036629427224))"

wkt2 = "POLYGON((-746457.3728116936981678 -1007413.26783435989636928,-747065.86045418574940413 -1012961.24339825788047165,-739978.7690886901691556 -1012889.65661678824108094,-740193.52943309908732772 -1007771.20174170809332281,-746457.3728116936981678 -1007413.26783435989636928))"

poly1 = ogr.CreateGeometryFromWkt(wkt1)
poly2 = ogr.CreateGeometryFromWkt(wkt2)

union = poly1.Union(poly2)

print poly1
print poly2
print union.ExportToJson()

How do I do that for the whole shapefile? What would be the easiest way? I'm looking for console-based solutions. Thanks for hints.

6

You can do that by first creating an empty geometry and "unioning" it with the polygons of the shapefile:

from osgeo import ogr
test = ogr.Open('polygons.shp')
layer = test.GetLayer()
print layer.GetGeomType()
3 # -> polygons
# empty geometry
union_poly = ogr.Geometry(ogr.wkbPolygon)
# make the union of polygons
for feature in layer:
      geom =feature.GetGeometryRef()
      union_poly = union_poly.Union(geom)

Example:

  • the shapefile:

enter image description here

  • result of the union:

enter image description here

Or you can also use Fiona, an other wrapper of the OGR library, and shapely which has unionbut also cascaded_union and unary_union, faster:

import fiona
from shapely.geometry import shape
polygons =[shape(feature['geometry') for feature in fiona.open("polygons.shp")]
from shapely.ops import cascaded_union, unary_union
union_poly = cascaded_union(polygons) # or unary_union(polygons)
  • Fiona example was very useful. FYI: minor typo in line: polygons =[shape(feature... – python1981 Dec 7 '15 at 23:18
3

With the OGR's function UnionCascaded()

from osgeo import ogr
shp= "yourfile"

shapef = ogr.Open(shpor)
lyr = shapef.GetLayer()
unionc = ogr.Geometry(ogr.wkbMultiPolygon)
for feat in xrange(lyr.GetFeatureCount()):
    fit= lyr.GetFeature(feat)
    geom= fit.GetGeometryRef()
    unionc.AddGeometry(geom)
union= unionc.UnionCascaded()

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.