在基于 GIS 的分析工作里，栅格数据（比如 DEM、土地利用、气候栅格）越来越常见。 把这些数据放进 PostGIS 之后，就可以用 SQL 做各种空间分析，还能和矢量数据无缝联动。

这篇文章以一个 DEM GeoTIFF 为例，介绍以下几点：

  1. 如何用 raster2pgsql 导入栅格；
  2. 如何用 ST_Value 在指定点上取栅格值；
  3. 如何把处理后的栅格导出成 GeoTIFF / PNG 等文件；
  4. 如何在 SQL 里做简单“栅格计算”（Map Algebra）。

中途会穿插官方文档和资料链接，方便继续深挖。

为什么要把栅格放进 PostGIS？

传统工作流里，栅格往往是单独放在文件系统里，由 QGIS、GDAL 等工具读取。PostGIS 的好处是：

  • 统一管理：矢量 + 栅格都在一个数据库里，方便做联动分析；
  • SQL 查询：用 ST_Value、ST_MapAlgebra 等函数，直接在数据库中完成复杂空间计算；
  • 服务化友好：后端只需要写 SQL，就能给前端提供切片、统计、分析接口。

用 raster2pgsql 导入栅格数据

PostGIS 提供了工具 raster2pgsql，把 GDAL 支持的栅格格式（GeoTIFF 等）转换为 SQL，再由 psql 导入数据库。

示例命令：

# 1. 把 GeoTIFF 转成 SQL
raster2pgsql -s 3857 -I -C -M GTI3857.tif -F -t 1000x1000 public.gti > elev.sql

# 2. 用 psql 导入数据库
psql postgres://postgres:ege%2523geoshare.cn@192.168.50.200:5432/gis -f elev.sql

逐个拆解一下几个关键参数（结合官方文档中的说明）：

  • -s 3857：指定输出栅格的 SRID，这里是 Web Mercator（EPSG:3857）。如果文件本身有投影信息但没有对应 SRID，可以在这里强制指定。
  • -I：在栅格列上创建 GiST 空间索引，后面做 ST_Intersects 之类查询会快很多。
  • -C：应用一批标准栅格约束（SRID、像元大小、对齐等），保证 raster_columns 视图能正确识别这张表。
  • -M：导入完成后对栅格表执行 VACUUM ANALYZE，让查询计划更合理。
  • -F：增加一个 filename 字段，记录每个 tile 对应的源文件名。
  • -t 1000x1000：把栅格按 1000×1000 像元切成瓦片，每瓦一行。这样一张大图就拆成很多 tile，读写更快。
  • public.gti：目标 schema 和表名。
  • > elev.sql：把生成的 SQL 保存到 elev.sql 里，再由 psql 导入。

官方手册中也给出了类似的示例，用 100×100 的 tile 导入多个 GeoTIFF：

raster2pgsql -s 4326 -I -C -M -F -t 100x100 *.tif public.demelevation > elev.sql
psql -d gisdb -f elev.sql

实战建议

  • DEM 等大栅格，通常都要切 tile，-t 256x256 / 512x512 / 1000x1000 都是常见配置；
  • 记得加 -I -C -M，这是比较「省心」的一套组合；
  • 如果想保持源文件名信息（比如便于追溯来源），务必加 -F。

在 PostGIS 中查询栅格值：ST_Value

栅格导入之后，数据表结构大致类似下面这种样子：

\d public.gti
-- rid | rast | filename | ...

我们想在某个坐标点（例如投影为 EPSG:3857 的坐标 11397776, 4456894）上取 DEM 高程值，并乘上一些系数，就可以用以下 SQL：

SELECT 302.061 * 0.855 *
       ST_Value(rast, ST_SetSRID(ST_MakePoint(11397776, 4456894), 3857)) AS val
FROM gti
WHERE ST_Intersects(rast, ST_SetSRID(ST_MakePoint(11397776, 4456894), 3857));

这里的核心函数是 ST_Value(raster, geometry)：

  • 根据给定的几何（点、线、面），返回该位置上栅格像元的值；([PostGIS][3])
  • 省略波段时，默认使用第 1 波段。

WHERE ST_Intersects(rast, geom) 则是用来过滤掉跟这个点毫无关系的 tile，减少不必要的扫描。

也可以把几何换成线、面，PostGIS 会根据几何和栅格重叠区域返回对应像元的值；比如沿线路径提取栅格值等。

也可以用诸如 302.061 * 0.855 * ST_Value(...) 这样的数学运算，可以理解为：

  • ST_Value(...)：拿到原始高程值（假设单位是米）；
  • 302.061 * 0.855：应用某种业务上的系数（比如高度换算、功率因子等）；
  • AS val：输出为最后的业务指标。

一旦把 DEM 放进 PostGIS，就可以轻松做到以下这些事：

  • 按 坐标点列表 批量取值；
  • 按 矢量图层（行政区、风机点位、线路）做交集统计；
  • 按 条件过滤（比如高程 > 2000m 的栅格区域）筛选区域。

导出栅格：从 PostGIS 回到 GeoTIFF / PNG

很多时候，你在 PostGIS 中做完一系列处理，也许会希望把结果再导出为 GeoTIFF 或 PNG，给 QGIS、GeoServer 或其他系统使用。

使用 GDAL gdal_translate 导出

PostGIS 官方 FAQ 推荐使用 GDAL 的 gdal_translate 工具从数据库里导出栅格：

PGHOST=localhost \
PGPORT=5432 \
PGUSER=postgres \
PGPASSWORD=your_password \
gdal_translate \
    -of GTiff \
    "PG:dbname=gis table=gti" \
    output.tif
  • PG:... 是 GDAL 的 PostgreSQL 连接字符串，可以指定数据库、表名甚至 WHERE 条件；
  • -of GTiff 指定导出格式，这里是 GeoTIFF，当然也可以换成 PNG、JPEG 等。

如果你只想导出某部分数据，可以在连接字符串里加入 where= 过滤，例如只导出某个文件名对应的 tile：

gdal_translate \
    -of PNG \
    "PG:dbname=gis table=gti where='filename=''GTI3857.tif''' " \
    gti3857.png

也可以用 ST_Intersects 做空间过滤，只导出特定范围：

gdal_translate \
    -of PNG \
    "PG:dbname=gis table=gti \
       where='ST_Intersects(rast, ST_SetSRID(ST_Point(11397776,4456894),3857))'" \
    clipped.png

使用 ST_AsTIFF / ST_AsPNG / ST_AsGDALRaster 导出

如果你在 应用里直接用 SQL 输出二进制栅格（比如 Web API 返回文件流），可以用 PostGIS 的导出函数：

  • ST_AsTIFF(rast, ...)：返回 TIFF（二进制 bytea）；
  • ST_AsPNG(rast, ...)：返回 PNG；
  • ST_AsJPEG(rast, ...)：返回 JPEG；
  • ST_AsGDALRaster(rast, 'GTiff')：更通用的导出函数，支持 GDAL 支持的多种格式。

一个简单示例（导出一张合成后的 DEM）：

SELECT ST_AsTIFF(rast, 'LZW') AS tiff_bin
FROM gti
WHERE rid = 1;

在后端应用中，把 tiff_bin 当作文件内容写出即可。

PostGIS 手册里还有用 ST_AsGDALRaster + PostgreSQL 大对象接口，直接在服务器上写文件的做法，适合批量导出，但需要更高权限。

在 PostGIS 中处理栅格：简单的 Map Algebra

比如需要对单个栅格值做系数相乘，如果要对整张栅格做类似处理，可以用 ST_MapAlgebra 这一类“地图代数（Map Algebra）”函数。

假设我们想把 DEM 的值整体乘以 0.855（例如做某种修正），可以写成：

CREATE TABLE gti_scaled AS
SELECT
  rid,
  ST_MapAlgebra(
    rast,
    1,                            -- 第 1 波段
    '32BF',                       -- 输出像元类型（float）
    '[rast] * 0.855',             -- 表达式：对每一个像元值做运算
    NULL
  ) AS rast
FROM gti;

这里 ST_MapAlgebra(rast, pixeltype, expression, ...) 的作用是，对每个像元执行一个表达式运算（这里就是 [rast] * 0.855），生成一张新的栅格。

常见用法包括：

  • 把 DEM 从米转换成其他单位；
  • 对风速、辐射量类栅格按系数缩放或叠加；
  • 按条件做重分类（例如：高程>2000 记为 1，否则为 0）。

做完这些运算之后，你就可以，用 gdal_translate 再导出成新的 GeoTIFF，或者用 ST_AsPNG、ST_AsTIFF 在 Web 接口中直接输出计算结果图层。

典型工作流小结

一个围绕 PostGIS 栅格的完整流程通常是：

  1. 导入栅格

    • 使用 raster2pgsql + psql 导入 GeoTIFF / 其他栅格；
    • 选择合适的 tile 大小、加上索引与约束。

  2. 查询与分析

    • 用 ST_Value(rast, geom) 在点、线、面上提取栅格值；
    • 用 ST_Intersects 等空间过滤函数提高查询效率；([Crunchy Data][10])
    • 利用矢量图层（行政区、项目范围、线路）与栅格做叠加分析。

  3. 栅格处理

    • 用 ST_MapAlgebra 做加减乘除、重分类等 Map Algebra 运算；
    • 结合 ST_Union、ST_Clip 等函数做裁剪、合并。

  4. 导出结果

    • 用 gdal_translate 从 PostGIS 表导出为 GeoTIFF / PNG 等；
    • 或在 Web 服务中使用 ST_AsTIFF / ST_AsPNG / ST_AsGDALRaster 直接返回图像二进制。

参考资料