【GDAL-Python】3-在Python中使用GDAL处理数字高程模型DEM

1-介绍

1.1 主要内容

（1）教程内容：使用GDAL处理数字高程模型（DEM)得到一些分析的结果，如坡度（slope)，坡向（aspect)，山体阴影（hillshade)
（2）四种技术路线：
1) 终端运行gdal可执行命令；
2）在Python脚本中调用gdal可执行命令；
3）在Python脚本中调用gdal.DEMProcessing
4）使用richdem库中的TerrainAttribute计算

PS:Python脚本中调用gdal可执行命令有两种方式，比如利用os.system(cmd），或者子进程库subprocess.call(cmd.split())

（3）视频教程：B站对应教程：3-在Python中使用GDAL处理数字高程模型DEM

1.2 坡度、坡向、山体阴影

（1）坡度：坡度是指地表单元陡缓的程度，通常通过坡面的垂直高度与水平距离的比值来表示，这个比值可以用百分比或角度（正切值）来衡量。数值越低，表示地形越平坦；数值越高，表示地势越陡峭
（2）坡向：坡向则是指地形坡面的朝向，定义为坡面法线在水平面上的投影的方向，它用于识别表面上某一位置处的最陡下坡方向。
（3）山体阴影：山体阴影是一种技术，用于可视化由光源和高程表面的坡度和坡向确定的地形，通过模拟光线和阴影落在表面上的方式，在平面地图上描绘景观的3D表面

2-代码实现

2.1 数据介绍

（1）影像信息

（2）经过global mapper软件渲染显示效果，左侧颜色条表示高程信息。

2.2 代码实现

#！/usr/bin/env python3
# -*- coding:utf-8 -*-
import os
import subprocess
from osgeo import gdal
import matplotlib.pyplot as plt
#======================================1.计算坡度===============================
# 处理dem.tif生产坡度文件slope2.tif,采用python脚本中调用gdaldem可执行文件的形式
cmd='gdaldem slope Tile_+017_+010\Tile_+017_+010_dem.tif slope2.tif -compute_edges'
# 形式1
os.system(cmd)
# 形式2
# subprocess.call(cmd.split()))
#显示效果
slp2=gdal.Open('slope2.tif')
slp2Array=slp2.GetRasterBand(1).ReadAsArray()

plt.figure()
plt.imshow(slp2Array)
plt.colorbar()
plt.show()

#======================================2.计算坡向===============================
cmd='gdaldem aspect Tile_+017_+010\Tile_+017_+010_dem.tif aspect2.tif -compute_edges'
os.system(cmd)
# subprocess.call(cmd.split()))
aspect2=gdal.Open('aspect2.tif')
aspect2Array=slp2.GetRasterBand(1).ReadAsArray()
plt.imshow(aspect2Array,cmap='magma')
plt.colorbar()
plt.show()
#======================================3.计算山体阴影===============================
# 采用技术路线3实现：gdal.DEMProcessing
dem=gdal.Open('Tile_+017_+010\Tile_+017_+010_dem.tif')
hill=gdal.DEMProcessing('hillshade.tif',dem,'hillshade',computeEdges=True)
hillshade=gdal.Open('hillshade.tif')
hillArray=hillshade.GetRasterBand(1).ReadAsArray()
plt.imshow(hillArray,cmap='gray')
plt.colorbar()
plt.show()

2.3 效果显示

（1）坡度显示效果

（2）坡向显示效果

（3）山体阴影显示效果

3.参考资料

3.1 使用richdem库中的TerrainAttribute计算坡度、坡向、山体阴影

# using richdem
dem = rd.LoadGDAL("dem.tif")
slp3 = rd.TerrainAttribute(dem, attrib="slope_degrees") # replace "slope_degrees" with "slope_riserun",  "aspect" ...
rd.SaveGDAL("slope3.tif", slp3) 

visualize (example)
plt.figure()
plt.imshow(slp2Array)
plt.colorbar()
plt.show()

TODO:上传实验数据DEM.tif