python - 溶解重叠多边形(使用 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 来处理 shapefile 通常是更好的方法。我在下面发布了我的解决方案。

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"