溶解重叠多边形(GDAL/OGR),同时保持非连接结果不同
Dissolve Overlapping Polygons (with GDAL/OGR) while keeping non-connected results distinct
是否有任何方法可以使用任何 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 通常是更好的方法。我已经在下面发布了我的解决方案。
因此,在多次尝试失败后,我放弃了 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}})
我有一个类似的问题并像这样解决了它,考虑了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
作为几何类型,这就是为什么我实现了重新创建单个多边形的选项。您还可以在命令行中使用 ogr2ogr
选项 -explodecollections
:
ogr2ogr -f "ESRI Shapefile" -explodecollections dissolved.shp input.shp -dialect sqlite -sql "select ST_union(Geometry) from input"
是否有任何方法可以使用任何 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 通常是更好的方法。我已经在下面发布了我的解决方案。
因此,在多次尝试失败后,我放弃了 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}})
我有一个类似的问题并像这样解决了它,考虑了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
作为几何类型,这就是为什么我实现了重新创建单个多边形的选项。您还可以在命令行中使用 ogr2ogr
选项 -explodecollections
:
ogr2ogr -f "ESRI Shapefile" -explodecollections dissolved.shp input.shp -dialect sqlite -sql "select ST_union(Geometry) from input"