30个PostGIS核心函数与pgRouting最短路径实战

做 GIS 开发这几年我越来越觉得 PostGIS 就是空间数据处理的地基。你可以在 MySQL 里存几个坐标点但只要一碰到“路网分析”“缓冲区计算”“最近邻查找”“最短路径规划”这类真需求最后基本都会回到地理空间数据库这套体系里来。尤其 PostGIS 配合 pgRouting一套 PostgreSQL 就能把空间存储、空间计算、路径规划全包了。这篇文章我不讲概念性的空话直接按我自己的使用经验把 30 个最常用的 PostGIS 核心空间函数拆开讲清楚再穿插 pgRouting 的最短路径和距离函数实战不管是应付考试、准备面试还是真要在项目里落地都能直接拿去参考。适合看这篇文章的人我猜大致有三类一是刚接触空间数据库、被各种函数搞晕的初学者二是已经用 PostGIS 做点查表业务但没深入做过路径分析的人三是需要在简历或项目里体现空间计算能力想把最短路径这块补上的开发者。无论哪类核心思路都是一样的先理解函数解决什么问题再记住使用场景最后上手跑一遍。1. 环境与安装把 PostGIS 装明白才算拿到入场券1.1 快速安装路线先说我个人最推荐的方式。除非你有特殊的离线部署要求否则本地开发一律用 Docker省掉一多半环境折腾。docker run -d --name gis-db \ -e POSTGRES_PASSWORDpostgres \ -e POSTGRES_DBgis \ -p 5432:5432 \ postgis/postgis:16-3.4这个镜像默认已经启用了 PostGIS 扩展进去之后直接执行CREATE EXTENSION IF NOT EXISTS postgis;如果是生产环境用 Ubuntu 部署通常用包管理器装sudo apt install postgresql-16 postgresql-16-postgis-3 postgresql-16-pgrouting装完同样先建库再启用扩展。注意 Windows 上别图省事一定要用 EnterpriseDB 的安装包配合 Stack Builder 选 PostGIS版本要和你安装的 PostgreSQL 大版本一致否则后面创建扩展时大概率会栽跟头。1.2 postgis 安装失败的几种典型原因这个坑太常见了我几乎每隔一阵就能在社区看到有人问“为什么 CREATE EXTENSION postgis 报错”。我自己踩过和帮人排查过的案例里无外乎下面几种版本不匹配。PostGIS 是跟着 PostgreSQL 大版本走的PostgreSQL 16 需要对应版本号的 PostGIS 3.4。用 apt 装的时候如果 PostgreSQL 的 apt 源没更新可能装到旧版 PostGIS然后扩展文件目录对不上。缺依赖库。PostGIS 编译安装需要 GEOS、GDAL、Proj 这些库少一个都会导致扩展无法加载。用官方包管理器一般会自动处理依赖但编译安装就要格外小心。扩展目录不对。报错信息里如果出现could not open extension control file十有八九是 PostgreSQL 的extension目录里没有 PostGIS 的控制文件。可以先检查SHOW shared_preload_libraries;再看扩展目录是否存在。在错误的数据库中执行。CREATE EXTENSION是在具体库上做的不是全局操作。有次同事在postgres默认库里建好了扩展结果连业务库之后一直提示找不到函数就是这个问题。解决之后可以用下面这段 SQL 做环境验证SELECT PostGIS_Full_Version(); SELECT name, default_version FROM pg_available_extensions WHERE name IN (postgis, pgrouting);如果能正常返回版本信息环境基本就绪。我习惯顺手跑一个最小测试SELECT ST_Distance( ST_SetSRID(ST_MakePoint(116.3, 39.9), 4326), ST_SetSRID(ST_MakePoint(116.4, 39.9), 4326) );这个结果会是一个以“度”为单位的数值不是米。刚接触的人很容易在这里被误导。真正的米制距离要么把坐标转投影坐标要么用geography类型。这个细节后面讲距离函数时会重点展开。2. 30 个核心空间函数一表读懂随查随用2.1 函数分类与完整清单我按自己的使用频率把 30 个函数分成四组构造与转换、空间关系判断、几何处理与聚合、测量计算。先给一张速查表分类函数一句话作用构造与转换ST_Point从经纬度构造点构造与转换ST_MakeLine把点集串成线构造与转换ST_GeomFromText从 WKT 文本构造几何构造与转换ST_SetSRID给几何设置坐标系编号构造与转换ST_Transform在不同坐标系间转换构造与转换ST_AsText把几何输出为 WKT 文本构造与转换ST_AsGeoJSON把几何输出为 GeoJSON构造与转换ST_Envelope获取几何外包矩形空间关系判断ST_Intersects判断两个几何是否相交空间关系判断ST_Contains判断是否包含空间关系判断ST_Within判断是否被包含空间关系判断ST_Touches判断是否仅边界接触空间关系判断ST_Crosses判断是否交叉穿过空间关系判断ST_Overlaps判断是否部分重叠空间关系判断ST_Covers判断是否覆盖比 Contains 更宽松空间关系判断ST_DWithin判断距离是否在某阈值内几何处理与聚合ST_Buffer生成缓冲区几何处理与聚合ST_Centroid求几何中心点几何处理与聚合ST_ConvexHull求凸包几何处理与聚合ST_Simplify简化几何顶点几何处理与聚合ST_Intersection求两个几何的交集几何处理与聚合ST_Difference求两个几何的差集几何处理与聚合ST_Union合并几何可聚合多行几何处理与聚合ST_Collect收集多行几何为一个集合测量计算ST_Distance求两个几何距离测量计算ST_Length求线的长度测量计算ST_Area求面的面积测量计算ST_DistanceSphere球面距离计算测量计算ST_ClosestPoint求两几何最近点测量计算ST_ShortestLine求两几何最短连线表里看着多其实真正要记的是每个函数“解决什么问题”。函数名本身已经很语义化了用一次比背十遍都管用。2.2 构造与转换类先让数据变成“空间数据”最基本的操作是从无到有造出几何对象或把已有坐标文本转成几何。其中最常用的是SELECT ST_Point(116.3, 39.9); SELECT ST_GeomFromText(POINT(116.3 39.9), 4326); SELECT ST_MakeLine(ARRAY[ ST_MakePoint(116.3, 39.9), ST_MakePoint(116.4, 40.0) ]::geometry[]);很多人容易忽略一个关键点ST_Point(x, y)的参数顺序是先经度后纬度也就是先 x 后 y。业务系统里如果一直用“纬度,经度”的字段排序到这儿就很容易把点写反画出来的位置完全不对。ST_SetSRID和ST_Transform两个函数经常被搞混。我打个比方ST_SetSRID只是给一个几何对象“贴标签”告诉数据库这些数字属于哪个坐标系并不会改变坐标数值ST_Transform才是真正的“坐标换算”会把一组坐标系下的数字转成另一组坐标系下的数字。-- 只设置 SRID不改变坐标值 SELECT ST_SetSRID(ST_MakePoint(116.3, 39.9), 4326); -- 从 WGS84 转 Web Mercator SELECT ST_Transform( ST_SetSRID(ST_MakePoint(116.3, 39.9), 4326), 3857 );输出格式也是日常非常常用的能力。尤其是前后端交互时GeoJSON 基本是事实标准SELECT ST_AsGeoJSON(ST_Transform(geom, 4326)) FROM buildings WHERE id 1;如果你只需要数据的“大概范围”ST_Envelope能快速拿到外包矩形很多地图缩放到某个要素集合的场景都会用到。2.3 空间关系判断类空间 SQL 的“WHERE 条件”这类函数本质上就是空间查询的过滤器。比如找出一个商圈 500 米范围内的所有餐厅或者判断某个小区是否在生态红线内。它们的返回值是布尔值直接用在WHERE里SELECT poi.name FROM poi, boundary WHERE ST_Intersects(poi.geom, boundary.geom);几个拓扑关系的关系很容易混淆我简单理一下ST_Contains(A, B)A 完全包裹 B并且 B 的任何部分不在 A 之外。ST_Within(A, B)是ST_Contains(B, A)的镜像。ST_Touches(A, B)两个几何只在边界上接触内部区域没有重叠。比如两个相邻地块共边就符合ST_Touches。ST_Crosses(A, B)两条线交叉穿过但交点没有内部重叠。常用于判断道路是否穿过一个区域。ST_Overlaps(A, B)两个几何有一部分内部重叠但又不完全互相包含。ST_Covers(A, B)跟ST_Contains类似但对边界处理更宽松B 在 A 的边界上也判断为 True。实际项目中我查“某个点是否落在面内”经常直接写ST_Covers因为它对边界情况的语义干扰最小。而ST_DWithin的实用性更强它不等价于ST_Distance N但配合索引时性能有巨大优势这个后续详细说。2.4 几何处理与聚合类从“单个对象”到“批量派生对象”做空间数据分析最爽的场景就是把一个复杂几何变成更有业务含义的几何对象。ST_Buffer是最典型的例子。比如给某个加油站点生成一个 3 公里服务半径面SELECT ST_Buffer(geom, 0.03) FROM stations WHERE id 10;注意如果几何是 43260.03 就代表 0.03 度不是 30 米。要做米制缓冲区最好先把几何转成投影坐标系或者使用geography。如果只需要知道“哪些点在多远以内”优先用ST_DWithin别动不动就ST_Buffer生成一个大面。ST_Centroid求几何质心常用来做点标注SELECT ST_Centroid(geom) FROM districts;ST_ConvexHull就是求一组点或几何的最小外凸多边形很多“覆盖范围分析”的粗糙版本会用它。ST_Simplify在数据瘦身时非常管用。路网数据几百万个顶点直接渲染只会拖垮前端。适当简化既能保持形状又能大幅降体积SELECT ST_Simplify(geom, 0.0001) FROM roads;容差参数取决于坐标系。如果坐标是度0.0001 大约就是几十米级别的简化幅度具体要根据比例尺反复试。ST_Intersection、ST_Difference、ST_Union和ST_Collect是几何“加减乘除”的四件套-- 两个面的交集 SELECT ST_Intersection(a.geom, b.geom) FROM parcels a, flood_zone b WHERE ST_Intersects(a.geom, b.geom); -- 面 A 减去面 B SELECT ST_Difference(a.geom, b.geom) FROM parcels a, red_line b; -- 将多行几何聚合为一个集合 SELECT ST_Union(geom) FROM counties; SELECT ST_Collect(geom) FROM points;ST_Union在这里有两种用法聚合函数写法会把多行合并成一个几何而二元函数写法会把两个几何合并。ST_Collect只做简单收集不做拓扑合并所以速度更快如果不需要消除内部边界优先考虑ST_Collect。2.5 测量计算类距离、长度、面积一网打尽测量类是日常查询里出镜率最高的。ST_Distance算距离ST_Length算线长ST_Area算面积但最容易出问题的点在于单位。直接对 4326 的geometry调用这些函数计算结果以“度”为单位。于是在中国地区一个“面积”算出来可能是零点几甚至更小的数字看起来完全不对。两条常规做法转投影坐标系后再算。比如全国范围用 3857 或者适合本地的 UTM 分区。直接使用geography类型让 PostGIS 自动按球面/椭球面计算。SELECT ST_Distance( ST_SetSRID(ST_MakePoint(116.3, 39.9), 4326)::geography, ST_SetSRID(ST_MakePoint(116.4, 39.9), 4326)::geography );ST_DistanceSphere和ST_ShortestLine这类函数更多用在“近似距离”“最近位置”的查询场景。比如判断两个多边形最近距离或者找到两个面之间的最短连线用于可视化或进一步分析。ST_Area同理SELECT ST_Area(geom::geography) AS area_m2 FROM land_parcels;这一节我不想每个函数都给一个例子因为后续“距离函数”章节会把ST_Distance、ST_DWithin、ST_DistanceSphere展开得更细。现在先把 30 个函数的地图建立起来真正用的时候你对号入座就行。3. 距离函数实战从“长度”到“业务权重”3.1 三种距离计算方式的取舍距离函数是空间数据库里最容易被误用的一个类别。同一个词“距离”在不同坐标系和类型下返回的单位、精度、性能都不同。第一种ST_Distance作用于geometry返回的是投影平面上的欧氏距离。如果几何是 4326返回单位是度如果是 3857返回单位是米。好处是计算快坏处是跨投影带时会有精度损失而且忘转投影就直接算的话结果基本不能用。第二种ST_Distance作用于geography。PostGIS 会按球面模型计算大圆距离返回单位是米。精度比平面计算高很多而且不用考虑投影选择。缺点是计算开销更大不适合在超大表上直接做全表两两距离。第三种ST_DistanceSphere。它绕过了椭球体精确计算直接用球体近似速度比geography快精度略低。对于打车、外卖这种对距离精度不敏感的场景完全够用。我平时给团队定的规范是数据入库统一存 4326geometry距离查询视精度要求决定是否 cast 成geography如果要做大量两两距离计算先建 GiST 索引再结合ST_DWithin做前置过滤。3.2 距离衰减函数的业务玩法很多人在 PostGIS 文档里找不到“distance decay function”这个词因为它不是一个内置函数名而是一种空间分析方法。简单说就是让“地理距离”转化为“业务权重”时距离越远权重越低衰减速度由模型决定。常见的三种衰减模型线性衰减weight 1 - d / D到达最大距离 D 时权重归零。指数衰减weight exp(-d / D)衰减速度先快后慢。高斯衰减weight exp(-0.5 * (d / D)^2)短距离内衰减缓慢过了阈值后快速下降。我在做商圈选址分析时常用这样的 SQLWITH dist AS ( SELECT poi.id, ST_Distance( shop.geom::geography, poi.geom::geography ) AS d FROM shop, poi WHERE ST_DWithin(shop.geom, poi.geom, 5000) ) SELECT id, CASE WHEN d 1000 THEN 1.0 WHEN d 5000 THEN 1.0 - (d - 1000) / 4000 * 0.5 ELSE 0.0 END AS linear_weight, EXP(-d / 2000.0) AS exp_weight FROM dist;注意ST_DWithin第一个参数是几何第二个也是几何阈值单位跟随坐标系。为了能在 5000 米的场景中使用我更习惯把几何 cast 成 geography 后比较但这时要确认索引匹配。最简单稳妥的做法是给几何建 GiST 索引然后查询里一开始就统一投影或统一成 geography。距离衰减函数本身不是 PostGIS 的计算瓶颈真正的瓶颈往往是你做了全表笛卡尔积式的距离计算。所以业务上一定要先圈定候选集再做权重计算。这也是ST_DWithin最大的价值所在。3.3 性能与索引为什么 ST_DWithin 比 ST_Distance 更讨喜很多新手写“5000 米内所有 POI”会写成SELECT * FROM poi WHERE ST_Distance(shop.geom, poi.geom) 5000;这个写法在数据量小的时候没毛病但在百万级 POI 表上会非常痛苦因为它无法有效利用 GiST 索引做边界框过滤。正确姿势是用ST_DWithinSELECT * FROM poi WHERE ST_DWithin(shop.geom, poi.geom, 0.05);配合 GiST 索引后PostGIS 会先做外包框相交判断快速筛掉绝大多数记录再精确计算少量候选。两者在业务效果上几乎等价但性能差距可能是几十倍。CREATE INDEX idx_poi_geom ON poi USING GIST(geom);这条规则记牢能用 ST_DWithin 判断邻近关系就别用 ST_Distance 比较阈值。4. pgRouting让路网真正“会算路”4.1 安装、建拓扑和成本模型pgRouting 是 PostGIS 的黄金搭档专门解决图论路径问题。扩展安装很简单CREATE EXTENSION IF NOT EXISTS pgrouting;难的从来不是安装而是把一份普通路网表变成可计算的拓扑图。pgRouting 的路径算法是基于“节点 边”的所以路网表里必须有source和target字段分别代表边的起点和终点顶点编号。假设路网表叫roads字段为id, geom接下来是标准流程ALTER TABLE roads ADD COLUMN source integer; ALTER TABLE roads ADD COLUMN target integer; ALTER TABLE roads ADD COLUMN cost double precision; ALTER TABLE roads ADD COLUMN reverse_cost double precision;接着用pgr_createTopology根据几何上的端点自动生成节点并填充 source/targetSELECT pgr_createTopology(roads, 0.0001, geom, id);这个 0.0001 是容差。如果路网坐标是度这个值大约对应十余米量级如果端点坐标因为数据采集问题有一点点缝隙小于容差的端点会被认为是同一个节点。容差设得太小拓扑会出现断头路设得太大又会把邻近但不相连的道路误连到一起。我在实际项目中一般先取数据精度的 2 到 3 倍再通过抽几条路线目视验证。最后填充成本UPDATE roads SET cost ST_Length(geom::geography), reverse_cost ST_Length(geom::geography);cost是正向行驶成本reverse_cost是反向行驶成本。如果允许双向通行且成本对称两个字段都设置同样的值。如果遇到单行道把逆向成本设为 -1 表示不可通过UPDATE roads SET reverse_cost -1 WHERE oneway true;pgRouting 约定负值表示该方向不连通这是新手最容易忽略的规则。4.2 pgr_dijkstra最短路径查询Dijkstra 算法是最经典的单源最短路径算法。pgRouting 中的标准调用方式如下SELECT * FROM pgr_dijkstra( SELECT id, source, target, cost, reverse_cost FROM roads, 1, 50, directed : false );参数说明第一个参数是边的 SQL第二个参数是起点顶点编号第三个参数是终点顶点编号directed决定是否按有向图处理。返回结果里比较重要的列是seq结果顺序编号。path_seq路径上边的顺序。node当前节点编号。edge当前边编号。如果 edge 为 -1表示已经是终点。cost当前边的行驶成本。只拿一条路径还不够业务上通常要连回原始路网几何把路径画出来SELECT r.seq, rd.geom FROM pgr_dijkstra( SELECT id, source, target, cost, reverse_cost FROM roads, 1, 50, directed : false ) r LEFT JOIN roads rd ON r.edge rd.id WHERE r.edge -1 ORDER BY r.seq;很多人在调用最短路径前会卡在一个问题上起点和终点是经纬度坐标但算法要的是节点编号。怎么找到最近的拓扑节点我用 KNN 查询SELECT id FROM roads_vertices_pgr ORDER BY geom - ST_SetSRID(ST_MakePoint(116.3, 39.9), 4326) LIMIT 1;-是 PostGIS 的 KNN 距离算子配合 GiST 索引可以快速找到最近节点。roads_vertices_pgr是pgr_createTopology自动创建的顶点表。求整条路径总成本也简单SELECT SUM(cost) FROM pgr_dijkstra( SELECT id, source, target, cost, reverse_cost FROM roads, 1, 50, directed : false );4.3 A* 与 KSP不同场景的差异化选型Dijkstra 虽然稳定但在大规模路网上计算较慢。pgr_aStar 在 Dijkstra 基础上加了启发式估算理论上可以减少搜索范围。它的 SQL 需要额外提供起点和终点的坐标字段SELECT * FROM pgr_aStar( SELECT id, source, target, cost, reverse_cost, x1, y1, x2, y2 FROM roads, 1, 50, directed : false );这里的 x1、y1、x2、y2 会由pgr_createTopology生成到顶点表中我们可以通过关联把坐标字段带进边数据里。实际项目中我通常默认用 Dijkstra只有路网特别大、对响应时间敏感时才切换 A*因为 A* 的启发式效果受数据质量影响很大网格状路网表现优秀但复杂路网可能并不比 Dijkstra 快多少。pgr_ksp则用于“前 K 条最短路径”适合需要给用户多个备选方案的场景SELECT * FROM pgr_ksp( SELECT id, source, target, cost, reverse_cost FROM roads, 1, 50, 3, directed : false );最后一个参数 3 表示返回前 3 条候选路径。返回结果里route_id可以区分不同备选路线。4.4 pgr_drivingDistance等时圈与可达范围最短路径解决的是“点到点”但很多业务要的是“从某点出发成本不超过某值能到哪些地方”。比如外卖配送的 30 分钟可达范围分析。这类问题用pgr_drivingDistance很合适SELECT * FROM pgr_drivingDistance( SELECT id, source, target, cost, reverse_cost FROM roads, 1, 30, directed : false );第三个参数 30 表示最大行驶成本。返回的是可达节点以及到达该节点的最小成本。把节点连回顶点表几何就能画出等时圈范围这比用圆面缓冲区更贴近真实路网约束。5. 高频踩坑与排查实录5.1 安装与扩展类问题现象原因解决办法CREATE EXTENSION 提示控制文件不存在PostGIS 包未安装或版本不匹配检查 PostgreSQL 版本与 PostGIS 版本对应关系重装扩展包执行 SQL 时提示 function st_distance(geometry, geometry) does not exist扩展建到了其他库在业务库执行 CREATE EXTENSION postgis;PostgreSQL 升级后扩展报错二进制不兼容重新安装匹配版本的 PostGIS并执行 ALTER EXTENSION ... UPDATE补充一个容易被忽略的问题如果数据库服务由云平台托管如 RDS 或托管 PostgreSQL部分云厂商默认不给普通用户创建扩展的权限需要到参数组或控制台开启postgis白名单。5.2 数据与坐标系类问题面积数值异常大或异常小。多半是坐标系问题。4326 直接算面积单位是度数值会异常3857 算面积不是真实地面面积在高纬度地区误差很大。最稳的办法是转geography或本地 UTM 投影。最短路径结果“绕远路”。先查拓扑是否完整source、target是否大量为空。如果pgr_createTopology容差设置偏小道路断开导致绕路是常态。用如下 SQL 快速检查悬空边数量SELECT count(*) FROM roads WHERE source IS NULL OR target IS NULL;cost 为度值导致路径结果失真。如果直接用ST_Length(geom)计算 cost而 geom 是 4326成本会以“度”为单位导致短距离路段和长距离路段在算法看来没有正确区分。务必用ST_Length(geom::geography)或者转换到投影坐标系后计算。5.3 性能与距离衰减调参实战做距离衰减分析时最常见的性能灾难是“两两全表扫描”。我见过一张 50 万行 POI 表和其他表做笛卡尔积距离计算一条 SQL 跑十几分钟。排查后发现WHERE ST_DWithin(...)里用了函数表达式包住索引列导致 GiST 索引完全失效。解决方式有两种一是保证参与索引的列是原始几何列不要在外面套ST_Transform或ST_SetSRID二是如果查询经常做投影转换直接增加一个生成列把转换后的几何存储下来并建索引。距离衰减模型中的阈值参数也需要根据业务反复调。比如指数衰减里的常数 D决定曲线在什么距离上衰减到约 36.8%。如果 D 设得太小稍微远一点的 POI 权重就趋近于零设得太大则距离对业务几乎没有区分度。我的经验是先从业务实际的核心服务半径出发比如外卖 3 公里、商场 5 公里把 D 设成半径的一半左右再用线上的转化率数据做 A/B 验证。我在实际项目中最后沉淀下来一套工作流程所有空间表先统一 SRID 并建 GiST 索引所有“邻近判断”优先 ST_DWithin所有距离成本计算统一转 geography 或投影坐标路径分析前先跑拓扑检查和少量可视化验证。这套流程虽然看起来平淡但能帮你在空间数据库这条路上少走很多弯路。我自己刚开始学空间 SQL 时也是先背函数名后来才发现真正值钱的不是记住 30 个函数而是知道在哪种业务场景下该用哪个函数、为什么用、单位是什么、索引怎么建。把这些底层逻辑搞清楚了学 pgRouting 时你会觉得顺理成章因为最短路和分析函数的组合应用本质还是空间数据思维的问题。