1. 为什么我推荐用 PSQ 处理 PostGIS 拓扑分析做空间数据的人应该都有这种体会PostGIS 的空间查询和空间分析能力非常强但一旦涉及拓扑比如检查多边形是否共享边界、构建路网连通关系、判断要素之间是否重叠或存在缝隙SQL 写起来就会变得特别绕。拓扑逻辑本身不复杂复杂的是把点、线、面之间的关系统一管理起来还要保证数据在编辑过程中始终保持一致。我最早做行政区划数据融合项目时需要频繁判断地物之间的邻接关系最常写的就是一串嵌套子查询先查出边缘几何再判断哪些面共享这条边最后还要去重、去自相交。这种 SQL 不是不能写但写出来以后维护成本很高换一个人来看基本读不懂。后来接触到 PostgreSQL 自带的拓扑扩展也就是postgis_topology底层已经帮我们管理好了节点、边、面之间的完整拓扑关系但直接用 SQL 去操作这些内部表依然很繁琐——而且容易一不小心就把拓扑结构搞坏。PSQ 这个开源项目解决的就是这个痛点。它本质上是一个针对 PostGIS Topology 的 Python 封装工具把拓扑库里的节点、边、面、图层这些概念转换成 Python 对象让我们可以通过代码去建拓扑、往拓扑里加要素、检查错误、查询邻接关系而不是整天面对一长串 SQL。对习惯用 Python 做数据处理的人来说用 PSQ 要比直接撸原生拓扑 SQL 顺手得多。这篇就当作一份 PSQ 使用教程来写从环境准备、核心概念、常用操作到实战案例一步一步展开同时把我实际使用中踩过的坑一并记录下来。如果你正在做地理围栏、路网分析、地块合并、地籍管理这类需要严格拓扑关系的项目这篇应该能帮到你。1.1 拓扑分析和普通空间分析本质上的差异很多刚接触 GIS 的人会把“空间分析”和“拓扑分析”混在一起觉得都是和几何打交道。其实两者关注点完全不同。普通的空间分析比如计算两个多边形是否相交、求缓冲区、算面积关心的是“几何本身”的坐标和形状而拓扑分析关心的是“几何对象之间的连接关系”——哪些边是共用的、哪些节点连接了哪些边、哪些面围成了闭合区域。举例来说一张地籍图中相邻两块宗地应该共享同一条边界线而不是各自画一条重叠的线。用普通空间数据模型存两块宗地的边界可能看起来是重合的但在数据库里其实是两条完全独立的线只是坐标值恰好接近。一旦数据更新一块地挪了边界另一块地不会自动跟着变中间就会产生缝隙或重叠。拓扑模型通过把边界抽象为统一的“边”让两块地引用同一条边从底层保证数据不会出现这种不一致。PostGIS 的拓扑扩展正是按这个思路设计的节点表管理所有顶点边表管理弧段面表管理由边围成的闭合区域。这三类元素之间的引用关系构成了完整的拓扑网络。理解了这个模型后面用 PSQ 操作时就很容易上手。1.2 PSQ 在整个技术栈里的定位PSQ 并不是要替代 PostGIS也不是替代 GeoPandas、Shapely 这类纯 Python 空间分析库。它的定位非常明确用 Python 的方式去操作 PostGIS 的拓扑模型。它的工作方式是在 Python 端建立数据库连接后通过 SQLAlchemy 和底层拓扑表进行交互。你不需要自己去写 INSERT 语句往edge_data表里插入记录也不需要手工维护next_left、next_right这些拓扑指针PSQ 会把简单的 API 暴露给你内部该调用哪个拓扑函数、该更新哪些元数据它自己处理。从项目集成角度看PSQ 适合已经使用 PostgreSQL 做空间数据管理、同时又希望把拓扑处理逻辑放进 Python 数据处理流水线的团队。它不是一个学习 GIS 概念的入门工具而是配合专业空间数据库使用的实操工具。2. 环境准备数据库和 Python 两端缺一不可开始用 PSQ 之前数据库侧和 Python 侧的环境都要配齐。很多新手第一次跑不通过不是代码写错而是数据库里压根没有开启拓扑扩展或者 Python 依赖装的版本不对。2.1 数据库侧启用 PostGIS 和拓扑扩展PSQ 的底层能力全部来自 PostGIS 的拓扑模块所以在 PostgreSQL 里必须先安装 PostGIS 并启用扩展。我用的是 PostgreSQL 14、PostGIS 3.x这里的操作步骤在更早版本上基本也适用。CREATE EXTENSION IF NOT EXISTS postgis; CREATE EXTENSION IF NOT EXISTS postgis_topology;这两条 SQL 执行完后数据库里会多出topology这个 schema里面存放着拓扑元数据表topology表记录每个拓扑的基本信息layer表记录图层信息node、edge_data、face表分别存节点、边和面。另外还有relation表负责将不同图层上的几何要素映射到拓扑元素上。验证是否安装成功可以执行SELECT name, srid, tolerance FROM topology.topology;如果能正常返回空结果集说明扩展已经可用只是还没有创建任何拓扑。如果这条查询报错说 relation 不存在多半是拓扑扩展没装成功。2.2 Python 侧安装 PSQ 及其依赖PSQ 本身是一个 Python 包建议在虚拟环境里安装。它的核心依赖包括 SQLAlchemy、GeoAlchemy2 和 psycopg2分别负责 ORM 映射、空间数据类型的支持和 PostgreSQL 连接。pip install psq sqlalchemy geoalchemy2 psycopg2-binary如果你用的是 PostgreSQL 12 以下的版本建议确认一下 psycopg2 和 libpq 的兼容性PostgreSQL 14 以上的环境正常直接用 psycopg2-binary 即可。注意不要把 psycopg2 和 psycopg2-binary 同时装两个混在一起容易出现一些莫名其妙的环境冲突。在部分环境下如果安装psq拉不到包可能是你的 PyPI 源没有同步到位换国内镜像源后重试即可基本都能正常装上。2.3 验证 Python 能连上拓扑模块环境装好后写一个小脚本验证连接。先建一个engine再让 PSQ 去读取数据库里已有的拓扑列表。我这种情况下建议先手动创建一个测试拓扑把环境链路整体跑通。SELECT topology.CreateTopology(demo_topo, 4326, 0.00001);然后在 Python 里执行from sqlalchemy import create_engine from psq import TopologyManager engine create_engine(postgresql://postgres:postgreslocalhost/gisdb) manager TopologyManager(engine) for topo in manager.topologies: print(topo.name, topo.srid, topo.tolerance)正常输出类似demo_topo 4326 1e-05说明数据库连接没问题拓扑扩展可用PSQ 也正确读取到了拓扑信息。到这一步最基础的运行环境就搭好了。3. 核心概念先把底层的表和对象搞清楚很多人用不好 PSQ 或者说用不好 PostGIS 拓扑根子在于不太理解底层的数据结构。PSQ 只是把底层结构包装成 Python 对象如果你不清楚节点、边、面分别对应什么代码照样写不明白。3.1 节点、边、面三张核心表一个拓扑里最核心的数据是三张表node、edge_data和face。节点表node记录拓扑中的顶点。每个节点有一个全局唯一的node_id以及对应的几何坐标。边表edge_data记录两个节点之间的弧段除了edge_id和几何外还记录左面 ID、右面 ID、起点、终点以及拓扑指针信息。这些指针是拓扑数据结构的关键它们把边串成环用来闭合面。面表face记录由边围成的闭合区域。注意一点任何一个拓扑中都有一个“宇宙面”代表整个拓扑覆盖范围之外的空间它始终存在面 ID 通常是 0。图层layer表则把一个拓扑里的要素按业务分组比如同一张表里的道路数据和地块数据可以分到不同图层。每个图层的要素通过relation表映射到具体的拓扑元素上。3.2 PSQ 中你会频繁碰到的几个对象在 PSQ 的封装下你不需要直接跟数据库表打交道而是通过对象去操作。我用过程中最常用的是下面这几个TopologyManager负责数据库连接管理多个拓扑实例。Topology对应一个具体拓扑可以查节点、边、面也可以给拓扑添加要素。Node、Edge、Face分别对应拓扑里的节点、边、面对象带有几何信息和拓扑属性。Layer对应图层关联到具体的业务表。读取拓扑下所有节点可以这样写topo manager.get_topology(demo_topo) nodes topo.nodes() for node in nodes.limit(10): print(node.node_id, node.geom)如果只是统计边和面的数量类似这样edge_count topo.edges().count() face_count topo.faces().count() print(fedges: {edge_count}, faces: {face_count})这些方法内部会映射到底层的SELECT查询所以如果你熟悉 SQLAlchemy 的 Query 对象会发现整个调用链非常顺畅。3.3 注意 schema 和 search_path 问题我在第一次用 PSQ 的时候踩过一个很低级的坑数据库里明明有拓扑数据但 Python 端查询时报 table not found。排查到最后发现是数据库连接串里的search_path没有包含topologyschema导致 SQLAlchemy 生成查询时找不到表。解决办法是在连接时显式指定 schemaengine create_engine( postgresql://postgres:postgreslocalhost/gisdb, connect_args{options: -c search_pathpublic,topology} )如果你在多个拓扑之间切换也可以先用SET search_path明确当前会话要用的 schema这样 PSQ 再执行内部查询时就不会迷路。4. 创建一个新拓扑并导入要素数据环境通了、概念也清楚了接下来就该动手建一个真实可用的拓扑。这里我用一个最简单的场景把两块相邻的多边形导入到一个新建拓扑里演示完整的操作流程。4.1 用 PSQ 创建拓扑在 SQL 里创建拓扑通常用CreateTopology函数PSQ 同样提供了对应的接口。参数需要指定拓扑名、坐标系 SRID 和容差值。manager.create_topology( nameland_topo, srid4326, tolerance0.00001 )容差值tolerance是拓扑建模里最重要的参数它决定坐标为多少范围内的两个点会被视为同一个点。0.00001在经纬度坐标系下大约相当于一米左右适合地块级别的数据。如果你处理的是高精度测量数据容差值可能要更小如果是大范围路网数据容差值可以适当放大。创建完后再用get_topology拿到这个拓扑对象后续操作都基于它。topo manager.get_topology(land_topo) print(topo.name, topo.srid, topo.tolerance)4.2 向拓扑添加要素把一块带几何的要素加入到拓扑中PSQ 通常提供与 PostGIS 拓扑函数对应的方法。以内置函数TopoGeo_AddPolygon为例它可以接收一个多边形几何自动将几何的边界切成拓扑边并创建对应的面。from shapely.geometry import Polygon # 两个相邻地块的边界坐标故意让中间线不完全重合 poly_a Polygon([ (120.0, 30.0), (120.1, 30.0), (120.1, 30.1), (120.0, 30.1) ]) poly_b Polygon([ (120.1, 30.0), (120.2, 30.0), (120.2, 30.1), (120.1, 30.1) ]) face_a topo.add_polygon(poly_a) face_b topo.add_polygon(poly_b)PSQ 在把 Shapely 几何转换成数据库里的 PostGIS 几何时会自动做坐标转换和格式匹配。如果传入的几何不是合法的多边形比如有自相交或太碎的边界这一步会直接报错。所以在大批量导入前务必先用 Shapely 的is_valid方法做一些前置过滤。4.3 检查导入结果导入完成后可以分别查看节点、边和面的数量确认拓扑结构是否建立起来。print(topo.nodes().count()) print(topo.edges().count()) print(topo.faces().count())两个相邻矩形一共会有 6 个节点每个矩形 4 个角重叠边共享 2 个角所以合计 6 个节点。边一共 7 条其中外边界 6 条、共享边界 1 条。面一共 3 个2 个实际地块加上 1 个宇宙面。如果你的结果和这个对不上比如边数变成了 8 条那基本可以断定是容差值设置得不够合理共享边界没有被正确合并两个矩形被当成了完全独立的空间。4.4 对已有图层做拓扑化处理实际项目里很少从一个空拓扑开始更多情况是已经有一张业务表里面存了几千几万个面要素想把这些要素变成拓扑结构。这种情况需要用toTopoGeom函数把普通几何转换成拓扑几何并写入指定图层。在 PSQ 中如果 API 直接支持图层式导入最好如果不支持可以回退到执行原生 SQLSELECT topology.TopoGeo_AddPolygon( land_topo, ST_GeomFromText(POLYGON((120.0 30.0, 120.1 30.0, 120.1 30.1, 120.0 30.1, 120.0 30.0)), 4326) );我的建议是大批量导入时不要一条条在 Python 里循环调用那样效率太低。正确做法是先把业务表里的要素批量推到一个临时表再用一条 SQL 把整个表注册到拓扑图层里这样数据库可以内部高效批量处理。5. 从实际问题出发如何查询相邻地块拓扑数据建好之后最常用到的查询就是“某个面的邻居是谁”。比如在土地分析里我需要知道某块宗地和周围哪些地块接壤在传统空间索引下要跑相交查询性能还不一定好在拓扑模型里因为共享边是同一份数据查询邻居只需要根据面 ID 找到它引用的边再根据边找到另一侧的面即可。5.1 根据面 ID 查询边一个面由多条边围成边表里记录了left_face和right_face两个字段。如果某个面的 ID 出现在某条边的左面字段里那这条边就是该面的边界边另一侧的面就是它的邻居。在 PSQ 里可以这样组织查询逻辑face_id 1 boundary_edges topo.edges().filter( (Edge.left_face face_id) | (Edge.right_face face_id) ) for edge in boundary_edges: neighbor_face edge.left_face if edge.right_face face_id else edge.right_face if neighbor_face ! 0: print(neighbor face:, neighbor_face)注意neighbor_face等于 0 的情况说明这条边没有邻居而是整个拓扑区域的外边界。5.2 获取邻居面的几何查到了邻居面 ID 后下一步自然是拿到邻居面的几何做面积计算或进一步分析。PostGIS 拓扑模块里ST_GetFaceGeometry就是专门干这个的。在 PSQ 里可以直接通过topo.face_geometry(face_id)获取对应的 Shapely 几何for face_id in neighbor_ids: geom topo.face_geometry(face_id) print(face_id, geom.area)如果只是想看拓扑内部结构topo.face_geometry()返回的是一个 MultiPolygon 或 Polygon可以直接转成 GeoJSON 做后端预览。5.3 分析结果的正确性验证每次都依赖拓扑关系查询最怕的是拓扑本身已经坏掉。比如边的left_face和right_face指向不存在的面 ID或者面几何被重复写入都会让查询结果出现脏数据。我通常会在查询前后做一次完整性校验用 PostGIS 自带的拓扑检查功能SELECT * FROM topology.ValidateTopology(land_topo);如果有问题会返回具体错误类型和相关的节点、边 ID。通过这个反馈再去定位并修复数据比凭感觉检查要高效得多。6. 工程化落地批量性能优化和避坑记录如果只是处理几十个面PSQ 怎么用都不会有大问题。一旦数据量上万性能、事务、术语约定这些工程化问题就会浮现出来。这一章把我踩过的一些坑和沉淀下来的优化经验写清楚。6.1 批量导入时不要在 Python 里循环最开始做批量导入时我直接写了个 for 循环把几万个多边形逐一调用topo.add_polygon()提交到数据库。结果速度慢到无法忍受每个多边形要经过 Python 对象构造、SQLAlchemy 查询、PostGIS 拓扑计算、事务提交这一整套流程。后来改成批量处理的方式先用 SQL 把数据整理到临时表再一次性执行 TopoGeo_AddPolygon 的批处理版本。数据量上万时这种做法的耗时只有循环方式的几十分之一。批量导入的基本流程是在数据库里建一张临时表结构包含 ID 和几何字段。用 Python 将数据批量写入这张临时表。用一条 SQL 调用TopoGeo_AddPolygon或toTopoGeom把临时表里的几何批量注册到拓扑中。导入完成后把relation表里对应的映射关系同步到业务图层。PSQ 的意义更多在后续查询和维护环节而不是替代数据库完成高频批量计算。6.2 事务和提交时机PSQ 依赖 SQLAlchemy 管理事务每一次拓扑操作都涉及多个底层表的变更。如果事务提交时机不对拓扑很可能会处于中间状态——节点加了边还没加完整面也没闭合。我建议把每批操作放在一个事务里全部做完后再统一提交。SQLAlchemy 的基本用法with engine.begin() as conn: # 在这里通过连接执行批量导入操作 pass如果中间某一步出错整个事务回滚拓扑对象不会被破坏。这个习惯能帮我避免很多“数据调到一半库里留下半成品”的尴尬情况。6.3 容差值设置不当引发的数据合并容差值太大会导致两个本不该合并的节点被合并比如相邻地块之间一个小小的转角可能会被吞掉导致边界走样容差值太小则会导致大量节点没被识别成同一个点共享边断成两半拓扑关系建立不起来。实际项目里容差值不是拍脑袋定的通常要和数据的源精度对齐。比如原始数据精度是 0.1 米容差值设置 0.05 米就比较合理如果是经纬度坐标要考虑不同纬度下一度对应的实际距离再换算成合适的容差值。必要的时候先在小范围内测试不同容差值对应的拓扑结果再全量执行。6.4 明确图层归属避免不同业务交错污染一个拓扑可以被多个图层共用路网数据一个图层、地块数据另一个图层这没问号。但如果每层都添加了互相覆盖的数据一旦某层做更新可能会影响到另一层的拓扑结构。我强烈建议不同业务数据不要混在同一个拓扑里。宁可多建几个拓扑也不要省事放在一起。原因很简单拓扑的维护是整体性的一个图层的数据变更可能导致整个拓扑网络重算影响范围不可控。分开拓扑至少一个拓扑坏了不会拖累全部业务。7. 用 SQLAlchemy 集成到现有数据管道PSQ 本身建立在 SQLAlchemy 之上这意味着它可以比较容易地嵌入到现有的 Python 数据管道中。如果你的项目已经用了 SQLAlchemy那么集成 PSQ 基本是零成本如果还没有也不必为了 PSQ 专门引入整套 ORM直接用它提供的接口即可。7.1 把 PSQ 的查询结果转成 GeoDataFrame很多空间分析最终都要落到 GeoPandas 上做可视化或后续计算。PSQ 查询出来的 Shapely 几何可以直接转成 GeoDataFrame不需要经过 GeoJSON 中转。import geopandas as gpd from shapely.geometry import mapping rows [] for face_id in face_ids: geom topo.face_geometry(face_id) rows.append({face_id: face_id, geometry: geom}) gdf gpd.GeoDataFrame(rows, crsEPSG:4326)这种方式在处理大量面要素时非常实用拓扑里算好的邻接关系可以直接跟业务属性关联再落到 GeoDataFrame 里做展示或进一步空间运算。7.2 配合 Alembic 做数据库版本管理如果你的项目已经用 Alembic 管理数据库结构变更建议把拓扑的创建过程也纳入迁移脚本。比如在迁移脚本里创建拓扑、注册图层这样部署到新环境时不需要手动在数据库里执行一遍 SQL。我实际执行时会写两个步骤先创建扩展再创建拓扑。扩展属于数据库级别操作通常放在初始化迁移里拓扑创建属于业务逻辑放在具体业务表的迁移里。这样职责清晰后期排查问题也容易定位。7.3 监控拓扑数据质量最后一个建议来自我自己的教训给拓扑加一个定期的质量巡检任务。用ValidateTopology函数扫描所有拓扑结构把异常结果写到另一张监控表里配合定时任务一旦发现拓扑错误就能第一时间收到通知。空间数据的质量问题往往是累积的小错误发现得越晚修复成本越高。尤其在使用 PSQ 这类工具简化了操作以后更需要在自动化层面兜底保证底层数据的健康状态。根据我的实际使用体验PSQ 最大的价值在于把 PostGIS 拓扑的操作门槛降了下来让我可以用熟悉的 Python 对象模型去理解和维护拓扑关系。但在真正投入生产之前还是建议先把 PostGIS 拓扑的底层原理吃透——工具替你省掉的是繁琐的 SQL 拼装过程省不掉的是对拓扑模型本身的理解和敬畏。