PostGIS实现DEM分析之坡度计算【ST_Slope】
全部文章内容请转公众号【时空智能计算】及交流群867520770，原创不易，求关注支持，更多开源GIS相关知识技能分享，免费提供学习问答交流。
一、函数概述
ST_Slope用于实现 DEM 坡度计算。它通过分析每个像素与其相邻像素的高程变化，返回地形的倾斜程度，结果可选择以度数（默认）、弧度或百分比表示。该函数支持自定义垂直缩放比例和 NoData 插值，是地形分析、城市规划及环境建模的基础工具。
二、参数解析
raster ST_Slope(raster rast,-- 输入DEM栅格integerband=1,-- 栅格波段（默认第1波段）textpixeltype='32BF',-- 输出像素类型（默认32位浮点型）textunits='DEGREES',-- 输出单位：DEGREES（度）、RADIANS（弧度）、PERCENT（百分比）doubleprecisionscale=1.0,-- 垂直缩放比例（默认1，用于调整垂直单位与水平单位的比例）booleaninterpolate_nodata=FALSE-- 是否插值处理NoData);参数说明：
1. units 参数：
- DEGREES（默认）：坡度以度数表示（0°-90°，0° 为水平，90° 为垂直）。
- RADIANS：以弧度表示（0-π/2）。
- PERCENT：以百分比坡度表示（上升高度 / 水平距离 ×100%，如 10% 表示每 100 米水平距离上升 10 米）。
2. 垂直缩放（scale）：
- 用于校正垂直单位与水平单位的差异：
- 经纬度数据（单位：度）：scale=111120（米 / 度，适用于米制垂直单位）。
- 英尺 / 度经纬度数据：scale=370400。
- 投影坐标系（如 UTM，单位：米）：scale=1.0（单位一致，无需缩放）。
3. NoData 处理：
- interpolate_nodata=TRUE时，使用反距离加权插值（ST_InvDistWeight4ma）填充空洞，提升结果连续性。
三、算法原理
ST_Slope 使用 3×3 窗口分析中心像素与 8 邻域像素的高程差异，计算公式为：
slope = arctan( (z_factor * dz) / dx )其中：
- dx = 像素水平间距（由栅格元数据的scalex确定）。
- dz = 中心像素与邻域像素的高程差。
- z_factor = scale / dx（将垂直单位转换为与水平单位一致）。
单位转换：
- 从弧度转换为度数：slope_degrees = slope_radians × (180/π)。
- 转换为百分比坡度：slope_percent = tan(slope_radians) × 100%。
四、典型示例
示例 1：默认参数计算坡度（度数）
WITHdemAS(SELECTST_SetValues(ST_MakeEmptyRaster(5,5,0,0,1,-1),1,1,1,ARRAY[[10,10,10,10,10],[10,15,15,15,10],[10,15,20,15,10],[10,15,15,15,10],[10,10,10,10,10]]::doubleprecision[][])ASrast)SELECTST_Slope(rast,1,'32BF','DEGREES')ASslope_degreesFROMdem;输出：中心像素（20）的坡度为约 26.57°，周围像素按高程变化计算对应坡度。
示例 2：经纬度数据垂直缩放
SELECTST_Slope(rast,band :=1,units :='DEGREES',scale :=111120,-- 米/度（经纬度数据垂直单位为米）interpolate_nodata :=TRUE)ASslope_rastFROMdem_latlon;场景：处理全球范围的经纬度 DEM 数据，正确转换垂直单位与水平单位比例。
示例 3：输出百分比坡度
SELECTST_Slope(rast,units :='PERCENT',scale :=1.0-- 投影坐标系（米/米）)ASslope_percentFROMdem_utm;效果：输出值为百分比形式，如坡度 20% 表示每 100 米水平距离上升 20 米。
五、输出值与地形等级对应
| 坡度范围（度数） | 坡度等级 | 典型场景 |
|---|---|---|
| 0°-3° | 极缓坡 | 平原、机场跑道 |
| 3°-10° | 缓坡 | 农田、高尔夫球场 |
| 10°-20° | 中坡 | 山地森林、乡村道路 |
| 20°-30° | 陡坡 | 山区、梯田 |
| >30° | 极陡坡 | 悬崖、冰川边缘 |
六、性能优化与注意事项
1. 边界处理：
- 栅格边缘像素因邻域不完整，结果可能不准确，建议通过ST_Expand扩展栅格边界：
WITHexpanded_demAS(SELECTST_Expand(rast,1)ASrast-- 扩展1个像素FROMdem_table)SELECTST_Slope(rast)FROMexpanded_dem;2. 数据类型选择：
- 输出类型推荐32BF（浮点型）以保留小数精度，若仅需等级划分可转为8BUI（整数型）：
ST_Slope(rast,pixeltype :='8BUI')3. 并行计算：
- 结合 PostgreSQL 并行查询加速大规模栅格处理：
SETmax_parallel_workers_per_gather=4;SELECTST_Slope(rast)FROMlarge_dem;4. NoData 值校验：
- 输入 DEM 需通过ST_SetBandNoDataValue正确设置 NoData 值，避免计算错误。
七、应用场景
1. 城市规划与工程：
- 识别陡坡区域（>25°），限制建筑开发或规划特殊工程措施。
2. 农业与林业：
- 坡度 > 15° 的区域易发生水土流失，需优先实施水土保持方案。
3. 灾害风险评估：
- 结合坡度（ST_Slope）与植被覆盖数据，预测滑坡、泥石流易发区。
4. 路径规划：
- 自动驾驶或无人机航线规划中避开陡坡区域（如 > 30°）。
八、总结
通过ST_Slope，可精确量化地形倾斜程度，为各类空间分析提供关键参数。合理设置垂直缩放和单位，结合插值处理 NoData 区域，能显著提升计算结果的准确性和适用性。