用 DuckDB 直接 SQL 查询 GeoJSON：不装 PostGIS 的轻量方案
--- title: 用 DuckDB 直接 SQL 查询 GeoJSON：不装 PostGIS 的轻量方案 ---
处理 GeoJSON 时有个尴尬的中间地带：文件小，jq 或 turf.js 就够了；数据量大、查询复杂，标准答案是 PostGIS。可为了查几次数据就搭一套 PostgreSQL 服务，实在不划算。DuckDB 的 spatial 扩展正好卡在这个空档上。它是个单文件的嵌入式分析数据库，没有服务进程，却能用标准 SQL 直接查 GeoJSON 文件。几百 MB 的全国区县边界文件，普通笔记本上跑起来也不费劲。
一行命令完成安装
DuckDB 本体是个几十 MB 的单文件程序，命令行版从官网下载就能跑，Python 用户执行 pip install duckdb。spatial 是官方扩展，首次使用装一次，之后每个会话加载即可：
INSTALL spatial;
LOAD spatial;spatial 扩展内置了 GDAL，所以它不只认 GeoJSON，Shapefile、GeoPackage、FlatGeobuf 这些常见格式都能直接读。
读文件就是一个函数
ST_Read 把 GeoJSON 文件当成一张表：properties 里的每个属性自动摊平成一列，几何体放在名为 geom 的列里。
SELECT COUNT(*) FROM ST_Read('quxian.geojson');
DESCRIBE SELECT * FROM ST_Read('quxian.geojson');第二条 DESCRIBE 用来看表结构，确认 name、adcode 这些属性列的类型。全程不建表、不导入，文件本身就是数据源。
三个高频查询
先说点面归属查询：给定一个经纬度，找出它落在哪个区县。这是行政区划数据最常见的用法。
SELECT name, adcode
FROM ST_Read('quxian.geojson')
WHERE ST_Contains(geom, ST_Point(116.397428, 39.90923));注意 ST_Point 的参数顺序是经度在前、纬度在后，与 GeoJSON 规范一致。这个顺序问题在规范四个坑那篇里详细讨论过，写反了不报错，只会默默查不到结果。
第二个是按行政区划代码筛选并导出子集。adcode 的前两位是省级代码，例如 44 代表广东省。从全国文件里切出广东所有区县，直接导出为新的 GeoJSON：
COPY (
SELECT * FROM ST_Read('quxian.geojson')
WHERE adcode LIKE '44%'
) TO 'guangdong.geojson'
WITH (FORMAT GDAL, DRIVER 'GeoJSON');COPY ... WITH (FORMAT GDAL) 是 spatial 扩展提供的导出通道，换个 DRIVER 参数就能导出 Shapefile 或 GeoPackage。
第三个：找出最"重"的几何体。用 ST_NPoints 统计每个面的顶点数，定位文件体积的主要贡献者。
SELECT name, ST_NPoints(geom) AS pts
FROM ST_Read('quxian.geojson')
ORDER BY pts DESC
LIMIT 10;排在前面的通常是海岸线复杂的沿海区县。确认了顶点大户之后，可以用Douglas-Peucker 抽稀或TopoJSON 压缩进行减重处理。
与 PostGIS、turf.js 如何选择
| 方案 | 安装成本 | 适合数据量 | 典型场景 |
|---|---|---|---|
| DuckDB spatial | 单文件程序或 pip 一行 | 百万级要素 | 本地分析、一次性 ETL |
| PostGIS | 需部署 PostgreSQL 服务 | 千万级以上、多人并发 | 生产数据库、长期服务 |
| turf.js | npm 引入 | 浏览器内存放得下的量 | 前端交互计算 |
| ogr2ogr | 需安装 GDAL | 任意（流式处理） | 格式转换批处理 |
需要长期运行的空间服务，还是选 PostGIS，配合 WKT 交换格式与其他系统对接。只是本地做一次性的筛选、统计、转换，DuckDB 的启动成本低得多。
两个容易踩的坑
第一个坑：ST_Area 在经纬度坐标下算出来的是"平方度"。直接对 4326 坐标的几何体调用 ST_Area，结果单位是度的平方，没有实际意义。想得到米制面积需要先投影：
SELECT name,
ST_Area(ST_Transform(geom, 'EPSG:4326', 'EPSG:3857', always_xy := true)) / 1e6 AS approx_km2
FROM ST_Read('quxian.geojson')
LIMIT 5;但要清楚 Web 墨卡托（EPSG:3857）不是等积投影，北纬 40° 处面积会被放大约 1.7 倍，这个数值只能用于同纬度地区的相对比较。需要精确面积时，应投影到 Albers 等积投影后再计算。
第二个坑：always_xy 参数不能省。GDAL/PROJ 按权威定义处理 EPSG:4326 时轴序是纬度在前，与 GeoJSON 的经度在前相反。ST_Transform 加上 always_xy := true 才能保持经度在前的顺序，这与 pyproj 里 always_xy=True 是同一个坑的两种写法。漏掉这个参数，转换结果的 X、Y 会整体对调，几何体直接飞出中国范围。
全国省、市、区县三级的行政区划边界 GeoJSON 数据，可以在 GeoJSONcn 主站按行政区划检索预览，下载后即可直接用上面的 SQL 进行查询分析。