上次介绍的河道管理线处理程序是完全独立于GIS软件的python程序，在实际使用中可能会面临操作系统和环境不同而导致的启动问题。因此，这回介绍一下如何在QGIS中利用python控制台实现批量处理河道管理线的自动闭合转面问题。
本文将分享基于QGIS3.44版Python控制台的轻量级解决方案，通过480行代码，实现从断点检测、中线过滤到自动封闭转面生成的全流程批处理，实测50余条河道线批量处理仅用时6.5秒完成。
一、问题聚焦与工具目标
审计人员拿到的河道管理范围线常出现三类典型问题：
1.边界不封闭：左右岸线之间存在未连接的缺口，无法直接构成多边形，导致面积无法计算。
2.冗余中线：部分数据将河道中心线混在边界线图层中，若不剔除，会严重干扰后续多边形化。
3.属性缺失：原始线段有河道名称、编码等属性，但修复后的面域往往丢失这些信息，需要人工回补。
我们的工具需在保证处理效率和结果准确的前提下，自动化解决上述问题，具体要求：
二、技术选型与设计理念
工具完全基于QGIS原生Python API开发，不依赖任何第三方库。QGIS具备强大的几何处理引擎，其内置的polygonize算法和空间索引，可以让我们在几行代码内完成复杂的拓扑重建。选择Python控制台作为运行环境，则大幅降低了使用门槛——无需安装插件，直接复制脚本即可运行。
设计上遵循三个原则：
轻量：代码紧凑，逻辑清晰，便于审计人员理解和二次修改。
务实：不追求大而全，聚焦解决“封闭性修复”这一核心痛点，输出成果直接可用于面积统计和进一步分析。
智能容错：提供多种备用算法（如内置多边形化失败时自动调用处理工具箱），并对中线识别采用全字段模糊匹配，适应不同字段命名。
三、核心算法分步解析
1. 河道中线智能过滤
中线通常会贯穿整个河道，且多数数据会有一个字段（如“类型”“类别”）标记其为中心线。脚本会遍历图层的所有字段，对每个字段的值进行正则匹配（关键词：中线|中心线|centerline），一旦命中即标记为中线并剔除。若属性中缺乏标识，则跳过该步骤，避免误删真实边界。分离出的中线另存为独立图层，供人工核验。
2. 断点检测与最近邻连接
将图层中的所有线（含多部件）拆解为单线，提取每条线的首尾端点，统计每个端点在全局出现的次数。出现次数为1的端点即为悬空断点。这些断点往往是数据采集或转换时产生的缺口。
连接策略采用贪心最近邻配对：计算所有断点两两之间的距离并升序排序，依次取出距离最小且均未被连接的断点对，用直线连接。该算法计算复杂度为O（n² log n），但实际断点数量通常不超过几百个，性能完全满足要求。极少数奇数个断点的情况会给出警告，不影响其他配对。
3. 双引擎多边形化
将原始边界线（已剔除中线）与自动生成的连接线合并，利用QGIS的polygonize函数生成面。但受浮点误差或特殊几何形态影响，有时内置方法会失败。此时工具会自动启用备选方案：将线图层保存为临时Shapefile，调用QGIS处理工具箱的native:polygonize算法再次尝试。双层保障使封闭面生成成功率大幅提升。
4. 属性继承与面积计算
生成的每一个多边形，通过空间索引查询其内部点（pointOnSurface）最近邻的原始线要素，将该线的所有属性（如河道名称、编码）赋给面，确保审计证据链完整。同时利用QgsDistanceArea自动计算面积，对于地理坐标系还考虑了椭球面积，直接输出结果。
四、批处理交互界面
为了让非编程人员也能轻松使用，脚本封装了一个PyQt5对话框（见图1，文中略），提供：
目录浏览：递归搜索子文件夹中所有支持的矢量文件。
四个开关选项：可分别控制是否生成中线图层、断点标注图层、修复后线图层，而“河道转面层”默认必选且不可取消。
实时日志窗口：显示每个文件的处理进度、中线数量、断点数量、成功/失败状态。
统计总结：所有文件处理完毕后汇总成功数、失败数、总断点数、总面数及耗时。
对话框采用非模态设计，处理过程中可自由切换QGIS主界面查看结果。
五、实测表现与适用场景
在QGIS 3.44.12版本下测试，我们使用一个包含53条河道的管理范围线数据进行测试，所有文件均为Shapefile格式，每条河包含数量不等的节点。工具自动识别并分离出12条中线，检测到108个断点，生成53个封闭面。总耗时6.5秒，平均每文件仅0.12秒。生成的面域属性完整，名称、编码与原始数据一致，面积计算准确。人工检查未发现明显拓扑错误。
该工具尤其适用于一次性处理大量标准化的河道界线数据，例如每年自然资源审计中需要对多个流域范围线进行快速验证和面积汇总的场景。对于少数自动封闭出现异常结果的情况（一边岸线延伸特别长、断点配对错误等复杂情况），工具会保留修复后的线图层，用户可手工微调后利用QGIS自带工具完成收尾。
六、总结
这套脚本源于实际审计项目中遇到的痛点，经过多个版本的迭代优化，最终形成了一套可复用的高效工具。它没有依赖任何外部包，全部代码在一个Python控制台中即可运行，充分体现了“用最简单方法解决实际问题”的理念。诚然，算法层面仍有可优化之处（例如断点连接可采用最小生成树避免潜在交叉），但目前的贪心策略在绝大多数河网较为简单的场景下已足够可靠。
欢迎同行试用、反馈，共同打磨更完善的审计数据处理工具集。
（代码全文附于文末，可直接复制到QGIS Python控制台试运行。）
完整代码如下：
from qgis.core import (
QgsProject, QgsVectorLayer, QgsField, QgsFeature,
QgsGeometry, QgsWkbTypes, QgsPointXY, QgsFields, Qgis,
QgsVectorFileWriter, QgsSpatialIndex, QgsDistanceArea
)
from qgis.utils import iface
from PyQt5.QtWidgets import (
QDialog, QVBoxLayout, QHBoxLayout, QLineEdit, QPushButton,
QCheckBox, QTextEdit, QLabel, QFileDialog
)
from PyQt5.QtCore import QVariant, Qt
from PyQt5.QtGui import QColor, QFont
from collections import defaultdict
import math, processing, tempfile, os, re, time, glob
# ------------------------------------------------------------
# 辅助函数：自动查找名称字段
# ------------------------------------------------------------
def find_name_field(layer):
name_pattern = re.compile(r'名称|name|河名|河道名|mc', re.IGNORECASE)
fields = layer.fields()
for i, field in enumerate(fields):
if name_pattern.search(field.name()):
return i, field.name()
for i, field in enumerate(fields):
if field.type() == QVariant.String:
return i, field.name()
return None
# ------------------------------------------------------------
# 中线识别
# ------------------------------------------------------------
def filter_centerlines(src_layer, log_func=print):
crs = src_layer.crs()
fields = src_layer.fields()
temp_layer = QgsVectorLayer(f"LineString?crs={crs.authid()}", "temp_nomidline", "memory")
pr = temp_layer.dataProvider()
pr.addAttributes(fields)
temp_layer.updateFields()
for f in src_layer.getFeatures():
pr.addFeature(QgsFeature(f))
temp_layer.updateExtents()
pattern = re.compile(r'中线|中心线|centerline|center\s*line', re.IGNORECASE)
midline_ids = []
for f in temp_layer.getFeatures():
for field in fields:
val = str(f.attribute(field.name()))
if pattern.search(val):
midline_ids.append(f.id())
break
if not midline_ids:
log_func(" 未发现中线标识。")
return temp_layer, None
log_func(f" 发现 {len(midline_ids)} 条可能的中线，正在分离...")
midline_layer = QgsVectorLayer(f"LineString?crs={crs.authid()}", "temp_midline", "memory")
pr_mid = midline_layer.dataProvider()
pr_mid.addAttributes(fields)
midline_layer.updateFields()
midline_feats = [QgsFeature(f) for f in temp_layer.getFeatures() if f.id() in midline_ids]
pr_mid.addFeatures(midline_feats)
midline_layer.updateExtents()
temp_layer.startEditing()
for fid in sorted(midline_ids, reverse=True):
temp_layer.deleteFeature(fid)
temp_layer.commitChanges()
return temp_layer, midline_layer
# ------------------------------------------------------------
# 单文件处理核心（返回统计字典或None）
# ------------------------------------------------------------
def process_single_file(file_path, options, log_func=print):
"""
options: dict with keys: 'gen_midline', 'gen_breakpoints', 'gen_repaired', 'gen_polygon'
"""
file_name = os.path.splitext(os.path.basename(file_path))[0]
log_func(f"\n开始处理文件：{file_path}")
# 加载图层
src_layer = QgsVectorLayer(file_path, file_name, "ogr")
if not src_layer.isValid() or src_layer.geometryType() != QgsWkbTypes.LineGeometry:
log_func(f" 错误：文件无效或不是线图层，跳过。")
return None
start_time = time.time()
crs = src_layer.crs()
src_fields = src_layer.fields()
# 中线过滤
filtered, midline_layer = filter_centerlines(src_layer, log_func)
midline_count = 0
if midline_layer:
midline_count = midline_layer.featureCount()
if options.get('gen_midline', False):
midline_layer.setName(f"{file_name}_中线")
QgsProject.instance().addMapLayer(midline_layer)
log_func(f" 已添加中线图层：{midline_layer.name()}，共 {midline_count} 条。")
else:
log_func(f" 中线未生成（选项关闭）。")
# 构建空间索引（使用原始图层用于属性传递）
src_index = QgsSpatialIndex()
src_features = list(src_layer.getFeatures())
for f in src_features:
src_index.addFeature(f)
name_field_info = find_name_field(src_layer)
name_idx = name_field_info[0] if name_field_info else None
# ---------- 提取所有线 ----------
all_lines = []
endpoint_counter = defaultdict(int)
for feat in filtered.getFeatures():
geom = feat.geometry()
if not geom:
continue
attrs = feat.attributes()
wkb = geom.wkbType()
if wkb == QgsWkbTypes.LineString:
coords = geom.asPolyline()
all_lines.append((QgsGeometry.fromPolylineXY(coords), attrs, True))
if len(coords) >= 2:
s, e = coords[0], coords[-1]
endpoint_counter[(s.x(), s.y())] += 1
endpoint_counter[(e.x(), e.y())] += 1
elif wkb == QgsWkbTypes.MultiLineString:
for part in geom.asMultiPolyline():
all_lines.append((QgsGeometry.fromPolylineXY(part), attrs, True))
if len(part) >= 2:
s, e = part[0], part[-1]
endpoint_counter[(s.x(), s.y())] += 1
endpoint_counter[(e.x(), e.y())] += 1
# ---------- 断点处理 ----------
break_points = [pt for pt, cnt in endpoint_counter.items() if cnt == 1]
break_count = len(break_points)
if break_count == 0:
log_func(" 未发现断点，线路已封闭。")
all_line_geoms_with_attrs = all_lines
else:
log_func(f" 检测到 {break_count} 个断点。")
# 是否生成断点图层
if options.get('gen_breakpoints', False):
pt_fields = QgsFields()
pt_fields.append(QgsField("id", QVariant.Int))
pt_fields.append(QgsField("河道名称", QVariant.String, len=254))
pt_fields.append(QgsField("经度", QVariant.Double, "double", 20, 8))
pt_fields.append(QgsField("纬度", QVariant.Double, "double", 20, 8))
pt_layer = QgsVectorLayer(f"Point?crs={crs.authid()}", f"{file_name}_断点", "memory")
pr_pt = pt_layer.dataProvider()
pr_pt.addAttributes(pt_fields)
pt_layer.updateFields()
pt_feats = []
for idx, (x, y) in enumerate(break_points):
f = QgsFeature()
f.setGeometry(QgsGeometry.fromPointXY(QgsPointXY(x, y)))
nearest_id = src_index.nearestNeighbor(QgsPointXY(x, y), 1)
river_name = ""
if nearest_id:
nearest_feat = src_features[nearest_id[0]]
if name_idx is not None:
river_name = str(nearest_feat.attribute(name_idx))
f.setAttributes([idx + 1, river_name, x, y])
pt_feats.append(f)
pr_pt.addFeatures(pt_feats)
pt_layer.updateExtents()
QgsProject.instance().addMapLayer(pt_layer)
# 样式
renderer = pt_layer.renderer()
symbol = renderer.symbol()
if symbol:
symbol.setColor(QColor(255, 0, 0))
symbol.setSize(3)
pt_layer.triggerRepaint()
log_func(f" 已添加断点图层：{pt_layer.name()}。")
else:
log_func(" 断点图层未生成（选项关闭）。")
# 连接线
n = break_count
pts = [QgsPointXY(x, y) for (x, y) in break_points]
edges = []
for i in range(n):
for j in range(i + 1, n):
d = math.sqrt((pts[i].x() - pts[j].x()) ** 2 + (pts[i].y() - pts[j].y()) ** 2)
edges.append((d, i, j))
edges.sort(key=lambda x: x[0])
used = [False] * n
connect_lines = []
empty_attrs = [None] * len(src_fields)
for d, i, j in edges:
if not used[i] and not used[j]:
used[i] = used[j] = True
geom = QgsGeometry.fromPolylineXY([pts[i], pts[j]])
connect_lines.append((geom, empty_attrs, False))
if any(not u for u in used):
log_func(f" 警告：{sum(not u for u in used)} 个断点未能配对。")
all_line_geoms_with_attrs = all_lines + connect_lines
# ---------- 修复线图层 ----------
if options.get('gen_repaired', False):
repaired_fields = QgsFields()
repaired_fields.append(QgsField("要素来源", QVariant.String, len=20))
for field in src_fields:
repaired_fields.append(field)
repaired_layer = QgsVectorLayer(f"LineString?crs={crs.authid()}", f"{file_name}_修复", "memory")
pr_line = repaired_layer.dataProvider()
pr_line.addAttributes(repaired_fields)
repaired_layer.updateFields()
line_feats = []
for geom, attrs, is_original in all_line_geoms_with_attrs:
f = QgsFeature()
f.setGeometry(geom)
if is_original:
new_attrs = ["原始边界"] + attrs
else:
new_attrs = ["连接线"] + ([None] * len(src_fields))
f.setAttributes(new_attrs)
line_feats.append(f)
pr_line.addFeatures(line_feats)
repaired_layer.updateExtents()
QgsProject.instance().addMapLayer(repaired_layer)
renderer = repaired_layer.renderer()
symbol = renderer.symbol()
if symbol:
symbol.setColor(QColor(0, 0, 0))
symbol.setWidth(0.5)
repaired_layer.triggerRepaint()
log_func(f" 已添加修复线图层：{repaired_layer.name()}。")
else:
log_func(" 修复线图层未生成（选项关闭）。")
# ---------- 多边形化（必选） ----------
snapped_geoms = [g.snappedToGrid(0.0001, 0.0001) for (g, _, _) in all_line_geoms_with_attrs]
polygon_geom = QgsGeometry.polygonize(snapped_geoms)
polygons = []
if polygon_geom.isEmpty():
log_func(" 内置 polygonize 失败，尝试工具箱...")
try:
# 需要临时写入修复后的线图层，如果用户没有生成修复图层，则写一个临时图层
tmp_layer = QgsVectorLayer(f"LineString?crs={crs.authid()}", "temp", "memory")
tmp_pr = tmp_layer.dataProvider()
tmp_pr.addAttributes([QgsField("dummy", QVariant.Int)])
tmp_layer.updateFields()
feats = []
for geom, _, _ in all_line_geoms_with_attrs:
f = QgsFeature()
f.setGeometry(geom)
feats.append(f)
tmp_pr.addFeatures(feats)
tmp_layer.updateExtents()
tmp_dir = tempfile.mkdtemp()
tmp_path = os.path.join(tmp_dir, "temp_lines.shp")
QgsVectorFileWriter.writeAsVectorFormat(tmp_layer, tmp_path, "UTF-8", crs, "ESRI Shapefile")
result = processing.run("native:polygonize", {
'INPUT': tmp_path,
'KEEP_FIELDS': False,
'OUTPUT': 'memory:'
})
poly_temp = result['OUTPUT']
if poly_temp.featureCount() > 0:
polygons = [f.geometry() for f in poly_temp.getFeatures()]
except Exception as e:
log_func(f" 多边形化错误：{e}")
return None
else:
if polygon_geom.isMultipart():
polygons = [g for g in polygon_geom.asGeometryCollection() if g.type() == QgsWkbTypes.PolygonGeometry]
else:
if polygon_geom.type() == QgsWkbTypes.PolygonGeometry:
polygons = [polygon_geom]
if not polygons:
log_func(" 无法生成封闭多边形。")
return None
polygon_count = len(polygons)
# 生成面图层（必选）
area_calc = QgsDistanceArea()
area_calc.setSourceCrs(crs, QgsProject.instance().transformContext())
if crs.isGeographic():
area_calc.setEllipsoid(crs.ellipsoidAcronym())
poly_fields = QgsFields()
for field in src_fields:
poly_fields.append(QgsField(field.name(), field.type(), field.typeName(), field.length(), field.precision()))
poly_fields.append(QgsField("面积", QVariant.Double, "double", 20, 2))
poly_layer = QgsVectorLayer(f"Polygon?crs={crs.authid()}", f"{file_name}_封闭面", "memory")
pr_poly = poly_layer.dataProvider()
pr_poly.addAttributes(poly_fields)
poly_layer.updateFields()
poly_feats = []
for poly in polygons:
f = QgsFeature()
f.setGeometry(poly)
area = area_calc.measureArea(poly)
point = poly.pointOnSurface()
nearest_ids = src_index.nearestNeighbor(point.asPoint(), 1)
attrs = []
if nearest_ids:
nearest_feat = src_features[nearest_ids[0]]
attrs = nearest_feat.attributes()
else:
attrs = [None] * len(src_fields)
attrs.append(area)
f.setAttributes(attrs)
poly_feats.append(f)
pr_poly.addFeatures(poly_feats)
poly_layer.updateExtents()
QgsProject.instance().addMapLayer(poly_layer)
renderer = poly_layer.renderer()
symbol = renderer.symbol()
if symbol:
symbol.setColor(QColor(255, 0, 0, 80))
symbol.setOpacity(0.5)
sl = symbol.symbolLayer(0)
sl.setStrokeColor(QColor(255, 0, 0))
sl.setStrokeWidth(0.8)
poly_layer.triggerRepaint()
log_func(f" 已添加封闭面图层：{poly_layer.name()}，共 {polygon_count} 个面。")
elapsed = time.time() - start_time
log_func(f" 处理成功，耗时 {elapsed:.2f} 秒。")
return {
'break_count': break_count,
'polygon_count': polygon_count,
'midline_count': midline_count,
'elapsed': elapsed,
'file': file_path
}
# ------------------------------------------------------------
# 批量处理对话框
# ------------------------------------------------------------
class BatchDialog(QDialog):
def __init__(self, parent=None):
super().__init__(parent)
self.setWindowTitle("河道管理线批处理工具 - 嘉峪关市审计局 孙永鑫")
self.resize(600, 500)
layout = QVBoxLayout()
# 目录选择行
dir_layout = QHBoxLayout()
dir_layout.addWidget(QLabel("源目录："))
self.dir_edit = QLineEdit()
dir_layout.addWidget(self.dir_edit)
btn_browse = QPushButton("浏览...")
btn_browse.clicked.connect(self.browse_dir)
dir_layout.addWidget(btn_browse)
layout.addLayout(dir_layout)
# 选项复选框
self.check_midline = QCheckBox("生成河道中线层")
self.check_breakpoints = QCheckBox("生成断点标注层")
self.check_repaired = QCheckBox("生成管理线修复层")
self.check_polygon = QCheckBox("生成河道转面层（必选）")
self.check_polygon.setChecked(True)
self.check_polygon.setEnabled(False) # 不可取消
layout.addWidget(self.check_midline)
layout.addWidget(self.check_breakpoints)
layout.addWidget(self.check_repaired)
layout.addWidget(self.check_polygon)
# 运行按钮
self.btn_run = QPushButton("开始批处理")
self.btn_run.clicked.connect(self.run_batch)
layout.addWidget(self.btn_run)
# 日志文本框
self.log_text = QTextEdit()
self.log_text.setReadOnly(True)
font = QFont("Courier New", 9)
self.log_text.setFont(font)
layout.addWidget(self.log_text)
self.setLayout(layout)
def browse_dir(self):
dir_path = QFileDialog.getExistingDirectory(self, "选择包含河道线文件的根目录")
if dir_path:
self.dir_edit.setText(dir_path)
def log(self, message):
self.log_text.append(message)
self.log_text.ensureCursorVisible()
def run_batch(self):
root_dir = self.dir_edit.text().strip()
if not root_dir or not os.path.isdir(root_dir):
self.log("错误：请选择有效的目录。")
return
self.btn_run.setEnabled(False)
self.log("=========================================")
self.log("开始批量处理...")
self.log(f"根目录：{root_dir}")
# 收集所有矢量文件（常见扩展名）
extensions = ['*.shp', '*.gpkg', '*.geojson', '*.kml', '*.kmz', '*.gml']
all_files = []
for ext in extensions:
all_files.extend(glob.glob(os.path.join(root_dir, '**', ext), recursive=True))
# 去重（可能大小写不同）
all_files = list(set(all_files))
if not all_files:
self.log("未找到任何矢量文件。")
self.btn_run.setEnabled(True)
return
self.log(f"共找到 {len(all_files)} 个文件。\n")
# 选项
options = {
'gen_midline': self.check_midline.isChecked(),
'gen_breakpoints': self.check_breakpoints.isChecked(),
'gen_repaired': self.check_repaired.isChecked(),
'gen_polygon': True # 必选
}
success_count = 0
fail_count = 0
total_break = 0
total_poly = 0
total_start = time.time()
for f in all_files:
# 跳过明显不是线图层的文件（快速检查）
try:
info_layer = QgsVectorLayer(f, "tmp", "ogr")
if not info_layer.isValid() or info_layer.geometryType() != QgsWkbTypes.LineGeometry:
self.log(f"跳过非线文件：{f}")
continue
except:
self.log(f"跳过无法识别的文件：{f}")
continue
stats = process_single_file(f, options, log_func=self.log)
if stats:
success_count += 1
total_break += stats['break_count']
total_poly += stats['polygon_count']
else:
fail_count += 1
total_elapsed = time.time() - total_start
self.log("\n============= 批处理完成 =============")
self.log(f"总计文件：{len(all_files)}，成功：{success_count}，失败：{fail_count}")
self.log(f"累计检测断点：{total_break} 个，生成封闭面：{total_poly} 个")
self.log(f"总耗时：{total_elapsed:.2f} 秒")
self.log(f"开发者：嘉峪关市审计局 孙永鑫")
self.btn_run.setEnabled(True)
# ------------------------------------------------------------
# 启动对话框
# ------------------------------------------------------------
def main():
dlg = BatchDialog()
dlg.exec_() # 模态运行，对话框不会消失
main()
