ARTICLE DETAIL

资讯详情

深耕郑州网站建设与运营推广的一线实战洞察。

QGIS Python脚本批量处理矢量数据:字段汇总、线段合并与坐标系实战

QGIS Python脚本批量处理矢量数据:字段汇总、线段合并与坐标系实战 做GIS的兄弟基本都遇到过这种尴尬同一套操作一个图层来一次能忍一百个图层还要一个个点鼠标那就真的坐不住了。QGIS里自带Python控制台就是拿来干这种事的。这篇文章是QGIS系列的第四十三篇核心就一件事怎么用Python在QGIS里批量、稳定、不返工地处理数据。不管你是刚接触脚本的小白还是已经在写插件的二次开发党环境配置、字段汇总、线段合并、可视化、批处理这块的内容都可以直接抄作业。你不需要一上来就啃QGIS源码先把Python这层窗户纸捅破日常工作至少快一半。1. 环境先行QGIS与Python怎么搭最顺手1.1 用自带Python控制台还是外部解释器很多人第一次接触PyQGIS是从“插件—Python控制台”开始的快捷键CtrlAltP左下角弹出一个交互窗口。QGIS 3.x版本内置的Python环境是经过打包隔离的目录和系统Python互不干扰所以你在系统里用conda装多少包都不影响QGIS控制台。反过来也一样控制台里pip list看到的库和命令行看到的可能是两套东西。内置控制台适合干三件事临时跑一段逻辑、调试图层处理结果、写短平快的小工具。我实际用得最多的场景是“选中一个图层写两行代码算点东西”比如统计字段和、按属性筛要素直接在控制台里回车就出结果比新建工程、拖算法模块快得多。但内置控制台有个毛病代码长了不好维护。你写了个几十行的脚本一旦缩进乱掉或者逻辑复杂起来交互式窗口体验就很差。我的建议是两种方式结合日常小操作直接在控制台里写一两行能解决的问题不要去开IDE。超过二十行、或者要反复调整参数的逻辑用外部编辑器写好再通过控制台的“加载脚本”按钮执行。1.2 PyCharm和VSCode对接QGIS环境不少朋友想用PyCharm或VSCode写PyQGIS脚本经常卡在import qgis这步。最省事的办法是直接用QGIS自带的Python解释器而不是系统里那个干净Python。Windows下QGIS安装目录的bin文件夹里通常有个python-qgis.bat它会帮你设置好所有QGIS相关的环境变量包括PYTHONPATH和QGIS_PREFIX_PATH。在PyCharm里如果你的解释器列表里能选到这个bat对应的Python环境就直接选如果选不上就手动创建一个解释器指向C:\Program Files\QGIS 3.28\bin\python.exe然后到环境变量里补两个关键项PYTHONPATH指向QGIS安装目录下的python目录以及python\plugins目录。QGIS_PREFIX_PATH指向QGIS安装根目录。VSCode的配置思路同理在.vscode/settings.json里设置python.defaultInterpreterPath再给终端配置环境变量。但说实话Windows上最稳定的还是那招——用python-qgis.bat启动脚本把你的代码从外部文件传进去。比如python-qgis.bat D:/project/qgis_script.py这样脚本里import qgis一定不会报模块找不到的错误。如果是Linux环境尤其是国产操作系统上部署QGIS开发环境原理一样只不过PYTHONPATH一般指向/usr/lib/qgis/python具体路径可以通过qgis --version或find命令来确认。1.3 CGCS2000坐标系在QGIS里的正确姿势关于“QGIS软件如何设置2000坐标系”是很多做国土、规划项目的同行问得最多的问题。CGCS2000是国家大地坐标系EPSG编码一般是4490地理坐标系但日常使用更多是带分带的高斯-克吕格投影比如3度分带经常用EPSG:4547或4548、6度分带对应EPSG:4491等。设置路径有三种项目层面菜单“项目—属性—坐标系”搜索4490或对应的投影EPSG码全局生效。图层层面右键图层“属性—源—设置坐标系”只改当前图层。快速切换右下角EPSG按钮点一下输入编号搜索回车即切换。用Python设置坐标系也很直接crs_2000 QgsCoordinateReferenceSystem(EPSG:4490) layer.setCrs(crs_2000)但注意setCrs只是强行修改图层的坐标系定义并不会把坐标数值重新算一遍。如果原始数据其实是WGS84经纬度你直接把图层CRS改成CGCS2000图面上看起来好像没问题一旦导出或者叠加其他数据就会偏。正确的做法是先确认原始坐标系再使用QgisCoordinateTransform或Processing里的重投影算法做转换。设置坐标系这事最怕的不是不会操作而是把“改定义”和“重投影”两件事混在一起。2. Python处理矢量数据的三个高频动作2.1 图层读取与字段汇总字段汇总是日常操作里频率极高的需求。如果只是加总一个字段用界面上的字段计算器或Basic Statistics工具也能做但Python的优势是可以把“汇总”和“下游处理”串起来一步到位。假设我手上有一个全国地级市人口分布SHP字段名是pop要快速算总人口layer iface.activeLayer() if not layer.isValid(): print(图层无效) total sum(feature[pop] for feature in layer.getFeatures()) print(f总人口: {total})这里有个细节feature[pop]这种写法的前提是字段名在图层里真实存在且pop字段的类型是整型或浮点型否则会拿到Nonesum直接报错。更稳妥的方式是加一层判断total 0 for feature in layer.getFeatures(): value feature[pop] if value is None: continue total float(value)如果要对某个分组字段做汇总比如按省汇总各地级市人口没必要自己写循环字典QGIS内置的聚合计算类更快from qgis.core import QgsAggregateCalculator calc QgsAggregateCalculator(layer) result calc.calculate(pop, QgsAggregateCalculator.Sum, province)返回结果是一个QgsAggregateComputation对象可以从中拿到分组键和值。但坦白讲聚合功能更强的是Processing里的statisticsbycategory算法你只需要指定分组字段和统计字段它会把均值、最大、最小、总和、标准差一次性全算出来。Python调用这个算法的写法是import processing processing.run(native:statisticsbycategory, { INPUT_LAYER: layer, CATEGORIES: [province], FIELDS: [pop], OUTPUT: D:/result_stat.csv })2.2 数据定义覆盖不只是界面功能搜“qgis 的数据定义覆盖”的人一部分是在图层样式面板里迷了路一部分是想用Python控制动态样式。数据定义覆盖是QGIS里特别强大但常被低估的功能。直白讲它允许你把某个属性值从“固定值”改成“一个字段或表达式动态生成的值”。比如做人口专题图你想让点符号大小随人口数量变化传统做法是先用字段计算器生成一个新字段再按字段分级。有了数据定义覆盖你可以直接在符号化面板里点击“大小”那一行右侧的蓝色/灰色小方块选择“按表达式”输入sqrt(pop) / 1000符号就会实时变化不需要新增字段。Python里动态设置数据定义覆盖的关键是拿到渲染器的符号层再用setDataDefinedProperty绑定属性from qgis.core import QgsSymbolLayer, QgsProperty renderer layer.renderer() symbol renderer.symbol() symbol_layer symbol.symbolLayer(0) symbol_layer.setSizeMode(QgsSymbolUnitType.RenderedSize) # 按图面尺寸控制 symbol_layer.setDataDefinedProperty( QgsSymbolLayer.PropertySize, QgsProperty.fromExpression(sqrt(\pop\) / 1000) ) layer.triggerRepaint()这段代码最常用的场景是点图层。如果你做的是面的透明度、线的宽度只需把PropertySize换成PropertyFillColor、PropertyStrokeWidth等对应枚举。踩坑提醒数据定义覆盖的表达式里如果字段名是中文记得用双引号包起来且表达式里不能再叠加双引号否则解析会报错。另外改了数据定义后一定要triggerRepaint()否则画布不会刷新你会以为代码没生效。2.3 多条线段合并成一条的实现方案“QGIS并行多条线段合并成一条线段”这个需求我在路网数据清洗里遇到过好多次。比如一段国道被行政边界切成几十段想把它们拼成一条完整路径。最简单的方案是Processing的mergelines算法调用方式import processing processing.run(native:mergelines, { INPUT: layer, OUTPUT: D:/merged_line.shp })这个算法会把所有输入线的节点串起来但它有一个隐含逻辑如果多条线段首尾不连续或者存在分支结果可能不是你想要的单线。更麻烦的是当线段方向混乱时合并后会出现“回头路”几何形状像一团乱麻。这时候就得手动控制合并顺序。我的做法是写一个处理函数读取所有要素的节点序列。按照起点和终点的相邻关系排序保证下一条线的起点尽量接近上一条线的终点。把节点坐标合并成一条QgsLineString。from qgis.core import QgsLineString, QgsGeometry, QgsFeature def merge_lines_to_single(lines): coords [] for geom in lines: g geom.constGet() if g.isMultipart(): for part in g.parts(): coords.extend(part.asPointList()) else: coords.extend(g.asPointList()) return QgsGeometry.fromPolylineXY(coords)这个方法的局限性是如果线段之间存在分叉或环状连接简单拼接会丢失拓扑关系。遇到河流支流、道路互通这种数据应该先做分叉处理或使用网络分析模块而不是强行合并。所以遇到这类需求先问自己一句我要的是“物理上连成一条线”还是“逻辑上是一条路径”前者可以用mergelines后者要用路径规划算法。3. 把属性数据拿出来做分析和可视化3.1 类型转换与数据清洗在QGIS里用Python处理数据最花时间的往往不是几何操作而是属性表里的脏数据。CSV导进来的数字可能带着空格Excel转出来的整型可能被读成浮点空值在不同来源里可能是空字符串、None、NULL字符串或者na。这些都要在进入分析前统一处理。PyQGIS里读取字段值时推荐用feature.attribute(field_name)它等价于feature[field_name]但可以传默认值value feature.attribute(area) or 0 # 空值回退为0 area float(value)不过要注意or会把0也当作空值处理所以更严谨的写法是判断Noneraw feature.attribute(area) area float(raw) if raw is not None else 0.0字段名称里的空格和特殊字符也要留个心眼。从Excel转过来的字段名可能有空格在QGIS里显示正常但在Python脚本里访问时可能会踩坑我习惯在脚本开头做一次字段名标准化for field in layer.fields(): new_name field.name().strip().replace( , _) if new_name ! field.name(): layer.renameAttribute(field.name(), new_name)字段重命名会导致外部引用的字段失效所以这个操作要在数据探索阶段做等确定字段名之后再写正式处理脚本。3.2 matplotlib画图时横坐标太密怎么办很多人喜欢把QGIS属性数据导出来再用matplotlib画图。最经典的翻车现场是横轴明明是日期或者序号画出来密密麻麻一条黑线刻度标签重叠成一根棍子。核心原因是matplotlib默认会自动给每个数据点生成刻度数据量一上去就炸。解决办法是主动设置刻度间隔import matplotlib.pyplot as plt import matplotlib.dates as mdates fig, ax plt.subplots(figsize(12, 5)) ax.plot(dates, values) # 每5个刻度显示一个 ax.xaxis.set_major_locator(mdates.DayLocator(interval5)) # 或者用 MultipleLocator 处理纯数字横轴 from matplotlib.ticker import MultipleLocator ax.xaxis.set_major_locator(MultipleLocator(10)) # 标签旋转避免重叠 plt.xticks(rotation45) plt.tight_layout() plt.show()如果横轴是字符串字段比如五十个行政区名称我会在画图前直接切片挑选关键点显示ticks list(range(0, len(names), 10)) ax.set_xticks(ticks) ax.set_xticklabels([names[i] for i in ticks], rotation45)还有一个中文显示问题。Windows上matplotlib默认中文字体可能缺图里全是不认识的方块。推荐在脚本顶部直接设置plt.rcParams[font.sans-serif] [Microsoft YaHei] plt.rcParams[axes.unicode_minus] FalseLinux环境没有微软雅黑就换成Noto Sans CJK SC。4. 批处理与性能优化从脚本到小工具4.1 批量遍历图层文件实际项目里很少只处理一个图层更多是几十上百个SHP或者GeoPackage。这个时候你会真正体会到“写一次循环省一周加班”的快感。批量处理的标准模板大致是这样from glob import glob shp_list glob(D:/data/*.shp) output_dir D:/output for index, shp in enumerate(shp_list, 1): layer QgsVectorLayer(shp, tmp, ogr) if not layer.isValid(): print(f[跳过] {shp} 加载失败) continue # 这里放你的核心处理逻辑 print(f[{index}/{len(shp_list)}] 处理 {shp})有几个细节值得注意。第一SHP文件是多个子文件组成的集合.shp、.dbf、.shx必须放在一起拷贝时别只拷一个。第二输出文件名建议加上处理时间戳或者序号避免覆盖前一轮结果。第三如果循环内部要写新字段记得在startEditing和commitChanges之间完成否则修改不会落盘。批量处理最怕的是中间某一步崩溃前面全都白干。我的习惯是每处理完一个文件就用CSV记录一行日志包含文件名、状态、处理耗时、要素数量。后面排查问题的时候一眼就能看出是哪个文件出的问题不用翻控制台日志。4.2 多进程加速的取舍图层数量大、单文件数据量也大的时候逐条循环确实让人着急。一台八核电脑只用一个核跑怎么看都亏。那多进程是不是无脑上这里必须先说一个硬性限制PyQGIS的对象不能直接通过multiprocessing传给子进程。你尝试pool.map(worker, layers)大概率会收到pickle序列化错误。正确思路是每个子进程独立创建一个QGIS应用实例然后自己加载文件、处理文件、写出临时文件最后主进程再汇总。我用过一段时间的模板如下from multiprocessing import Pool from qgis.core import QgsApplication def process_one(shp_path): qgs QgsApplication([], False) qgs.initQgis() try: layer QgsVectorLayer(shp_path, tmp, ogr) # 处理逻辑 return {path: shp_path, count: layer.featureCount()} finally: qgs.exitQgis() if __name__ __main__: paths glob(D:/data/*.shp) with Pool(processes4) as pool: results pool.map(process_one, paths)多进程提升主要在CPU密集型计算上有效。比如对每个要素做复杂的几何运算、缓冲区计算、空间关系判断。如果你的瓶颈其实是磁盘读取或数据库连接多进程反而可能把磁盘IO拖垮速度不升反降。我的建议是先用单进程跑一遍用time记录耗时再用cProfile看热点函数。如果时间确实大量消耗在几何计算上再上多进程。如果主要耗时在读取文件优先考虑把SHP转成GeoPackage或者做空间索引而不是盲目开进程。5. 常见问题排查与这半年的实操心得5.1 环境安装与导入失败的排查PyQGIS最常见的问题是ModuleNotFoundError: No module named qgis。不管是在PyCharm、VSCode还是命令行里遇到基本就是PYTHONPATH没指对。排查思路可以按顺序来问题现象可能原因解决办法import qgis 报模块不存在PYTHONPATH未设置追加QGIS的python目录到sys.pathimport qgis.core 报DLL加载失败bin目录不在PATH把QGIS安装目录的bin加入PATHQGIS版本和Python版本不一致使用了系统Python改用python-qgis.bat或QGIS内置解释器控制台里能跑外部脚本不行缺少环境变量在脚本中设置QGIS_PREFIX_PATH插件安装后报错插件兼容问题确认QGIS版本或重新安装插件外部脚本最稳妥的启动模板我放在每次写脚本的头部防止现场翻车import sys from pathlib import Path if qgis not in sys.modules: qgis_prefix rC:\Program Files\QGIS 3.28 sys.path.append(str(Path(qgis_prefix) / python)) sys.path.append(str(Path(qgis_prefix) / python / plugins)) os.environ[QGIS_PREFIX_PATH] qgis_prefix from qgis.core import QgsApplication5.2 数据处理过程中容易踩的坑第一个坑是坐标系混乱。之前说过修改图层CRS定义和重投影是两回事。处理跨坐标系数据时我习惯先用QgsCoordinateReferenceSystem打印出两个坐标系的EPSG编码再用QgsCoordinateTransform做转换from qgis.core import QgsCoordinateReferenceSystem, QgsCoordinateTransform source_crs QgsCoordinateReferenceSystem(EPSG:4326) target_crs QgsCoordinateReferenceSystem(EPSG:4547) transform QgsCoordinateTransform(source_crs, target_crs, QgsProject.instance())第二个坑是字段名被截断。Shapefile的字段名最多10个字符从GeoDatabase转过来时字段名会被截断这时候脚本里写的字段名怎么都对不上。处理这种问题先打印所有字段列表print([field.name() for field in layer.fields()])第三个坑是几何无效。某些从第三方软件导出的数据几何拓扑有问题缓冲区或者相交计算直接报错或者结果为空。遇到这种情况先跑一遍修复几何算法processing.run(native:fixgeometries, {INPUT: layer, OUTPUT: fixed.gpkg})第四个坑是改了样式忘了刷新画布。脚本修改了渲染器或者数据定义覆盖后画布上没反应别怀疑代码先检查有没有调用iface.mapCanvas().refresh()或者layer.triggerRepaint()。5.3 Python控制台写脚本的一些习惯QGIS系列写了四十多篇中间踩了很多坑也总结出一些写PyQGIS脚本的小习惯这里分享给新人。第一先小步跑通再批量。处理一百个文件之前先用一个文件调试完整流程确认输出完全符合预期再套循环。第二中间结果及时落盘。连续操作三个小时跑了五个中间处理步骤最后发现源头一个参数错了如果没有中间文件只能全部重来。我一般会在每个关键节点输出一个带步骤编号的文件宁可多占点磁盘也绝不让前面算好的东西白费。第三学会用错误日志。QGIS控制台会输出Python异常信息但如果你用的是外部脚本最好把异常写进日志文件import traceback try: # 你的逻辑 pass except Exception: with open(D:/error.log, a, encodingutf-8) as f: f.write(traceback.format_exc())第四善用QgsProject.instance()。很多刚接触PyQGIS的新人搞不明白为什么QgsProject很重要。简单理解它管着当前工程的所有图层注册表。你手动创建的QgsVectorLayer如果不加进去界面上不会显示运行processing算法也找不到它。所以创建完图层后记得QgsProject.instance().addMapLayer(layer)写到这突然想起前几天一个朋友问我为什么我每次在QGIS里跑一段同样的Python脚本结果都和手动画的不一样聊到最后发现他手动操作时用的是当前选中要素而脚本里遍历的是全部要素。这种“选中范围”和“全部范围”的差异其实比任何语法错误都容易坑人。所以脚本开头我习惯打印一下执行范围print(f当前选中要素: {layer.selectedFeatureCount()}, 总要素: {layer.featureCount()})QGIS里用Python处理数据说到底不是拼API记得多熟而是先把数据本身想清楚。坐标系、字段类型、几何拓扑、要素范围任何一个前提变了脚本都可能跑出完全不同的结果。把重复劳动交给脚本把判断留给自己这才是用QGIS写Python处理数据最舒服的状态。
返回列表