Python简化Polygon Shapefile：解决邻接多边形平滑对齐问题

问题背景 #

我有一个栅格矢量化得到的Polygon Shapefile，包含数千个要素，仅含一个代表多边形等级（1-5）的属性列，文件体积过大。尝试用Shapely的simplify工具减小文件大小，目标是达到ArcGIS Pro简化多边形工具的效果（可将文件大小减小80%）。目前Shapely处理后的文件大小相近，但邻接多边形无法平滑对齐，推测需要平滑或插值函数，但未找到合适方案。

现有代码 #

from shapely.geometry import shape, mapping
from shapely.ops import transform
from shapely.validation import make_valid
import fiona
import pyproj

# Define the projections
project_to_utm = pyproj.Transformer.from_crs("EPSG:4326", "EPSG:32633", always_xy=True).transform
project_to_wgs84 = pyproj.Transformer.from_crs("EPSG:32633", "EPSG:4326", always_xy=True).transform

# Open your shapefile
with fiona.open(r"shapefile.shp", 'r') as source:
    schema = source.schema
    crs = source.crs
    features = list(source)

    # Reproject to UTM, simplify, then reproject back to WGS84
    simplified_features = []
    for feature in features:
        # Reproject to UTM
        geom = shape(feature['geometry'])
        geom_utm = transform(project_to_utm, geom)

        # Simplify in UTM
        simplified_geom_utm = geom_utm.simplify(tolerance=0.5, preserve_topology=True)

        # Fix any invalid geometries (self-intersections)
        if not simplified_geom_utm.is_valid:
            simplified_geom_utm = make_valid(simplified_geom_utm)

        # Reproject back to WGS84
        simplified_geom_wgs84 = transform(project_to_wgs84, simplified_geom_utm)

        # Append the simplified and validated geometry
        simplified_features.append({
            'geometry': mapping(simplified_geom_wgs84),
            'properties': feature['properties']
        })

# Save the simplified polygons to a new shapefile
with fiona.open(r"shapfile_simplified.shp", 'w', driver='ESRI Shapefile', schema=schema, crs=crs) as output:
    for feature in simplified_features:
        output.write(feature)

效果对比 #

• Python输出矢量结果（邻接多边形错位）
• 栅格矢量化后的原始矢量
• ArcGIS Pro简化多边形工具输出结果（邻接多边形平滑对齐）

解决方案 #

Shapely单独简化每个要素会破坏拓扑一致性，因为邻接边的简化结果无法同步。要解决这个问题，需采用基于拓扑的批量简化方案，确保共享边保持一致：

方案1：使用TopoJSON进行拓扑简化 #

TopoJSON通过共享边存储拓扑关系，简化时同步处理共享边界，完美解决错位问题：

1. 安装依赖：pip install topojson geopandas
2. 代码实现：

import geopandas as gpd
from topojson import TopoJSON

# 读取原始Shapefile并转换为UTM投影（保证简化单位统一）
gdf = gpd.read_file("shapefile.shp")
gdf_utm = gdf.to_crs("EPSG:32633")

# 生成TopoJSON并保留拓扑简化
topo = TopoJSON(gdf_utm, topology=True)
simplified_topo = topo.simplify(tolerance=0.5, preserve_topology=True)

# 转换回GeoDataFrame并重新投影到WGS84
simplified_gdf = simplified_topo.to_gdf()
simplified_gdf = simplified_gdf.to_crs("EPSG:4326")

# 保存结果
simplified_gdf.to_file("shapefile_simplified_topojson.shp")

方案2：用缓冲操作修复拓扑错位 #

若不想引入TopoJSON依赖，可先简化再通过缓冲闭合邻接缝隙：

import geopandas as gpd

# 读取数据并投影到UTM
gdf = gpd.read_file("shapefile.shp").to_crs("EPSG:32633")

# 批量简化要素
gdf['geometry'] = gdf['geometry'].simplify(tolerance=0.5, preserve_topology=True)

# 缓冲修复拓扑：先微小扩张再收缩，闭合邻接缝隙
gdf['geometry'] = gdf['geometry'].buffer(0.1).buffer(-0.1)

# 修复无效几何并过滤异常要素
gdf = gdf[gdf.is_valid]
gdf['geometry'] = gdf['geometry'].make_valid()

# 投影回WGS84并保存
gdf.to_crs("EPSG:4326").to_file("shapefile_simplified_fixed.shp")

方案3：调用ArcGIS Pro原生工具 #

如果已安装ArcGIS Pro，可直接调用其简化工具，原生支持拓扑一致性：

import arcpy

# 设置工作空间
arcpy.env.workspace = "your_workspace_path"
# 调用简化工具，开启拓扑保持参数
arcpy.cartography.SimplifyPolygon(
    in_features="shapefile.shp",
    out_feature_class="shapefile_simplified_arcpy.shp",
    algorithm="POINT_REMOVE",
    tolerance=0.5,
    minimum_area="0 SquareMeters",
    simplify_connected_features=True,  # 关键：保持连接要素的拓扑一致性
    collapsed_point_option="NO_KEEP"
)

关键说明 #

• 拓扑一致性的核心是共享边同步处理，单独简化每个要素必然导致错位；
• 优先推荐TopoJSON方案，其简化效果最接近ArcGIS Pro；
• 缓冲修复方案适合轻量场景，但需根据数据精度调整缓冲值，避免过度改变多边形面积。

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