Search code examples
pythongdalshapelyogrfiona

Dissolve Overlapping Polygons (with GDAL/OGR) while keeping non-connected results distinct


Is there any way to dissolve (merge) overlapping polygons, using any GDAL/OGR API or command line tool, while keeping the resulting non overlapping areas distinct? I have searched a lot and I cannot find anything resembling what need. However, I think it is very unlikely that this problem has not been solved yet.

Here's a more detailed description of what I need:

  • My input consists of a single shape file (ESRI Shapefile) with a single layer.
  • This layer contains polygons which are not distinguishable by attributes. (All have the same attributes).
  • Many of them are overlapping and I would like to get the union of those who are overlapping.
  • Areas that are not connected should result in separate polygons.

It is the last point which is causing troubles. I basically get what I need except for the last point. If I run the typical solution for dissolving a shape file

$ ogr2ogr -f "ESRI Shapefile" dissolved.shp input.shp -dialect sqlite -sql "select ST_union(Geometry) from input"

I end up with a single polygon which includes everything, even if the areas are not connected.

Update: I solved the problem by ditching GDAL entirely. As many sources point out it is generally a better approach to use fiona and shapely to work with shapefiles. I have posted my solution below.


Solution

  • So, after many unsuccessful attempts I have ditched gdal/ogr and went on with shapely and fiona. This does exactly what I need. The filtering was necessary becuase my dataset contains self-intersecting polygons which need to be filtered out before calling cascaded_union.

    import fiona                                                                                                       
    from shapely.ops import cascaded_union                                                                             
    from shapely.geometry import shape, mapping  
    
    with fiona.open(src, 'r') as ds_in:                                                                                                                                                                                                                   
        crs = ds_in.crs 
        drv = ds_in.driver 
    
        filtered = filter(lambda x: shape(x["geometry"]).is_valid, list(ds_in))                                                                                   
    
        geoms = [shape(x["geometry"]) for x in filtered]                                                   
        dissolved = cascaded_union(geoms)                                    
    
    schema = {                                                                                                     
        "geometry": "Polygon",                                                                                     
        "properties": {"id": "int"}                                                                                
    }  
    
    with fiona.open(dst, 'w', driver=drv, schema=schema, crs=crs) as ds_dst:                                       
        for i,g in enumerate(dissolved):                                                                           
            ds_dst.write({"geometry": mapping(g), "properties": {"id": i}})