当前位置：首页 > news >正文

使用 Python 自动化检查矢量面数据的拓扑错误（含导出/删除选项）

news 2025/9/1 6:30:39

在地理信息系统（GIS）数据处理中，矢量面数据的拓扑关系至关重要。如果拓扑错误没有及时发现和修复，可能会影响后续空间分析、建模与可视化效果。本文介绍如何使用 GeoPandas + Shapely 进行拓扑检查，并提供扩展功能（打印、导出、删除问题要素）。

一、拓扑错误常见类型

在多边形矢量数据中，常见的拓扑问题包括：

非法几何（Invalid Geometry）
多边形自相交或几何结构错误。

带洞要素（Polygons with Holes）
多边形内部包含空洞（通常需要视业务是否允许）。

极小面（Small Polygons）
面积过小的碎片，通常是误差或无意义的边角。

完全重复几何（Duplicate Geometries）
两个要素几何完全相同。

同层重叠（Self-Overlap / Overlap within Layer）
两个多边形在同一层里相互叠压。

二、Python 拓扑检查脚本

# -*- coding: utf-8 -*-
"""
矢量面数据拓扑检查与清洗
记得先把数据加上投影奥
- 检查项：非法几何、带洞、极小面、重复几何、同层重叠
- 可选操作：打印 / 导出 / 删除
"""import geopandas as gpd
from shapely.validation import explain_validity
from shapely.geometry import Polygon, MultiPolygon, box
import shapely
import os# ========= 基本参数 =========
SRC = r""
MIN_AREA_M2 = 10.0            # 极小面阈值（平方米）
FIX_INVALID_BEFORE = True     # 检查前修复非法几何
EXPORT_PROBLEMS = True        # 是否导出问题要素
DELETE_PROBLEMS = False       # 是否直接删除问题要素（谨慎使用）OUT_DIR = os.path.dirname(SRC)# ========= 读取数据 =========
gdf = gpd.read_file(SRC)
print(f"读取完成：{len(gdf)} 要素；原CRS = {gdf.crs}")if gdf.crs is None:gdf = gdf.set_crs(epsg=4326, allow_override=True)print("未检测到CRS，已强制设置 EPSG:4326")# ========= 构造米制副本 =========
try:utm_crs = gdf.estimate_utm_crs()gdf_m = gdf.to_crs(utm_crs)print(f"用于量测的米制CRS：{utm_crs}")
except Exception as e:print("UTM估算失败，仍使用原CRS：", e)gdf_m = gdf.copy()# ========= 几何修复函数 =========
def fix_geom(geom):if geom is None or geom.is_empty:return geomif not geom.is_valid:try:return shapely.make_valid(geom)except Exception:return geom.buffer(0)return geomif FIX_INVALID_BEFORE:gdf["geometry"] = gdf["geometry"].apply(fix_geom)gdf_m["geometry"] = gdf_m["geometry"].apply(fix_geom)# ========= 检查1：非法几何 =========
invalid_idx = [i for i, geom in gdf.geometry.items() if not geom.is_valid]# ========= 检查2：带洞要素 =========
def has_holes(geom):if isinstance(geom, Polygon):return len(geom.interiors) > 0if isinstance(geom, MultiPolygon):return any(len(p.interiors) > 0 for p in geom.geoms)return Falseholes_idx = [i for i, geom in gdf.geometry.items() if has_holes(geom)]# ========= 检查3：极小面 =========
small_idx = [i for i, geom in gdf_m.geometry.items() if geom.area < MIN_AREA_M2]# ========= 检查4：完全重复几何 =========
dup_idx = list(gdf.index[gdf.geometry.duplicated(keep=False)])# ========= 检查5：同层重叠 =========
print("\n[同层重叠] 正在扫描 ...")
overlap_pairs, overlap_ids = [], set()
sindex = gdf_m.sindexfor i, geom in enumerate(gdf_m.geometry):if geom.is_empty:continuecandidate_idx = sindex.query(box(*geom.bounds))for j in candidate_idx:if j <= i:continueg2 = gdf_m.geometry.iloc[j]inter = geom.intersection(g2)if not inter.is_empty and inter.area > 0:overlap_pairs.append((i, j, float(inter.area)))overlap_ids.update([i, j])# ========= 汇总 =========
print("\n—— 拓扑检查汇总 ——")
print(f"非法几何：{len(invalid_idx)}")
print(f"带洞要素：{len(holes_idx)}")
print(f"极小面(<{MIN_AREA_M2} m²)：{len(small_idx)}")
print(f"完全重复几何：{len(dup_idx)}")
print(f"同层重叠：{len(overlap_pairs)} 对，涉及 {len(overlap_ids)} 个要素")# ========= 导出 / 删除选项 =========
problem_idx = sorted(set(invalid_idx + holes_idx + small_idx + dup_idx + list(overlap_ids)))if problem_idx:print(f"\n⚠️ 检测到 {len(problem_idx)} 个问题要素：{problem_idx[:20]}{' ...' if len(problem_idx)>20 else ''}")if EXPORT_PROBLEMS:out_path = os.path.join(OUT_DIR, "problems.shp")gdf.loc[problem_idx].to_file(out_path, encoding="utf-8")print(f"已导出问题要素到: {out_path}")if DELETE_PROBLEMS:gdf_clean = gdf.drop(problem_idx)out_clean = os.path.join(OUT_DIR, "cleaned.shp")gdf_clean.to_file(out_clean, encoding="utf-8")print(f"已删除问题要素，结果保存到: {out_clean}")
else:print("\n✅ 数据没有发现拓扑问题！")