#妙笔生花创作挑战#

数字高程模型（Digital Elevation Model)，简称DEM，是描述地表高程和形态的数字数据集，也就是用一组有序数值阵列形式表示地面高程起伏形态的一种实体地面模型。DEM数据在测绘、气象、地质、军事、土地资源规划、应急管理、地形分析、洪水模拟等领域都有较深的应用。Python中有多个库可以用于处理GeoTIFF格式的数据，常见的有GDAL和Rasterio等。下面我们介绍这两种方法的使用方法。

一、DEM数据常见的格式

DEM数据可以使用多种文件格式存储，常见的格式包括GeoTIFF（.tif）、ASCII Grid（.asc）和Esri Grid（.adf）等，当然不同的格式有不同的存储结构和编码方式。GeoTIFF是一种常见的地理信息格式，广泛应用于DEM数据的存储和传输，本教程重点谈如何读GeoTIFF格式文件的方法。

二、使用GDAL模块

GDAL（Geospatial Data Abstraction Library）是一个用于读取和写入地理空间数据的开源库。它支持多种DEM数据的格式，包括GeoTIFF。

安装GDAL库。

pip install gdal或者conda install gdal -c conda-forge或者mamba install gdal -c conda-forge

使用GDAL库打开GeoTIFF格式的DEM数据，并读取其中的高程信息。示例代码如下：

from osgeo import gdal# 打开DEM数据dataset = gdal.Open('srtm_58_06.tif')# DEM数据的宽度和高度width = dataset.RasterXSizeheight = dataset.RasterYSize# DEM数据的地理空间参考信息geotransform = dataset.GetGeoTransform()projection = dataset.GetProjection()# 数据的平面四角，计算数据的起止经纬度Xmin = geotransform[0] Ymin = geotransform[3]Xmax = geotransform[0] + width * geotransform[1] + width * geotransform[2]Ymax = geotransform[3] + height * geotransform[4] + height * geotransform[5]# DEM数据的高程信息band = dataset.GetRasterBand(1)elevation = band.ReadAsArray(0, 0, width, height)

使用gdal.Open函数，打开了文件名为：srtm_58_06.tif的GeoTIFF格式的DEM数据。

使用RasterXSize和RasterYSize属性，得到宽度和高度。

使用GetGeoTransform和GetProjection方法，得到地理空间参考信息。

使用geotransform的数据，计算起止经纬度。

使用GetRasterBand和ReadAsArray方法读取高程信息。只有一个通道，所以参数为1。注意数据的通道的标号是从1开始的。

三、使用Rasterio模块

Rasterio是一个基于GDAL库的Python库，也就是Rasterio的运行需要GDAL的库的支持。

安装Rasterio库。

pip install rasterio或者conda install rasterio -c conda-forge或者mamba install rasterio -c conda-forge

使用Rasterio库打开GeoTIFF格式的DEM数据，并读取其中的高程信息。示例代码如下：

import rasterio# 打开DEM数据dataset = rasterio.open('srtm_58_06.tif'')# DEM数据的宽度和高度width = dataset.widthheight = dataset.height# DEM数据的地理空间参考信息transform = dataset.transformcrs = dataset.crs   # DEM数据的空间边界框,平面四角           bounds = dataset.boundsXmin = bounds.leftYmin = bounds.topXmax = bounds.rightYmax = bounds.bottom# DEM数据的高程信息，elevation = dataset.read(1)

使用rasterio.open函数，打开文件名为srtm_58_06.tif的GeoTIFF格式的DEM数据。

使用width和height属性，得到宽度和高度。

使用transform和crs属性，得到地理空间参考信息。

使用bounds性，得到的空间边界框,也就是平面四角坐标。

使用read方法读取高程信息。只有一个通道，所以参数为1。注意数据的通道的标号是从1开始的。

从上可以看出，rasterio库提供比GDAL库更简单和友好的接口来处理栅格数据。

四、简单的可视化

获取到DEM数据之后，我们可以对其进行各种处理和分析。下面演示将DEM数据可视化。

DEM数据，地形图示例代码如下：

# 加载 gdalfrom osgeo import gdalimport matplotlib.pyplot as plt# 打开DEM数据dataset = gdal.Open('srtm_58_06.tif')# DEM数据的宽度和高度width = dataset.RasterXSizeheight = dataset.RasterYSize# DEM数据的地理空间参考信息geotransform = dataset.GetGeoTransform()projection = dataset.GetProjection()print(geotransform)# 数据的平面四角，计算数据的起止经纬度Xmin = geotransform[0] Ymin = geotransform[3]Xmax = geotransform[0] + width * geotransform[1] + width * geotransform[2]Ymax = geotransform[3] + height * geotransform[4] + height * geotransform[5]# DEM数据的高程信息band = dataset.GetRasterBand(1)elevation = band.ReadAsArray(0, 0, width, height)#创建名为elevation的图窗plt.figure('elevation')#绘制海拔图，配色为rainbowplt.imshow(elevation, extent = [Xmin, Xmax, Ymax, Ymin], cmap = 'rainbow')#绘制海拔图例plt.colorbar()plt.show()

运行后，地形图如下：

有了DEM数据，上例展示了地形图的展示。我们可以计算DEM数据的坡度和坡向，生成等高线图、、洪水模拟、淹没模型展示等。这些知识和技巧我们后面有机会进一步讨论和研究。