PostGIS 深度解析：从 OGC Simple Features 到 PG 18 时代的 GIS 引擎

本文是「PostgreSQL 扩展系列」GIS 篇。同系列前文：

• PostgreSQL 前世今生：从 1986 Berkeley 实验室到 2025 全球基础设施
• PostgreSQL 核心特性全景：12 大功能的设计动因、实现原理与版本演化
• PostgreSQL 元数据存储机制

PostGIS 不是简单的”PG 加几何函数”，而是 **PG 历史上第一个非贡献者、内核之外的”超级扩展”**——它有自己的数据类型、自己的索引、自己的函数库、自己的类型转换、自己的栅格子系统、自己的拓扑子系统、自己的点云子系统。

它 24 年还在持续迭代，已经在 ArcGIS、QGIS、Mapbox、OpenStreetMap、NASA、欧洲航天局、Uber 等无数 GIS / 位置服务的底层。

本文回答 4 个问题：

1. PostGIS 是怎么从 OGC Simple Features 变成 PG 的内置能力的？
2. geometry / geography / raster / topology / pointcloud 五大子系统各自干什么？
3. GiST 索引 + PROJ 投影 + GEOS 拓扑 + SFCGAL 三维的几何运算链路？
4. PostGIS 3.x → 4.x 的未来走向？

全文 8 大章节，25+ 张架构图，120+ 个 SQL/SQL-MM 示例。

一、PostGIS 在 PG 扩展生态的位置

1.1 PostGIS 是什么

1.2 PostGIS 的 4 个”非 OGC”能力

1.3 PostGIS vs Oracle Spatial vs SQL Server Geography

二、PostGIS 历史：从 2001 到 2025

2.1 24 年大事记

2.2 关键人物

三、PostGIS 架构：6 大子系统

3.1 子系统依赖图

依赖原则：geometry 是核心依赖，其他模块都依赖它。geography 用 geometry 实现球面计算（CAST AS geometry）。

四、PostGIS 数据类型：OGC Simple Features

4.1 OGC SFA 7 种几何

OGC SFA（Simple Features Access）定义了 7 种基础几何类型 + 14 种继承子类型：

4.2 WKT / WKB 编码

SELECT 'POINT(2 5)'::geometry;
SELECT 'LINESTRING(0 0, 1 1, 2 0)'::geometry;
SELECT 'POLYGON((0 0, 1 0, 1 1, 0 1, 0 0))'::geometry;
SELECT 'MULTIPOINT((0 0), (1 1))'::geometry;
SELECT 'GEOMETRYCOLLECTION(POINT(2 0), LINESTRING(0 0, 1 1))'::geometry;

WKB（Well-Known Binary）：

WKB Point:
- byte order: 01 (little endian)
- type: 01000000 (Point)
- x: 0000000000000040 (2.0)
- y: 0000000000001440 (5.0)

4.3 Geometry 内部存储

源码 src/postgis.h（简化）：

typedef struct {
    uint32_t type;        // geometry type
    uint32_t flags;       // hasZ | hasM | hasSRID | ...
    uint32_t srid;        // spatial reference ID
    uint8_t data[1];     // 实际坐标数据
} GSERIALIZED;

五、Geometry 子系统

5.1 Geometry 函数概览

5.2 OGC SFA SQL/MM 函数分类

5.3 实际示例

-- 创建空间表
CREATE TABLE cities (
    id SERIAL PRIMARY KEY,
    name TEXT NOT NULL,
    location GEOMETRY(Point, 4326) NOT NULL
);

-- 插入几何
INSERT INTO cities (name, location) VALUES
  ('Beijing', ST_SetSRID(ST_Point(116.4074, 39.9042), 4326)),
  ('Shanghai', ST_SetSRID(ST_Point(121.4737, 31.2304), 4326)),
  ('Tokyo', ST_SetSRID(ST_Point(139.6917, 35.6895), 4326));

-- 查询上海 1000km 内的城市
SELECT name, ST_Distance(location, ST_SetSRID(ST_Point(121.4737, 31.2304), 4326)) / 1000 AS km
FROM cities
WHERE ST_DWithin(
    location::geography,
    ST_SetSRID(ST_Point(121.4737, 31.2304), 4326)::geography,
    1000000  -- 1000km
);

六、GiST 索引：R-Tree + SP-GiST + BRIN

6.1 PostGIS 支持的 4 种索引

6.2 GiST R-Tree 原理

bbox (Bounding Box) 是几何的外接矩形（4 个 double），作为索引键。

6.3 索引使用

-- GiST 索引（默认）
CREATE INDEX idx_cities_location ON cities USING GIST (location);

-- SP-GiST 索引
CREATE INDEX idx_cities_location_spgist ON cities USING SPGIST (location);

-- BRIN 索引（适合时序）
CREATE INDEX idx_cities_location_brin ON cities USING BRIN (location);

-- 查询触发索引
EXPLAIN ANALYZE
SELECT name FROM cities
WHERE location && ST_MakeEnvelope(120, 30, 125, 35, 4326);
-- && 是 "bbox 共享" 操作符，触发 GiST / SP-GiST

6.4 GiST 接口

PostGIS 通过 PG 的 GiST extension API 实现：

// src/backend/access/gist/gistproc.c
Datum gist_geometry_consistent(PG_FUNCTION_ARGS) {
    GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0);
    StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2);
    // strategy = 1: <<(left), 2: &>(right), 3: && (overlap), ...
}

PG 提供给扩展的接口：

七、PROJ 投影：10000+ 坐标参考系统

7.1 什么是 SRID

SRID（Spatial Reference ID） 是 EPSG 数据库里的唯一 ID，例如：

7.2 PG 的 spatial_ref_sys 表

-- 安装 PostGIS 后
SELECT srid, auth_name, auth_srid, srtext FROM spatial_ref_sys WHERE srid = 4326;
-- 4326 | EPSG | 4326 | GEOGCS["WGS 84",DATUM["WGS_1984",...]

-- 列出 PG 支持的 10000+ SRID
SELECT count(*) FROM spatial_ref_sys;

7.3 投影转换

-- WGS 84 → Beijing-54 (高斯克吕格)
SELECT ST_AsText(
    ST_Transform(
        ST_SetSRID(ST_Point(116.4074, 39.9042), 4326),  -- WGS 84
        2415                                           -- Beijing-54
    )
);
-- 输出：POINT(5011999.9998 3938599.9999)

-- 中国境内常用：
-- 4326 → 4490 (CGCS2000)
-- 4326 → 4513 (CGCS2000 / 3-degree Gauss-Kruger CM 117E)

7.4 PROJ pipeline

PostGIS 3.0+ 使用 PROJ pipeline 表达投影链：

+proj=pipeline
  +step +proj=longlat +datum=WGS84
  +step +proj=tmerc +lat_0=0 +lon_0=117 +k=1 +x_0=500000 +y_0=0 +datum=WGS84

八、Geography：球面几何

8.1 Geometry vs Geography

8.2 何时用 Geography

8.3 内部实现

源码 src/geography.c：

Datum geography_distance(PG_FUNCTION_ARGS) {
    GSERIALIZED *g1 = PG_GETARG_GSERIALIZED_P(0);
    GSERIALIZED *g2 = PG_GETARG_GSERIALIZED_P(1);
    
    // 1. cast 到 geometry
    Geometry *ga = lwgeom_from_gserialized(g1);
    Geometry *gb = lwgeom_from_gserialized(g2);
    
    // 2. 转球面 (geographic lib / spheroid)
    POINT3D A = ..., B = ...;
    
    // 3. 大圆距离 (Haversine / Vincenty)
    return geography_spheroid_distance(A, B);
}

九、GEOS 几何运算

9.1 GEOS 是什么

9.2 GEOS 操作分类

9.3 索引触发 GEOS

-- ST_Intersects 自动用索引
CREATE INDEX idx_parcels_geom ON parcels USING GIST (geom);
SELECT * FROM parcels
WHERE ST_Intersects(geom, ST_MakeEnvelope(...));
-- 流程：先 bbox 索引过滤，再用 GEOS 精确判断

-- ST_DWithin 自动用索引
SELECT * FROM cities
WHERE ST_DWithin(location, target::geography, 1000);

9.4 R-Tree 索引加速示例

-- 准备 100 万条路网数据
CREATE TABLE roads (id SERIAL, geom GEOMETRY(LineString, 4326));
INSERT INTO roads SELECT ST_GeneratePoints(...) FROM generate_series(1, 1000000);

-- 不加索引
EXPLAIN ANALYZE
SELECT count(*) FROM roads
WHERE ST_DWithin(geom::geography, ST_Point(...), 5000);
-- 30+ 秒

-- 加索引
CREATE INDEX idx_roads_geom ON roads USING GIST (geom);
-- 同一查询：0.05 秒（1000x 加速）

十、Raster 子系统

10.1 Raster 类型

10.2 Raster vs Geometry 集成

-- 加载 GeoTIFF
CREATE TABLE dem (
    rid SERIAL,
    rast RASTER
);
INSERT INTO dem (rast) VALUES
    (ST_FromGDALRaster('/path/to/dem.tif'));

-- 高程查询
SELECT ST_Value(rast, ST_SetSRID(ST_Point(...), 4326))
FROM dem;

-- Raster → Polygon (轮廓提取)
SELECT (ST_DumpAsPolygons(rast)).geom, (ST_DumpAsPolygons(rast)).val
FROM dem
WHERE ST_Intersects(rast, ST_MakeEnvelope(...));

10.3 3.5 新增 ST_AsMVT

-- 生成 Mapbox Vector Tile (MVT)
SELECT ST_AsMVT(q, 'cities', 4096, 'geom')
FROM (
    SELECT name, location AS geom
    FROM cities
) AS q;

十一、Topology 子系统

11.1 Topology 概念

11.2 Topology 优势 vs Geometry

11.3 Topology 用例

-- 创建一个 topology
SELECT topology.CreateTopology('my_topo', 4326);

-- 加一个面（自动维护拓扑）
INSERT INTO my_topo.face (topo_id, mbr)
SELECT topology.createFaceFromPolygon('MH', 4326, ...);

-- 查所有面
SELECT * FROM my_topo.face;

-- 移动边（拓扑自动更新关联 face）
SELECT topology.ST_MoveIsoNode('my_topo', ..., ...);

十二、PostGIS 性能优化

12.1 10 条优化建议

12.2 详细例子

-- 1. 总是加 GiST 索引
CREATE INDEX idx_parcels_geom ON parcels USING GIST (geom);

-- 2. simplify (Douglas-Peucker)
UPDATE parcels SET geom = ST_Simplify(geom, 0.5);
-- 50% 大小缩减，几何相似度 > 95%

-- 5. 复合索引
CREATE INDEX idx_parcels_district_geom ON parcels USING GIST (geom) INCLUDE (district_id);

-- 6. ST_SnapToGrid
UPDATE addresses SET geom = ST_SnapToGrid(geom, 0.0001);
-- 减少 7 位小数精度

-- 7. Cluster
CLUSTER parcels USING idx_parcels_geom;

12.3 常见 5 个陷阱

十三、PostGIS 3.x → 4.x

13.1 4.0 路线

13.2 4.0 计划特性

1. PointCloud 重写 — 基于 PLY 4.0 协议
2. Parallel ST_AsMVT — 并行生成 MVT tile
3. 3D GeoParquet — Geoparquet 1.1 IO
4. BRIN 压缩改进 — 更小 BRIN 索引
5. PG 18 + columnar 试验集成

十四、PostGIS 实战 5 大场景

14.1 位置查询（POI 找附近）

-- "北京 5km 内所有咖啡店"
SELECT name, address
FROM cafes
WHERE ST_DWithin(
    location::geography,
    ST_SetSRID(ST_Point(116.4074, 39.9042), 4326)::geography,
    5000
)
ORDER BY ST_Distance(
    location::geography,
    ST_SetSRID(ST_Point(116.4074, 39.9042), 4326)::geography
)
LIMIT 20;

14.2 多边形包含（区域查询）

-- "包含在某区域的所有路"
SELECT *
FROM roads
WHERE ST_Intersects(geom, (SELECT geom FROM regions WHERE name='Beijing'));

14.3 GeoJSON 输出

-- 输出 GeoJSON FeatureCollection
SELECT json_build_object(
    'type', 'FeatureCollection',
    'features', json_agg(json_build_object(
        'type', 'Feature',
        'geometry', ST_AsGeoJSON(geom),
        'properties', json_build_object('id', id, 'address', address))
    ))
FROM cafes
WHERE ST_DWithin(
    location::geography,
    ST_SetSRID(ST_Point(116.4074, 39.9042), 4326)::geography,
    5000
);

14.4 MVT 地图瓦片

-- 生成矢量瓦片（XYZ schema）
SELECT ST_AsMVT(q, 'cities', 4096, 'geom')
FROM (
    SELECT
        name,
        ST_AsMVTGeom(
            location,
            ST_TileEnvelope(z, x, y),
            4096, 64, true
        ) AS geom
    FROM cities
    WHERE location && ST_TileEnvelope(z, x, y)
) AS q;

14.5 路由分析（pgRouting 集成）

-- pgRouting 扩展
CREATE EXTENSION pgrouting;

-- 最短路径
SELECT * FROM pgr_dijkstra(
    'SELECT id, source, target, cost, reverse_cost FROM ways',
    -- 起点 / 终点
    (SELECT id FROM ways_vertices WHERE the_geom = ST_SetSRID(ST_Point(...), 4326)),
    (SELECT id FROM ways_vertices WHERE the_geom = ST_SetSRID(ST_Point(...), 4326)),
    directed := true
);

十五、PostGIS 5 大生产案例

十六、PostGIS 设计哲学：6 个原则

十七、PostGIS 学习路径

源码引用索引

• src/postgis.h — Geometry 数据结构
• src/geometry.c — Geometry 操作
• src/lwgeom.c — lwgeom 内部
• src/geography.c — Geography 实现
• src/raster/rt_api.c — Raster API
• src/topology/ — Topology 子系统
• src/sfcgal/ — SFCGAL 三维
• externals/proj/ — PROJ C 库
• externals/geos/ — GEOS C 库

同系列前文

• PostgreSQL 前世今生：39 年演化史
• PostgreSQL 核心特性全景：12 大功能
• PostgreSQL 元数据存储机制