使用GDAL/OGR溶解重叠多边形，同时保持非连接结果不同。

Question

是否有任何方法可以使用GDAL/OGR API或命令行工具溶解（合并）重叠的多边形，同时保持结果中非重叠区域的明显性？我已经搜索了很多，但我找不到类似于我需要的东西。不过，我认为这个问题很不可能还没有被解决。

以下是我需要的更详细的描述：

- 我的输入包括一个单一的形状文件（ESRI Shapefile），其中包含一个图层。 - 该图层包含的多边形在属性上无法区分（所有多边形具有相同的属性）。 - 其中许多是重叠的，我想获取那些重叠部分的联合。 - 没有连接的区域应该生成单独的多边形。

正是最后一点引起了麻烦。基本上，我获得了我需要的，除了最后一点。如果我运行溶解形状文件的典型解决方案

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

我最终得到了一个单一的多边形，其中包括所有区域，即使这些区域并不相连。

更新： 我通过完全放弃GDAL来解决了这个问题。正如许多来源指出的那样，使用fiona和shapely处理shapefiles通常是更好的方法。我在下面发布了我的解决方案。

- StarSheriff

Answer 1

我曾经遇到类似的问题，通过参考 Taking Union of several geometries in GEOS Python? 的所有先前答案来解决它：

import os
from osgeo import ogr

def createDS(ds_name, ds_format, geom_type, srs, overwrite=False):
    drv = ogr.GetDriverByName(ds_format)
    if os.path.exists(ds_name) and overwrite is True:
        deleteDS(ds_name)
    ds = drv.CreateDataSource(ds_name)
    lyr_name = os.path.splitext(os.path.basename(ds_name))[0]
    lyr = ds.CreateLayer(lyr_name, srs, geom_type)
    return ds, lyr

def dissolve(input, output, multipoly=False, overwrite=False):
    ds = ogr.Open(input)
    lyr = ds.GetLayer()
    out_ds, out_lyr = createDS(output, ds.GetDriver().GetName(), lyr.GetGeomType(), lyr.GetSpatialRef(), overwrite)
    defn = out_lyr.GetLayerDefn()
    multi = ogr.Geometry(ogr.wkbMultiPolygon)
    for feat in lyr:
        if feat.geometry():
            feat.geometry().CloseRings() # this copies the first point to the end
            wkt = feat.geometry().ExportToWkt()
            multi.AddGeometryDirectly(ogr.CreateGeometryFromWkt(wkt))
    union = multi.UnionCascaded()
    if multipoly is False:
        for geom in union:
            poly = ogr.CreateGeometryFromWkb(geom.ExportToWkb())
            feat = ogr.Feature(defn)
            feat.SetGeometry(poly)
            out_lyr.CreateFeature(feat)
    else:
        out_feat = ogr.Feature(defn)
        out_feat.SetGeometry(union)
        out_lyr.CreateFeature(out_feat)
        out_ds.Destroy()
    ds.Destroy()
    return True

UnionCascaded 需要 MultiPolygon 作为几何类型，这就是为什么我实现了重新创建单个多边形的选项。您也可以在命令行中使用选项 -explodecollections 的 ogr2ogr：

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

- s6hebern

Answer 2

经过多次尝试，我放弃了gdal/ogr，转而使用shapely和fiona。这正好满足我的需求。进行过滤是必要的，因为我的数据集包含自相交的多边形，在调用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}}) 

- StarSheriff

---

你好。这个脚本需要什么输入？我如何将其用于给定的shapefile？谢谢。 - Loonuh