Monday, 3 October 2016

python - Union multiple shapefiles with overlapping shapes


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.



Answer



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)

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...