postgis-spatial-sql

PostGIS 与空间 SQL

用途：编写正确且快速的空间 SQL。常见的两类问题是语义问题（混淆 geometry 与 geography、SRID 不匹配 → 结果错误）和性能问题（未使用索引 → 连接耗时数小时）；此技能可防范这两类问题。

何时数据库才是合适的工具

出现以下任一情况时，请从文件/GeoPandas 转用 PostGIS：要素数量超过几百万、需要并发读写、需要反复进行临时查询、需要在其上提供服务 API，或需要事务完整性。对于 GeoParquet 的一次性分析扫描，DuckDB Spatial 通常是最快的免安装方案——采用同样的 SQL 思路，无需服务端。

需求不完整时，不要把这一经验法则当作最终建议。应先了解当前及预测的数据量、并发情况、交付和变更模式、延迟/SLA、服务需求和运维责任（包括备份与恢复）。为两个可行的后端定义有代表性的摄取、连接和读取查询；只有在行数、连接基数、SRID、几何有效性和样本输出均一致后，才比较运行时间和资源使用情况。在当前回复中包含此基准测试和正确性验证计划；不要只提出稍后再起草。

模式基础

此可运行示例假设数据位于 UTM 33N 分区。请将 EPSG:32633 替换为针对实际关注区域验证过的投影 CRS。

CREATE TABLE parcels (
  id          bigint GENERATED ALWAYS AS IDENTITY PRIMARY KEY,
  parcel_no   text NOT NULL,
  landuse     text,
  area_m2     double precision,          -- unit in the name, always
  geom        geometry(MultiPolygon, 32633) NOT NULL
);
CREATE INDEX parcels_geom_gix ON parcels USING gist (geom);
ANALYZE parcels;

• 完整指定 geometry 列的类型：geometry(MultiPolygon, SRID)——未指定类型的 geometry 列会接受各种混杂数据。

• 加载时转换为 Multi* 类型（ST_Multi），避免 Polygon/MultiPolygon 混用造成问题。

• geometry 与 geography：区域分析时，使用投影 SRID 下的 geometry（速度快、函数集完整）；范围覆盖全球/多个分区且希望直接使用米作为单位、不想选择投影时，使用 geography（SRID 4326）（速度较慢、函数集较少）。切勿将数据存储为 4326 geometry，却调用 ST_Area 并期望得到平方米——结果单位是平方度。

• 切勿使用 EPSG:3857/Web Mercator 测量面积或长度。如果尚未确定分析 CRS，可使用 4326 geography 获得大地测量结果，或暂停并选择经过验证的本地/等积 CRS；不要将已知存在失真的 CRS 作为可运行的测量方案。

• 建议使用的任何已存储 geometry 列都必须指定类型和 SRID。 如果建议为重复测量添加“第二个投影 geometry 列”，就必须将其写成 geometry(<Type>, <SRID>)，并提供索引和用于填充该列的 ST_Transform，否则建议并不完整。将未指定类型的列作为解决方案，反而会重新引入原本要解决的混合 SRID 问题：

  ALTER TABLE parcels ADD COLUMN geom_32633 geometry(MultiPolygon, 32633);
  UPDATE parcels SET geom_32633 = ST_Transform(geom, 32633);
  CREATE INDEX parcels_geom_32633_gix ON parcels USING gist (geom_32633);
  

• 为每个 geometry 列创建 GiST 索引，并在批量加载后运行 ANALYZE；BRIN 仅适用于超大型、按空间排序且只追加数据的表。

• 加载方式：ogr2ogr -f PostgreSQL、shp2pgsql 或 GeoPandas to_postgis（小型/中型数据）。COPY 的速度比 INSERT 快几个数量级。

正确的空间谓词

• 用 ST_Intersects 表示“有任何接触”，用 ST_Contains/ST_Within 表示包含关系，用 ST_DWithin(a, b, dist) 表示邻近关系——切勿使用 ST_Distance(a,b) < dist（这种形式无法使用索引）。
• 经典的点在多边形连接：

SELECT p.id, a.district
FROM points p
JOIN admin a ON ST_Intersects(a.geom, p.geom);   -- GiST on both sides

• 使用距离运算符进行 KNN 最近邻查询（可借助索引）：

SELECT h.id, h.name
FROM hospitals h
ORDER BY h.geom <-> (SELECT geom FROM incident WHERE id = 42)
LIMIT 3;

现代 PostGIS 中，<-> 对 geometry 会按真实距离排序；表非常大时，可用 ST_DWithin 限定搜索范围。

性能实战指南

1. 先运行 EXPLAIN (ANALYZE, BUFFERS)——确认使用了 GiST 索引 （查找 "Index Scan ... _gix"）；大型空间连接中出现 Seq Scan 意味着需要重写查询，而不是换一台更大的服务器。
2. 每个谓词两侧的 SRID 必须相同——在连接谓词中使用 ST_Transform 会导致索引失效；应改为存储转换后的索引副本。
3. 大型多边形问题：国家／流域大小的几何体会让索引边界框失去作用 → 使用 ST_Subdivide 拆分到工作表中（与这类几何体连接时通常可提速 10-100×）。

以下示例假设存在 countries(country_id, geom)。

CREATE TABLE country_parts AS
SELECT c.country_id, part.geom
FROM countries AS c
CROSS JOIN LATERAL ST_Subdivide(c.geom, 256) AS part(geom);

CREATE INDEX country_parts_geom_gix ON country_parts USING gist (geom);
ANALYZE country_parts;

ST_Subdivide 是一个返回集合的函数；不要通过 (ST_Subdivide(...)).geom 访问其结果。

4. 在数据库中验证有效性：使用 ST_IsValid 审计，使用 ST_MakeValid 修复；如果写入方不可信，则添加 CHECK (ST_IsValid(geom))。
5. 为服务场景简化，而非用于分析：保留全分辨率几何体； 为 Web 层生成 ST_SimplifyPreserveTopology 副本或矢量瓦片 （ST_AsMVT）。
6. 在事务中批量更新；大量数据变更后运行 VACUUM ANALYZE。

常见分析模式

-- Area-weighted aggregation (e.g., population into custom zones)
SELECT z.zone_id,
       SUM(b.pop * ST_Area(ST_Intersection(z.geom, b.geom)) / ST_Area(b.geom)) AS pop_est
FROM zones z JOIN blocks b ON ST_Intersects(z.geom, b.geom)
GROUP BY z.zone_id;

-- Dissolve with attribute
SELECT landuse, ST_Multi(ST_Union(geom))::geometry(MultiPolygon, 32633) AS geom
FROM parcels GROUP BY landuse;

面积加权插值假设源单元内的密度均匀——报告结果时应说明这一假设。有效性修复应使用 ST_MakeValid， 绝不要使用 ST_Buffer(geom, 0)。

DuckDB Spatial 快速路径

INSTALL spatial; LOAD spatial;
SELECT a.name, count(*)
FROM 'admin.parquet' a, 'points.parquet' p
WHERE ST_Intersects(a.geom, p.geom)
GROUP BY a.name;

可直接读取 GeoParquet／Shapefile／GPKG，默认并行——非常适合无需服务器的一次性大型连接和管道步骤。没有 GiST；它会自行规划连接——请进行基准测试，不要想当然。

验证流程

1. 每个连接／叠加 CTE 后都运行行数核算查询。
2. 对所有涉及的表运行 SELECT DISTINCT ST_SRID(geom), GeometryType(geom)——一个查询即可排除两类经典错误。
3. 将 5 个输出要素叠加到底图上渲染查看（QGIS 可直接连接）——即使数值正确，几何体也可能毫无意义。
4. 将所有作为可运行代码呈现的 sql 代码块视为语法和别名边界：在给定的架构假设下，必须能从上到下执行。绝不要在其中放入尖括号占位符、省略号、伪代码、弃用的连接或不完整的别名。如果 SRID 等架构值未知，请询问，或将模板放在标注为 text 的代码块中。

常见问题检查清单

• 对 4326 几何体使用 ST_Area／ST_Length（单位为平方度）。
• 使用 EPSG:3857／Web Mercator 测量面积或长度（存在系统性失真）。
• 使用 ST_Distance < x 而不是 ST_DWithin（无法使用索引）。
• 在连接谓词中使用 ST_Transform。
• 几何列未指定类型且包含混合 SRID。
• 连接国家大小的多边形时未使用 ST_Subdivide。
• 用 buffer(0) 修复有效性（会悄然丢失部分内容）——应使用 ST_MakeValid。
• 向 Web 客户端提供全分辨率几何体。

执行约定

• **工作流：**检查架构、SRID、几何类型、数据量和查询目标；选择谓词和索引；编写便于审计的 CTE；检查执行计划；核对结果；安全地投入运行。
• **决策规则：**并发、重复或事务性空间工作负载使用 PostGIS；对于范围明确的一次性转换，如果服务器没有额外价值，则使用文件管道或 DuckDB Spatial。
• **验证流程：**断言 SRID 和几何不变量，核算每次连接的行数，比较使用索引的执行计划及耗时，在地图上抽样检查几何体，并测试边界语义。
• **失败模式：**遇到 SRID 混合、意外的多对多数据膨胀、无效几何体、无法使用索引的谓词、geography／geometry 单位混淆或无法解释的执行计划退化时，阻止发布。
• **交付物：**自包含的参数化 SQL 或迁移脚本，包含一致的 CTE／表别名、索引及其依据、查询计划证据、行数核算、抽样验证、预期架构、性能说明和回滚指导。
• **来源时效性：**在选择函数或执行计划前，请查阅权威来源注册表，了解已部署的数据库和扩展版本。