使用 Python 访问和可视化数字高程模型

数字高程模型 (DEM)代表地形的 3D 表面模型。它通过一系列像元表示连续的地形高程表面，每个像元代表其所在位置（X 和 Y）的海拔 (Z)。数字高程模型仅包含有关地质特征（例如山谷、山脉和滑坡）的海拔信息，不包含任何有关植被或建筑物等特征的高程数据。

准确的高分辨率DEM数据对于灾害地图绘制至关重要，因为它能够提供地形的详细表示，这对于评估自然灾害带来的潜在风险至关重要。这些数据可以帮助科学家测量气温、降水和其他气候相关因素变化对不同海拔地表的影响，从而更好地为气候变化如何影响不同地表的预测模型提供信息。DEM数据还可用于识别面临洪水、山体滑坡和其他极端天气事件风险的区域，这有助于制定关于如何减轻气候变化影响的政策决策。此外，DEM数据有助于创建详细的地图，用于制定疏散计划、警报系统和其他风险管理策略。

美国地质调查局 (USGS) 的EarthExplorer是全球数字高程模型 (DEM) 数据集最受欢迎的来源之一。USGS 提供的数据集分辨率从 10 米到 1 弧秒不等，每 5 到 10 年更新一次。其他全球 DEM 数据集来源包括哥白尼陆地监测服务、国际热带农业中心和全球土地覆盖机构。

本博客总结了如何使用 Python 下载和可视化任何感兴趣的地理区域的公开全球数字高程模型

首先，我们将从人道主义数据交换中心 ( HDX )下载所选区域的 Shapefile 。然后，我们将使用Python 中的Elevation包，该包可以轻松下载、缓存和访问由 NASA 和 NGA 制作并托管在 Amazon S3 上的全球数据集 SRTM 30 米全球 1 弧秒 V003，以及由 CGIAR-CSI 制作的 SRTM 90 米数字高程数据库 v4.1。

此软件包允许我们将选定地理边界的 DEM 下载到本地系统。然后，我们可以使用GDAL或rasterio等库将此 DEM 导入 Python。之后，我们将使用matplotlib探索和可视化 DEM ，并计算坡度和等高线等导数。

通过 HDX 从 OCHA 提取一个国家任意行政区域的多边形

人道协调厅实地信息服务科（FISS）隶属于联合国人道主义事务协调办公室，总部位于瑞士日内瓦。该组织提供通用操作数据集 ，这些数据集是支持人道主义紧急情况初期响应行动和决策的权威参考数据集。其中一个数据集是国家级以下行政边界，该数据集通过人道主义数据交换平台（HDX）根据知识共享政府间组织署名许可提供。

下面的代码可视化了使用folium从 HDX 下载的亚美尼亚的行政边界（本教程以此为例）。

df_armenia_shp = gpd.read_file('arm_adm_2029_shp/arm_admbnda_adm1_2019.shp')
    
m = fl.Map(zoom_start=10, tiles="OpenStreetMap")

for _, r in df_armenia_shp.iterrows():
    sim_geo = gpd.GeoSeries(r["geometry"]).simplify(tolerance=0.001)
    geo_j = sim_geo.to_json()
    geo_j = fl.GeoJson(data=geo_j, style_function=lambda x: {"fillColor": "orange"})
    fl.Popup(r["ADM1_EN"]).add_to(geo_j)
    geo_j.add_to(m)
m

通过 GADMDownloader 从 GADM 中提取的国家行政边界（图片由作者提供）

要下载 DEM，让我们选择 Gegharkunik 省之一。

from shapely.ops import unary_union

selected_province = 'Gegharkunik'
gdf_selected_province = df_armenia_shp[df_armenia_shp['ADM1_EN']==selected_province]
combined_polygon = unary_union(gdf_selected_province['geometry'].values)
combined_polygon

为选定的 ADM 边界提取的多边形（图片由作者提供）

使用高程包将 DEM 下载到本地路径

该包允许下载两个可由属性产品（SRTM1 或 SRTM3）指定的数据集。

• SRTM1（30米分辨率）——NASA JPL（2013）。NASA航天飞机雷达地形测绘任务全球1弧秒数据[数据集]。NASA EOSDIS陆地过程DAAC。访问日期：2023年3月6日，网址：https://doi.org/10.5067/MEaSUREs/SRTM/SRTMGL1N.003
• SRTM3（90 米分辨率）— Jarvis A.、HI Reuter、A. Nelson、E. Guevara，2008 年，填孔无缝 SRTM 数据 V4，国际热带农业中心（CIAT），可从http://srtm.csi.cgiar.org获取。

在本教程中，我们使用 SRTM1 数据，数据引用和有关使用该数据的政策可在土地过程分布式活动档案中心找到。

import os
import elevation
import rasterio as rio
import rioxarray as riox
from rasterio.plot import show

dem_path = '/DEM_Data/image_original.tif'
output = os.getcwd() + dem_path

#Extract dem based on bounds extracted from the ADM boundary polygon

bounds_combined = combined_polygon.bounds
west_c, south_c, east_c, north_c = bounds_combined
elevation.clip(bounds=bounds_combined, output=output, product='SRTM1')
dem = rio.open(output)
show(dem)

# The bounds does not extract in the same shape as the polygon but using the bounds
# The following script clips based on the polygon

raster = riox.open_rasterio(output)
# Shapely Polygon  to clip raster
geom = combined_polygon

# Use shapely polygon in clip method of rioxarray object to clip raster
clipped_raster = raster.rio.clip([geom])
 
dem_path1 = '/DEM_Data/image_clipped.tif'
output1 = os.getcwd() + dem_path1

# Save clipped raster
clipped_raster.rio.to_raster(output1)

dem = rio.open(output1)
show(dem)

请注意，我们需要提供绝对文件路径才能将 DEM 下载为 TIF 文件；否则将导致错误。

使用 Python 中的Elevation包提取 DEM ，并裁剪感兴趣的多边形（图片由作者提供）

使用 matplotlib 和 RichDEM 可视化 DEM 地形属性

在本节中，我们将使用matplotlib和richdem可视化数字高程模型的一些地形属性。

DEM 中的等高线是连接等高点的线。这些等高线有助于显示地形的三维形状，并可用于测量距离、计算分水岭和确定坡度。matplotlib 提供将等高线绘制为线或填充区域的选项，如以下代码所示。

import matplotlib.pyplot as plt

# Plot out data with Matplotlib's 'contour' 
fig = plt.figure(figsize = (12, 8))
ax = fig.add_subplot(111)
plt.contour(dem_array, cmap = "viridis", 
            levels = list(range(0, 5000, 100)))
plt.title("Elevation Contours of "+selected_province+" province")
cbar = plt.colorbar()
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

# For filled Contour plots, use contourf 
fig = plt.figure(figsize = (12, 8))
ax = fig.add_subplot(111)
plt.contourf(dem_array, cmap = "viridis", 
            levels = list(range(0, 5000, 500)))
plt.title("Elevation Contours of "+selected_province+" province")
cbar = plt.colorbar()
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

使用 matplotlib 可视化海拔轮廓（图片来自作者）

DEM 的直方图以条形图的形式呈现，每条柱状图代表特定高程值的频率。每条柱状图的高度对应于该特定高程值在数据集中出现的次数。DEM 直方图可用于识别高程极值，以及检测高程异常值或确定数据集内的整体高程范围。

dem_im = riox.open_rasterio(output1,masked=True)

f, ax = plt.subplots(figsize=(10, 6))
dem_im.plot.hist(ax=ax,
       color="purple")
ax.set(title="Distribution of DEM Elevation Values",
       xlabel='Elevation (meters)',
       ylabel='Frequency')
plt.show()

DEM 高程值分布（图片由作者提供）

坡度是一种地形属性或测量值，用于描述给定区域内高程的变化率。它衡量的是地表的陡度或倾斜度，通常以百分比或角度表示。

DEM 的坡度通常以度为单位，范围从 0°（平坦）到 90°（垂直）。亚美尼亚总体上是一个多山的国家，大部分地形为陡坡。该国也有一些平坦或缓坡的地区，主要分布在阿拉斯河和阿拉克斯河的河谷。

import richdem as rd

dem_slope = rd.TerrainAttribute(beau, attrib='slope_degrees')
rd.rdShow(dem_slope, axes=True, cmap='YlOrBr', figsize=(16, 10));

地形属性坡度（以度为单位）（图片由作者提供）

当坡度以百分比表示时，它是海拔（上升）与给定距离（运行）的比率，称为坡度上升-运行。14 % 的坡度上升-运行意味着水平方向每上升 100 个单位，坡度就会垂直上升 14 个单位。例如，最高的坡度上升-运行位于亚美尼亚的高加索山脉，靠近格鲁吉亚和阿塞拜疆的边界。亚美尼亚的高加索山脉是欧亚山脉系的一部分，位于该国南部。该山脉的平均海拔约为 4,500 米（14,764 英尺）。亚美尼亚高加索山脉的平均坡度为 15%。

import richdem as rd

beau  = rd.rdarray(dem_array, no_data=-9999)
slope = rd.TerrainAttribute(beau, attrib='slope_riserun')
rd.rdShow(slope, axes=True, cmap='jet', figsize=(16,10))

地形属性坡度上升-平移百分比（作者提供的图片）

坡向（以度为单位）是地球表面特定位置的坡度方向的数值测量单位。坡向测量范围为 0 至 360 度，其中 0° 代表北，90° 代表东，180° 代表南，270° 代表西。

aspect = rd.TerrainAttribute(beau, attrib='aspect')
rd.rdShow(aspect, axes=False, cmap='jet', figsize=(8,5.5))

地形属性坡向（以度为单位）（图片由作者提供）

结论

数字高程模型的可视化和分析对于气候变化、灾害风险缓解和基础设施规划等诸多用例至关重要，因为它能够分析地形特征，而这些特征对于理解气候变化动态、评估自然灾害风险以及设计有效的基础设施项目至关重要。DEM 数据也是生成数字地形模型的基础，可用于绘制易受洪水、山体滑坡和其他极端天气事件影响的区域。此外，DEM 数据还可用于创建详细的景观和地形 3D 模型，用于模拟潜在的基础设施项目（例如道路和桥梁），并评估其环境影响。

我们希望您现在能够更好地理解如何使用 Python 下载、探索和可视化数字高程模型 (DEM) 及其属性。掌握这些知识后，您可以使用 Python 创建详细的地形地图并分析不同区域的高程。此外，您还可以使用 Python 比较不同的 DEM 及其属性，深入了解地形随时间的变化。借助这些工具，您将能够深入探索和分析景观。

最后，让我们以亚美尼亚阿拉加茨索滕省的DEM可视化效果图来结束这篇博文。阿拉加茨山是亚美尼亚最高峰阿拉加茨山的所在地。阿拉加茨山是一座死火山，由四座独立的山峰组成，每座山峰都可以攀登，因此是徒步旅行者和登山者的梦想之地。