第12章-栅格数据处理
第12章：栅格数据处理
12.1 PostGIS Raster 概述
PostGIS Raster 是 PostGIS 的栅格数据扩展，提供了在 PostgreSQL 数据库中存储和分析栅格数据的能力。
12.1.1 栅格数据特点
┌─────────────────────────────────────────────────────────────┐
│ 栅格数据特点 │
├─────────────────────────────────────────────────────────────┤
│ 组成元素 │
│ ├── 像素（Pixel）：最小数据单元 │
│ ├── 波段（Band）：存储不同类型的数据值 │
│ └── 元数据：空间参考、范围、分辨率等 │
│ │
│ 常见格式 │
│ ├── GeoTIFF (.tif) │
│ ├── JPEG2000 (.jp2) │
│ ├── NetCDF (.nc) │
│ └── HDF (.hdf) │
│ │
│ 应用领域 │
│ ├── 卫星影像、航拍影像 │
│ ├── 数字高程模型（DEM） │
│ ├── 气象数据、环境数据 │
│ └── 土地覆盖、土地利用 │
└─────────────────────────────────────────────────────────────┘
12.1.2 启用 PostGIS Raster
-- 创建 raster 扩展
CREATE EXTENSION postgis_raster;
-- 查看版本
SELECT PostGIS_Raster_Version();
-- 查看支持的栅格格式
SELECT * FROM ST_GDALDrivers();
-- 启用外部栅格支持
ALTER SYSTEM SET postgis.enable_outdb_rasters = true;
ALTER SYSTEM SET postgis.gdal_enabled_drivers = 'ENABLE_ALL';
SELECT pg_reload_conf();
12.2 栅格数据导入
12.2.1 使用 raster2pgsql
# 基本导入
raster2pgsql -s 4326 -I -C -M dem.tif public.dem | psql -d gis_db
# 参数说明
# -s 4326 : SRID
# -I : 创建空间索引
# -C : 添加约束
# -M : 使用 COPY 而非 INSERT
# -t 256x256 : 瓦片大小
# -R : 引用外部文件（不导入数据）
# -r : 常规块
# 分块导入（推荐用于大文件）
raster2pgsql -s 4326 -I -C -M -t 256x256 large_dem.tif public.dem | psql -d gis_db
# 导入目录下所有文件
raster2pgsql -s 4326 -I -C -M -t 256x256 *.tif public.imagery | psql -d gis_db
# 追加数据
raster2pgsql -s 4326 -a -t 256x256 new_data.tif public.dem | psql -d gis_db
# 外部引用（数据保留在文件系统）
raster2pgsql -s 4326 -R -I -C /path/to/dem.tif public.dem_ref | psql -d gis_db
12.2.2 使用 SQL 导入
-- 从文件创建栅格
SELECT ST_FromGDALRaster(pg_read_binary_file('/tmp/dem.tif'));
-- 创建空栅格
SELECT ST_MakeEmptyRaster(
100, -- 宽度（像素）
100, -- 高度（像素）
0, -- 左上角 X
100, -- 左上角 Y
1, -- X 像素大小
-1, -- Y 像素大小
0, -- X 倾斜
0, -- Y 倾斜
4326 -- SRID
);
-- 添加波段
SELECT ST_AddBand(
ST_MakeEmptyRaster(100, 100, 0, 100, 1, -1, 0, 0, 4326),
'32BF'::text, -- 像素类型
0 -- 初始值
);
12.2.3 像素类型
| 类型 | 说明 | 范围 |
|---|---|---|
| 1BB | 1位布尔 | 0-1 |
| 2BUI | 2位无符号整数 | 0-3 |
| 4BUI | 4位无符号整数 | 0-15 |
| 8BSI | 8位有符号整数 | -128 到 127 |
| 8BUI | 8位无符号整数 | 0-255 |
| 16BSI | 16位有符号整数 | -32768 到 32767 |
| 16BUI | 16位无符号整数 | 0-65535 |
| 32BSI | 32位有符号整数 | ±2^31 |
| 32BUI | 32位无符号整数 | 0 到 2^32 |
| 32BF | 32位浮点 | 单精度浮点 |
| 64BF | 64位浮点 | 双精度浮点 |
12.3 栅格元数据访问
12.3.1 基本元数据
-- 获取栅格宽度和高度
SELECT ST_Width(rast), ST_Height(rast) FROM dem LIMIT 1;
-- 获取波段数量
SELECT ST_NumBands(rast) FROM dem LIMIT 1;
-- 获取像素大小
SELECT ST_PixelWidth(rast), ST_PixelHeight(rast) FROM dem LIMIT 1;
-- 获取 SRID
SELECT ST_SRID(rast) FROM dem LIMIT 1;
-- 获取边界框
SELECT ST_Envelope(rast) FROM dem LIMIT 1;
SELECT Box2D(rast) FROM dem LIMIT 1;
-- 获取完整元数据
SELECT ST_MetaData(rast) FROM dem LIMIT 1;
-- 返回: (upperleftx, upperlefty, width, height, scalex, scaley, skewx, skewy, srid, numbands)
-- 获取波段元数据
SELECT ST_BandMetaData(rast, 1) FROM dem LIMIT 1;
12.3.2 波段属性
-- 获取像素类型
SELECT ST_BandPixelType(rast, 1) FROM dem LIMIT 1;
-- 获取无数据值
SELECT ST_BandNoDataValue(rast, 1) FROM dem LIMIT 1;
-- 获取波段统计
SELECT (ST_SummaryStats(rast, 1)).* FROM dem LIMIT 1;
-- 返回: count, sum, mean, stddev, min, max
-- 获取直方图
SELECT (ST_Histogram(rast, 1, 10)).* FROM dem LIMIT 1;
12.4 像素值操作
12.4.1 读取像素值
-- 读取单个像素值（通过行列号）
SELECT ST_Value(rast, 1, 50, 50) FROM dem LIMIT 1;
-- 读取单个像素值（通过坐标）
SELECT ST_Value(rast, ST_SetSRID(ST_MakePoint(116.4, 39.9), 4326)) FROM dem;
-- 读取多个点的值
SELECT
p.name,
ST_Value(d.rast, p.geom) AS elevation
FROM dem d, poi p
WHERE ST_Intersects(d.rast, p.geom);
-- 读取沿线的像素值
SELECT (ST_DumpValues(
ST_Clip(rast, ST_Buffer(line_geom, 100))
)).* FROM dem;
12.4.2 设置像素值
-- 设置单个像素值
UPDATE dem
SET rast = ST_SetValue(rast, 1, 50, 50, 100)
WHERE rid = 1;
-- 设置多个像素值
UPDATE dem
SET rast = ST_SetValues(
rast, 1,
ARRAY[
ROW(10, 10, 100)::geomval,
ROW(11, 10, 101)::geomval,
ROW(12, 10, 102)::geomval
]
)
WHERE rid = 1;
-- 使用几何设置值
UPDATE dem
SET rast = ST_SetValue(rast, 1, polygon_geom, 999)
WHERE ST_Intersects(rast, polygon_geom);
12.4.3 像素值统计
-- 整体统计
SELECT (ST_SummaryStats(rast, 1)).*
FROM dem;
-- 按区域统计
SELECT
d.name,
(ST_SummaryStatsAgg(ST_Clip(r.rast, d.geom), 1, true)).*
FROM dem r, districts d
WHERE ST_Intersects(r.rast, d.geom)
GROUP BY d.name;
-- 获取最大最小值
SELECT
MAX((ST_SummaryStats(rast, 1)).max) AS global_max,
MIN((ST_SummaryStats(rast, 1)).min) AS global_min
FROM dem;
12.5 栅格分析
12.5.1 栅格代数（Map Algebra）
-- 单栅格表达式
SELECT ST_MapAlgebra(rast, 1, '[rast] * 2') FROM dem;
-- 双栅格表达式
SELECT ST_MapAlgebra(
r1.rast, 1,
r2.rast, 1,
'[rast1] + [rast2]'
)
FROM raster1 r1, raster2 r2
WHERE ST_Intersects(r1.rast, r2.rast);
-- 使用回调函数
CREATE OR REPLACE FUNCTION custom_callback(
value double precision[][][],
pos integer[][],
VARIADIC userargs text[]
)
RETURNS double precision AS $$
BEGIN
RETURN value[1][1][1] * 2; -- 简单示例：值乘以2
END;
$$ LANGUAGE plpgsql IMMUTABLE;
SELECT ST_MapAlgebra(
rast, 1,
'custom_callback(double precision[], integer[], text[])'::regprocedure,
'32BF'
) FROM dem;
-- NDVI 计算示例（假设波段4是近红外，波段3是红光）
SELECT ST_MapAlgebra(
rast, ARRAY[4, 3],
'([rast.4] - [rast.3]) / ([rast.4] + [rast.3])::float'
) AS ndvi
FROM satellite_imagery;
12.5.2 邻域分析
-- 焦点统计（窗口分析）
-- 计算3x3窗口的平均值
SELECT ST_MapAlgebra(
rast, 1,
'32BF',
'st_mean4ma(double precision[][][], text[], text[])'::regprocedure,
NULL,
1, 1 -- 窗口大小：左右各1像素
) FROM dem;
-- 坡度计算
SELECT ST_Slope(rast, 1, '32BF') FROM dem;
-- 坡向计算
SELECT ST_Aspect(rast, 1, '32BF') FROM dem;
-- 山体阴影
SELECT ST_HillShade(rast, 1, '32BF', 315, 45, 1) FROM dem;
-- 参数：栅格，波段，输出类型，方位角，高度角，高程夸张因子
12.5.3 重分类
-- 离散重分类
SELECT ST_Reclass(
rast, 1,
'0-100:1, 100-200:2, 200-300:3, 300-:4', -- 重分类规则
'8BUI', -- 输出类型
0 -- 无数据值
) FROM dem;
-- 连续重分类
SELECT ST_Reclass(
rast, 1,
'0-1000:0-100', -- 线性映射
'8BUI'
) FROM dem;
-- 创建土地覆盖分类
SELECT ST_Reclass(
rast, 1,
'0-10:1, 10-20:2, 20-30:3, 30-40:4, 40-:5',
'8BUI',
0
) FROM land_cover;
12.6 栅格裁剪与合并
12.6.1 ST_Clip（裁剪）
-- 使用几何裁剪栅格
SELECT ST_Clip(rast, geom)
FROM dem, districts
WHERE ST_Intersects(rast, geom) AND districts.name = '东城区';
-- 带无数据值的裁剪
SELECT ST_Clip(rast, geom, -9999, true) FROM dem, districts;
-- 裁剪并保留整个瓦片
SELECT ST_Clip(rast, geom, true) FROM dem, districts;
-- 裁剪多个波段
SELECT ST_Clip(rast, ARRAY[1, 2, 3], geom) FROM multi_band_raster, districts;
12.6.2 ST_Union（合并）
-- 合并栅格
SELECT ST_Union(rast) FROM dem WHERE region = 'north';
-- 使用不同的合并策略
SELECT ST_Union(rast, 'MAX') FROM dem; -- 取最大值
SELECT ST_Union(rast, 'MIN') FROM dem; -- 取最小值
SELECT ST_Union(rast, 'MEAN') FROM dem; -- 取平均值
SELECT ST_Union(rast, 'FIRST') FROM dem; -- 取第一个值
SELECT ST_Union(rast, 'LAST') FROM dem; -- 取最后一个值
SELECT ST_Union(rast, 'SUM') FROM dem; -- 求和
SELECT ST_Union(rast, 'COUNT') FROM dem; -- 计数
-- 创建镶嵌影像
CREATE TABLE mosaic AS
SELECT ST_Union(rast, 'MEAN') AS rast
FROM imagery_tiles;
12.6.3 ST_Tile（分块）
-- 将大栅格分割成小瓦片
SELECT ST_Tile(rast, 256, 256) FROM large_raster;
-- 分割并保存
CREATE TABLE dem_tiles AS
SELECT row_number() OVER () AS tid,
ST_Tile(rast, 256, 256) AS rast
FROM dem;
12.7 矢量栅格转换
12.7.1 ST_AsRaster（矢量转栅格）
-- 多边形转栅格
SELECT ST_AsRaster(
geom,
100, -- 宽度
100, -- 高度
'8BUI', -- 像素类型
1, -- 值
0 -- 无数据值
) FROM districts WHERE name = '东城区';
-- 使用参考栅格
SELECT ST_AsRaster(
geom,
ref_rast, -- 参考栅格（用于获取分辨率和范围）
'8BUI',
1,
0
) FROM districts, dem LIMIT 1;
-- 点密度栅格化
SELECT ST_AsRaster(
ST_Collect(geom),
0.01, -- 像素宽度
0.01, -- 像素高度
'16BUI',
0,
0
) FROM poi;
12.7.2 ST_Polygon/ST_DumpAsPolygons（栅格转矢量）
-- 栅格转多边形
SELECT ST_Polygon(rast, 1) FROM dem;
-- 分别转换每个值
SELECT (ST_DumpAsPolygons(rast, 1)).* FROM classified_raster;
-- 等值线
SELECT ST_Contour(rast, 1, interval := 100) FROM dem;
-- 提取特定值的区域
SELECT ST_Polygon(ST_Reclass(rast, 1, '500-:1, 0-500:0', '8BUI'), 1)
FROM dem;
12.8 栅格导出
12.8.1 导出为文件
-- 导出为 GeoTIFF
SELECT ST_AsTIFF(rast) FROM dem WHERE rid = 1;
-- 带压缩的导出
SELECT ST_AsTIFF(rast, 'LZW') FROM dem WHERE rid = 1;
-- 导出为 JPEG
SELECT ST_AsJPEG(rast) FROM imagery WHERE rid = 1;
-- 导出为 PNG
SELECT ST_AsPNG(rast) FROM imagery WHERE rid = 1;
-- 使用 GDAL 驱动导出
SELECT ST_AsGDALRaster(rast, 'GTiff', ARRAY['COMPRESS=LZW']) FROM dem;
-- 保存到文件（需要超级用户权限）
COPY (SELECT ST_AsTIFF(ST_Union(rast)) FROM dem) TO '/tmp/output.tif' WITH (FORMAT binary);
12.8.2 导出为其他格式
-- 导出为 WKB
SELECT ST_AsBinary(rast) FROM dem LIMIT 1;
-- 导出为十六进制
SELECT ST_AsHexEWKB(rast) FROM dem LIMIT 1;
-- 获取支持的驱动
SELECT * FROM ST_GDALDrivers() WHERE can_write = true;
12.9 性能优化
12.9.1 分块存储
-- 创建分块表
CREATE TABLE dem_tiled AS
SELECT row_number() OVER () AS rid,
ST_Tile(rast, 256, 256) AS rast
FROM dem_original;
-- 创建空间索引
CREATE INDEX idx_dem_tiled_rast ON dem_tiled USING GIST (ST_ConvexHull(rast));
-- 添加约束
SELECT AddRasterConstraints('dem_tiled', 'rast');
12.9.2 外部栅格
-- 引用外部文件（数据不存储在数据库中）
CREATE TABLE dem_outdb AS
SELECT rid, ST_AddBand(
ST_MakeEmptyRaster(rast),
'/path/to/dem.tif'::text,
ARRAY[1]
) AS rast
FROM dem_metadata;
-- 检查是否为外部栅格
SELECT ST_BandIsNoData(rast, 1), ST_BandPath(rast, 1) FROM dem_outdb;
12.10 本章小结
本章详细介绍了 PostGIS Raster：
- 基础概念：栅格数据结构、像素类型
- 数据导入：raster2pgsql、SQL 方法
- 元数据访问：栅格和波段属性
- 像素操作：读取、设置、统计
- 栅格分析：代数运算、邻域分析、重分类
- 裁剪合并：ST_Clip、ST_Union、ST_Tile
- 格式转换：矢量栅格互转
- 性能优化：分块存储、外部栅格
12.11 下一步
在下一章中，我们将学习拓扑处理，包括：
- 拓扑数据模型
- 拓扑构建
- 拓扑编辑
- 拓扑验证
相关资源：
浙公网安备 33010602011771号