Растворите перекрывающиеся многоугольники (с помощью GDAL/OGR), оставив несвязанные результаты отличными

Есть ли способ слить (слить) перекрывающиеся полигоны, используя любой 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 для работы с шейп-файлами. Я разместил свое решение ниже.

Источник

2 ответа

Решение

Итак, после многих неудачных попыток я отбросил gdal/ogr и продолжил стройную и фиону. Это именно то, что мне нужно. Фильтрация была необходима, потому что мой набор данных содержит самопересекающиеся многоугольники, которые необходимо отфильтровать перед вызовом. 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}})
Источник

У меня была похожая проблема, и я решил ее следующим образом, взяв все предыдущие ответы из книги "Объединение нескольких геометрий в 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 как тип геометрии, поэтому я реализовал опцию для повторного создания отдельных многоугольников. Вы также можете использовать ogr2ogr из командной строки с опцией -explodecollections:

ogr2ogr -f "ESRI Shapefile" -explodecollections dissolved.shp input.shp -dialect sqlite -sql "select ST_union(Geometry) from input"
Источник
Другие вопросы по тегам