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:
- result of the
union
:
Or you can also use Fiona, an other wrapper of the OGR library, and shapely which has union
but 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