14.7. 处理数字地形模型¶
gdaldem 命令致力于分析数字地形模型，建议生成七种类型的产品：山体阴影、斜坡、面貌、色彩浮雕、地形崎岖指数（TRI）、
地形位置指数（TIP）和粗糙度地图。
14.7.1. 数字地形模型的准备¶
SRTM数据以WGS 84形式传输，坐标单位为度（纬度、经度），高度单位为米。默认情况下， gdaldem 应用程序假设x、y和z中的
单位相同。因此，建议将此栅格重新投影到Lambert 93中。
QGIS处理步骤
1.打开 DTM
在QGIS中：
仅打开以下文件：
srtm_37_03.tif
2.重投影
这里的目标是通过指定重新采样方法（立方）和输出空间分辨率（90 m）来扭曲 DTM 。
在菜单栏中：
在“Warp”窗口中单击 Raster > Projections > Warp (Reproject)… :
设置如下：
输入文件: srtm_37_03
输出文件: srtm_37_03_l93.tif
来源SRS: EPSG:4326
目标SRS: EPSG:2154
单击编辑按钮，
在命令版本字段中插入以下选项：-tr 90 90
单击“OK”
表2.13. DTM 准备
14.7.2. 斜坡和坡向¶
斜坡是数字地形模型的一阶求导，是像素在x和y方向上的变化率的估计，作为其邻近区域的函数（见图2.8）。它以指数或百分比表示。 该方位对应于斜坡方向。实现了两种方法，Zevenbergen和Thorne [ZEVENBERGEN87] 和 Horn [HORN81] 。后者默认提供。
ZEVENBERGEN L.W., THORNE C.R., “地表地形的定量分析”， 地球表面过程和地貌，第12卷，第1期，第47-56页，1987年。
HORN B.K., “山丘阴影和反射率图，” 《IEEE学报》，第69卷，第1期，第14–47页，1981年。
x和y两个方向的变化率Dx和DY重新计算如下：
Dx = ((c + 2f + i) - (a + 2d + g)/(8 x res_pixel)
Dx =（（c + 2f + i）-（a + 2d + g）/（8 x res_pixel）
Dy = ((g + 2h + i) - (a + 2b + c))/(8 x res_pixel)
DY =（（g +2 h + i）-（a + 2b + c））/（8 x res_pixel）
其中res_pixel是栅格的空间分辨率。
在示例中，Dx和DY等于：
Dx = ((376 + 2 x 398 + 410) - (368 + 2 x 377 + 400))/(8 x 90)
Dx =（（376 + 2 x 398 + 410）-（368 + 2 x 377 + 400））/（8 x 90）
Dy = ((400 + 2 x 406 + 410) - (368 + 2 x 369 + 376))/(8 x 90)
DY =（（400 + 2 x 406 + 410）-（368 + 2 x 369 + 376））/（8 x 90）
斜坡计算：
以弧度为单位：
slope_rad = atan((Dx² + Dy²)0.5)
slope_rad = atan（（Dx² + Dy²）0.5）
slope_rad = 0.2085
slope_rad = 0.2085
以度为单位：
slope_deg = pente_rad x 180/π
slope_Degree = pente_rad x 180/pi
slope_deg = 11.94°
slope_Degree = 11.94°
坡向计算：
以弧度为单位： expo_rad = atan2(Dy, -Dx) expo_rad = 1.976
以度（三角形）为单位：
expo_deg = expo_rad x 180 / π expo_deg = 113.198
expo_Degree = expo_rad x 180 /pi expo_Degree = 113.198
转换为方位角：
如果expo_dev> 90，则：
expo_azimut = 450 - expo_Dege 其他：
expo_azimut = 90 - expo_dev 因此：expo_azimut = 450 - 113.198_azimut = 336.80°，即N-N-W暴露
图2.8.斜坡和长宽计算示例
QGIS处理步骤
1.打开 DTM，在QGIS中：
仅打开以下文件：
srtm_37_03_l93.tif
2.斜率计算
在菜单栏中：
在“ DEM ”窗口中单击 Raster > Analysis > DEM (Terrain Models)… :
设置如下：
输入文件： srtm_37_03_l93
输出文件： srtm_37_03_l93_slope.tif
模式：坡度
单击“OK”
3.坡向计算
在菜单栏中：
单击 Raster > Analysis > DEM (Terrain Models)…
在“ DEM ”窗口中：
设置如下：
输入文件： srtm_37_03_l93
输出文件： srtm_37_03_l93_aspect.tif
模式： Aspect
单击“OK”
4.结果分析
srtm_37_03_l93_slope.tif 栅格是以指数表示的斜坡地图。 比例参数（默认情况下，等于1）在水平和垂直单位不相同的情况下很有用， 例如地理坐标中海拔值以米为单位的图像。
srtm_37_03_l93_aspect.tif 栅格是曝光图，其值在0和360之间，代表斜坡方向的方位角。 默认情况下，平面区域编码为 无数据 （-9999）。
DTM 斜坡方面
5.对应GDAL 命令上面使用的斜坡工具相当于在终端中启动以下命令：
gdaldem slope srtm_37_03_l93.tif srtm_37_03_l93_slope.tif -s 1.0 -of
GTiff
其中：
[slope]: 表示期望的乘积，这里是斜率，
[-s ]: 设置垂直单位与水平单位的比率，
[-of]: 设置输出格式。
一些有用的选项：
[-p]: 将斜率表示为百分比（默认情况下，以度为单位）
上面使用的方面工具相当于在终端中启动以下命令：
gdaldem aspect srtm_37_03_l93.tif srtm_37_03_l93_aspect.tif -of GTiff
其中：
[aspect]: 表示预期的产品，此处为曝光量，
[-of]: 设置输出格式。
一些有用的选项：
[-trigonometric]: 返回三角形(0° for East, 90° North, 180° for West and 270° for South),
[-zero_for_flat]: 对于平坦区域，返回0。
有关更多信息，请参阅以下网页： http://www.gdal.org/gdaldem.html 。
表2.14.从 DTM 计算坡度和坡向
14.7.3. 山体阴影¶
山体阴影通常用于制图中，以给人一种深度的印象并促进浮雕的解释。该技术向像素指定与假设光源相关的照明值 [BURROUGH98] 。 它需要四个输入参数：光源的方位角和仰角、斜率和纵横比（见图2.9）。 结果显示灰色阴影（8位编码，即0到255之间的256个值），反映直接阴影和光线。
BURROUGH P.A., MCDONNELL R., MCDONNELL R.A. 等人, 地理信息系统原理, 牛津大学出版社，1998年。
考虑光源的以下参数：
方位角：az_deg = 315°
仰角： el_deg = 45°
方位角为：
径向： az_rad = (450 – az_deg) × π / 180 az_rad = 2.356
天角为：
以度为单位：
ze_deg = 90 - el_deg = 45°
ze_Degree = 90 - el_Degree = 45°
以弧度为单位： ze_rad = ze_deg × π / 180 ze_rad = 0.7854
山体阴影计算：
hillshade = 255.0 x（（cos（ze_rad）x cos（slope_rad））+（sin（ze_rad）x sin（slope_rad）x cos（az_rad - expo_rad） hillshade = 255.0 x（（cos（0.7854）x cos（0.2085））+（sin（0.7854）x sin（0.2085）x cos（2.356- 1.976） hillshade = 211.05，即8位编码的211。
图2.9.山体阴影计算示例
QGIS处理步骤
1.山体计算
在菜单栏中：
单击 Raster > Analysis > DEM (Terrain Models)…
在“ DEM ”窗口中：
设置如下：
输入文件： srtm_37_03_l93
输出文件： srtm_37_03_l93_hillshade.tif
模式： Hillshade
单击“OK”
2.结果分析
srtm_37_03_l93_hillshade.tif 栅格通过应用位于315°方位角和45°俯仰角的光源的修辞照明来突出显示浮雕。 因此，面向东的斜坡是有阴影的。
DTM
山体阴影
2.相应的GDAL 命令
上面使用的hillshade工具相当于在终端中启动以下命令：
gdaldem hillshade srtm_37_03_l93.tif
srtm_37_03_l93_hillshade.tif -z 1.0 -s 1.0 -az 315.0 -alt 45.0 –of GTiff
其中：
[hillshade]: 表示期望的产品，这里是hillshade，
[-z]: 设置垂直夸张（默认为1）
[-s ]: 设置垂直单位与水平单位的比率，
[-az]: 指定光源的方位角，
[-alt]: 指定光源的仰角，
[-of]: 设置输出格式。
有关更多信息，请参阅以下网页： http://www.gdal.org/gdaldem.html 。
表2.15.从 DTM 计算Hillshade
14.7.4. 粗糙度和地形位置指数¶
TRI是地形不均匀性的衡量标准 [RILEY99] , [WILSON07] 。该方法计算每个像素与其八个邻居的垂直差的平均值（见图2.10）。 与TRI不同，粗糙度指数将邻近像素的最大垂直差异指定。
地形位置指数（TIP）是从中心像素值中减去八个邻居的平均值 [WEISS01] , [WILSON07]_(图2.10) 。 正值代表升高的区域（顶峰、峰等）与邻近地区相关；负值代表覆盖区（山谷等）； 接近零的值表示平坦区域或具有恒定斜坡的区域。
RILEY S.J.,“量化地形异质性的指数”，《山间科学杂志》，第5卷，第1-4期，第23-27页，1999年。
WEISS A.,“地形位置和地形分析”，ESRI用户会议，加利福尼亚州圣地亚哥，第200卷，2001年7月。
WILSON M.F., O’CONNELL B., BROWN C等人, “用于大陆斜坡栖息地测绘的多波束测深数据的多尺度地形分析”，海洋大地测量学，第30卷， 第1–2号，第3–35页，2007年。
TRI计算如下：
TRI = (|a - e|+|b - e|+|c - e|+|d - e|+|f - e|+|g - e|+|h - e|+|i - e|) / 8
在示例中，中心像素TRI为：
TRI = (|205-196|+|195-196|+|197-196|+|209-196|+|192-196|+|212196|+|195-196|+|185-196|)/8
TRI = 56/8
TRI = 7 meters
在示例中，中心像素粗糙度为：
Roughness = g – i
Roughness = 212-185
Roughness = 27 meters
TPI is calculated as follows:
TPI = e - ((a + b + c + d + f + g + h + i) / 8)
在示例中，中心像素TIP为：
TPI = 196-((205+195+197+209+192+212+195+185)/8)
TPI = 196-198.75
TPI = -2.75
图2.10.粗糙度和地形位置指数计算示例
QGIS处理步骤
TRI计算
在菜单栏中：
在" DEM "窗口中单击 Raster > Analysis > DEM (Terrain Models)… :
设置如下：
输入文件： srtm_37_03_l93
输出文件： srtm_37_03_l93_TRI.tif
模式： TRI (Terrain Ruggedness Index)
单击“OK”
注意：要计算粗糙度指数，请选择"粗糙度"模式。
TIP计算
在菜单栏中：
在" DEM "窗口中单击 Raster > Analysis > DEM (Terrain Models)… :
设置如下：
输入文件： srtm_37_03_l93
输出文件： srtm_37_03_l93_TPI.tif
模式：: 地形位置指数(TPI)
单击“OK”
3.结果分析
TRI值与斜坡强相关。它们在平坦地区地势较低，在陡峭地区地势较高。这种方法的缺点是没有考虑与不同斜坡方向相关的变化性。 为了纠正此缺陷，通常建议使用矢量坚固性测量（VRM） [SAPPINGTON07] 。不幸的是，GDAL 中不存在后者。 不过，它可以在Grass和Saga中使用。至于TIP值，它们使识别地形形式和地形突变成为可能。这两个指数通常用于生态学和地貌学， 特别是用于河流过程的研究。
SAPPINGTON J.M., LONGSHORE K.M., THOMPSON D.B., “量化景观崎岖度以进行动物栖息地分析：莫哈韦沙漠大角羊的案例研究”，《野生动物管理杂志》， 第71卷，第5期，第1419-14262007页。
DTM TRI TPI
4.相应的GDAL 命令
上面使用的TRI工具相当于在终端中启动以下命令：
gdaldem TRI srtm_37_03_l93.tif srtm_37_03_l93_TRI.tif -of GTiff
其中：
[TRI]: 表示预期产品，此处为TRI
上面使用的TIP工具相当于在终端中启动以下命令：
gdaldem TPI srtm_37_03_l93.tif srtm_37_03_l93_TPI.tif -of GTiff where:
[TPI]: 表示预期产品，此处为TPI
有关更多信息，请参阅以下网页： http://www.gdal.org/gdaldem.html 。
表2.16.从 DTM 计算TRI和TPI
14.7.5. 浮雕着色¶
GDAL提供了一个最终工具，根据表格，该工具将地形与海拔相关联。与前几节描述的工具相反，结果具有制图文档等美学目标。
QGIS处理步骤
1.颜色表创建
颜色表是一个文本文件，通常每行包含四个或五个列：海拔值、红色、绿色和蓝色的相应颜色（值在0和255之间）， 以及可选的Alpha通道的值（透明度）。
可以为 无数据 像素分配颜色。在这种情况下，海拔值将被 “nv” 替换。海拔也可以用百分比表示：0%对应于 DTM 的最小值， 100%对应于最大值。
使用文本编辑器（例如Notepad++）：
根据以下示例（左列）构建海拔高度的颜色表：
将结果保存为：color_terrain.txt
2.浮雕上色
在菜单栏中：
在 DEM 窗口中单击 Raster > Analysis > DEM (Terrain Models)… :
Set as follows:
输入文件： srtm_37_03_l93
输出文件： srtm_37_03_l93_color.tif
模式：彩色浮雕
颜色配置文件： color_terrain.txt
注意：有一个名为 terrain.txt 的默认文件。然而，颜色范围不适合本示例中处理的SRTM文件。
单击“OK”
3.结果分析
输出图像包含红色、绿色和蓝色三个分量的三个频段，而不是海拔。 默认情况下，color_terrain.txt 文件中提到的海拔之间的颜色进行插值以获得渐变。
DTM 色彩浮雕
4.颜色和浮雕阴影
在制图中，海拔颜色梯度和灰化（山体阴影上的阈值）的结合是一种经常用于视觉恢复地形的技术。 可以在QGIS中生成此类图像。
颜色浮雕明暗处理浮雕贴图
5.相应的GDAL 命令
终端中的命令：
gdaldem color-relief srtm_37_03_l93.tif color_terrain.txt
srtm_37_03_l93_color.tif -of GTiff
其中：
[color-relief]: 表示预期的产品，这里是彩色浮雕，
[-of]: 设置输出格式。
一些有用的选项：
[-alpha]: 将alpha通道添加到输出栅格数据，
[-exact_color_entry]: 在颜色表中搜索时使用严格匹配。如果没有找到匹配的颜色条目则高度以黑色显示，
[-nearest_color_entry]: 使用与颜色表中最近的条目相对应的颜色。
有关更多信息，请参阅以下网页：http://www.gdal.org/gdaldem.html。
表2.17.浮雕上色