使用Geopandas在Python中分割重叠的Multipolygon多边形

需求说明 #

将所有位于主多边形内部的子多边形合并为主多边形的内洞，外部多边形保留独立状态，最终得到带内洞的主多边形+独立外部多边形的结果。

示例数据 #

data_01 = {
    "type": "FeatureCollection",
    "features": [
        {
            "type": "Feature",
            "geometry": {
                "type": "Polygon",
                "coordinates": [
                    [[2, 2], [2, 22], [22, 22], [22, 2], [2, 2]]
                ]
            },
            "properties": {"z": 1412.5, "la": "ba"}
        },
        {
            "type": "Feature",
            "geometry": {
                "type": "Polygon",
                "coordinates": [
                    [[12, 16], [7, 10], [17, 10], [12, 16]]
                ]
            },
            "properties": {"z": 1412.5, "la": "ba"}
        },
        {
            "type": "Feature",
            "geometry": {
                "type": "Polygon",
                "coordinates": [
                    [[27, 15], [24, 12], [29, 12], [27, 15]]
                ]
            },
            "properties": {"z": 1412.5, "la": "ba"}
        }
    ]
}

现有问题 #

原代码仅能处理单个内部多边形，无法适配任意数量的内部/外部多边形场景。

通用解决方案代码 #

import geopandas as gpd
import json

# 加载数据（示例直接使用给定数据，实际可替换为gpd.read_file("data_pol.json")）
gdf = gpd.GeoDataFrame.from_features(data_01["features"])

# 自动识别主多边形：取面积最大的多边形
main_idx = gdf.geometry.area.idxmax()
main_geometry = gdf.loc[main_idx, "geometry"]
main_properties = gdf.loc[main_idx, "properties"].to_dict()

inner_polygons = []
outer_features = []

# 遍历所有多边形，区分内部/外部
for idx, row in gdf.iterrows():
    if idx == main_idx:
        continue
    if row.geometry.within(main_geometry):
        inner_polygons.append(row.geometry)
    else:
        # 构造外部Feature
        outer_features.append({
            "type": "Feature",
            "geometry": row.geometry.__geo_interface__,
            "properties": row.properties.to_dict()
        })

# 批量减去内部多边形，生成带内洞的主多边形
for poly in inner_polygons:
    main_geometry = main_geometry.difference(poly)

# 构造最终结果
final_features = outer_features.copy()
final_features.append({
    "type": "Feature",
    "geometry": main_geometry.__geo_interface__,
    "properties": main_properties
})

final_data = {
    "type": "FeatureCollection",
    "features": final_features
}

# 验证结果（可选）
print(json.dumps(final_data, indent=2))

代码说明 #

• 自动识别主多边形：通过面积最大值判断，无需手动指定索引，适配不同数据场景
• 批量处理内部多边形：循环执行difference操作，支持任意数量的内洞生成
• 保留外部多边形：直接收集外部多边形并加入最终结果，不做修改
• 输出格式合规：最终输出完全符合要求的GeoJSON FeatureCollection格式

内容的提问来源于stack exchange，提问作者user20567574