LXX
读完需要
速读仅需 1 分钟
前言：
最近一直加班出差拖更太久.......最近看到小伙伴询问如何使用 Python 绘制地形剖面图，其实在前面有一期就分享过在 GEE 中提取高程点并绘制剖面图的方法，使用 Python 绘制的原理其实也差不多，下面是具体的实现方法。
希望各位同学点个关注，点个小赞，这将是更新的动力，不胜感激❥(^_-)
1
绘制地形剖面图（也称为纵断面图）是将某一路径上的高程变化显示为二维图形的过程。通常用于显示沿路径的地形起伏情况，如山脉、谷地等。
1.1
具体步骤
1.读取 DEM（数字高程模型）数据：从 DEM 文件中读取高程数据。这些数据通常存储为栅格（像素）格式，每个像素代表一个地理位置的高程值。
2.选择路径：确定剖面路径的起点和终点。这条路径可以是直线，也可以是曲线。
3.提取路径上的高程数据：沿路径从 DEM 数据中提取对应位置的高程值。
4.计算距离：计算路径上各点之间的距离，以便绘制剖面图的横轴。
5.绘制剖面图：将路径上的距离和对应的高程值绘制成图。
1.2
路径图绘制
在绘制了地形剖面图后，我们也想看一下路径在地图中的位置信息，因此采用地图神器 folium 包来生成。
Folium 是建立在 Python 生态系统的数据整理 Datawrangling 能力和 Leaflet.js 库的映射能力之上的开源库。用 Python 处理数据，然后用 Folium 将它在 Leaflet 地图上进行可视化。Folium 能够将通过 Python 处理后的数据轻松地在交互式的 Leaflet 地图上进行可视化展示。它不单单可以在地图上展示数据的分布图，还可以使用 Vincent/Vega 在地图上加以标记。
这个开源库中有许多来自 OpenStreetMap、MapQuest Open、MapQuestOpen Aerial、Mapbox 和 Stamen 的内建地图元件，而且支持使用 Mapbox 或 Cloudmade 的 API 密钥来定制个性化的地图元件。Folium 支持 GeoJSON 和 TopoJSON 两种文件格式的叠加，也可以将数据连接到这两种文件格式的叠加层，最后可使用 color-brewer 配色方案创建分布图。
Folium 可以让你用 Python 强大生态系统来处理数据，然后用 Leaflet 地图来展示。Folium 内置一些来自 OpenStreetMap、MapQuest Open、MapQuest Open Aerial、Mapbox 和 Stamen 的地图元件(tilesets)，并且支持用 Mapbox 或者 Cloudmade API keys 来自定义地图元件。Folium 支持 GeoJSON 和 TopJSON 叠加(overlays)，绑定数据来创造一个分级统计图(Choropleth map)。但是，Folium 库绘制热点图的时候，需要联网才可显示。
# 可添加不同的底图folium.TileLayer('openstreetmap', name='OpenStreetMap').add_to(map_)folium.TileLayer('Stamen Terrain', name='Stamen Terrain', attribution="Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.").add_to(map_)folium.TileLayer('Stamen Toner', name='Stamen Toner', attribution="Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.").add_to(map_)folium.TileLayer('Stamen Watercolor', name='Stamen Watercolor', attribution="Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.").add_to(map_)
location：经纬度，list 或者 tuple 格式，顺序为 latitude, longitude zoom_start：缩放值，默认为 10，值越大比例尺越小，地图放大级别越大 control_scale：Bool型，控制是否在地图上添加比例尺，默认为 False 即不添加 tiles：显示样式，默认 “OpenStreetMap”，也就是开启街道显示 crs：地理坐标参考系统，默认为 “EPSG3857”
2
1. 读取DEM数据
我们使用rasterio库来读取DEM文件，DEM文件通常是GeoTIFF格式，包含栅格数据和地理空间信息。
import rasteriodef read_dem(dem_file):with rasterio.open(dem_file) as src:data = src.read(1) # 读取第一波段（通常是唯一波段），即高程数据transform = src.transform # 获取地理变换参数，用于像素坐标到地理坐标的转换return data, transform
2. 选择路径
选择路径上的起点和终点。路径默认是直线。
start_point = (122.00159416410716, 40.722738546800684) # 起点 (经度, 纬度)end_point = (121.97962150785716, 41.31952738477779) # 终点 (经度, 纬度)
3. 提取路径上的高程数据
沿路径从DEM数据中提取高程值。我们可以使用线性插值的方法，沿着路径以均匀间隔取样，然后从DEM中提取每个点的高程。
import numpy as npdef get_profile(dem_data, transform, start_point, end_point, num_points=100):# 计算起点和终点在DEM栅格上的行列号（像素坐标）start_row, start_col = ~transform * start_pointend_row, end_col = ~transform * end_point# 生成路径上num_points个点的行列坐标row_coords = np.linspace(start_row, end_row, num_points)col_coords = np.linspace(start_col, end_col, num_points)# 提取每个点的高程值elevations = []for row, col in zip(row_coords, col_coords):row = int(np.round(row))col = int(np.round(col))if 0 <= row < dem_data.shape[0] and 0 <= col < dem_data.shape[1]:elevations.append(dem_data[row, col])else:elevations.append(np.nan) # 超出范围的地方用NaN表示return elevations, row_coords, col_coords
4. 计算路径上各点之间的距离
使用geopy库计算路径上每个采样点之间的距离，以便作为剖面图的横轴。
from geopy.distance import geodesicdef calculate_distance(lat_lon_coords):distances = [0] # 起始距离为0for i in range(1, len(lat_lon_coords)):dist = geodesic(lat_lon_coords[i - 1], lat_lon_coords[i]).kilometersdistances.append(distances[-1] + dist)return distances
5. 绘制剖面图
使用matplotlib将路径上的距离和对应的高程绘制成图。
import matplotlib.pyplot as pltdef plot_profile(distances, elevations, output_file):plt.figure(figsize=(10, 6))plt.fill_between(distances, elevations, color='lightgreen', alpha=0.7)plt.plot(distances, elevations, color='green', lw=2)plt.xlabel('Distance (km)')plt.ylabel('Elevation (m)')plt.title('Terrain Profile')plt.grid(True)plt.savefig(output_file, format='png', dpi=300)plt.show()
3
高程剖面图如下：
路径图如下：（代码运行后会生成一个HTML链接，打开即可）
【完整代码：】
import rasterioimport numpy as npimport matplotlib.pyplot as pltfrom geopy.distance import geodesicfrom rasterio.transform import xyimport folium# 获取地形剖面数据def get_profile(dem_file, start_point, end_point, num_points=100):with rasterio.open(dem_file) as src:transform = src.transformdata = src.read(1)start_row, start_col = ~transform * start_pointend_row, end_col = ~transform * end_pointrow_coords = np.linspace(start_row, end_row, num_points)col_coords = np.linspace(start_col, end_col, num_points)elevations = []for row, col in zip(row_coords, col_coords):row = int(np.round(row))col = int(np.round(col))if 0 <= row < data.shape[0] and 0 <= col < data.shape[1]:elevations.append(data[row, col])else:elevations.append(np.nan)return elevations, row_coords, col_coords, transform# 计算距离def calculate_distance(lat_lon_coords):distances = [0] # 起始距离为0for i in range(1, len(lat_lon_coords)):dist = geodesic(lat_lon_coords[i - 1], lat_lon_coords[i]).kilometersdistances.append(distances[-1] + dist)return distances# 绘制地形剖面图def plot_profile(elevations, distances, output_file):plt.figure(figsize=(10, 6)) # 调整图的宽高比例plt.fill_between(distances, elevations, color='lightgreen', alpha=0.7) # 填充下方区域plt.plot(distances, elevations, color='green', lw=1) # 绘制剖面线# 标注最高点max_elev = max(elevations)max_idx = elevations.index(max_elev)plt.scatter([distances[max_idx]], [max_elev], color='purple', zorder=5)plt.text(distances[max_idx], max_elev, f'{int(max_elev)}m', fontsize=10, ha='center')# 标注起点和终点plt.scatter([distances[0], distances[-1]], [elevations[0], elevations[-1]], color='purple', zorder=5)plt.text(distances[0], elevations[0], f'{int(elevations[0])}m', fontsize=10, ha='center')plt.text(distances[-1], elevations[-1], f'{int(elevations[-1])}m', fontsize=10, ha='center')# 添加坐标轴标签plt.rcParams['font.sans-serif'] = ['STSong']plt.rcParams['axes.unicode_minus'] = Falseplt.xlabel('Distance (km)')plt.ylabel('Elevation (m)')plt.title('高程剖面图')# 调整坐标轴范围plt.xlim(0, distances[-1])plt.ylim(min(elevations) - 10, max(elevations) + 10)# 显示网格线plt.grid(True)# 保存图片plt.savefig(output_file, format='png', dpi=300)plt.close()# 绘制路径地图def create_map(start_point, end_point, map_output_file):# 创建一个以起点为中心的地图map_center = [(start_point[1] + end_point[1]) / 2, (start_point[0] + end_point[0]) / 2]map_ = folium.Map(location=map_center, zoom_start=10)# 添加不同的底图folium.TileLayer('openstreetmap', name='OpenStreetMap').add_to(map_)folium.TileLayer(tiles='https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}',name='Esri World Imagery',attr='Tiles © Esri — Source: Esri, i-cubed, USDA, USGS, AEX, GeoEye, Getmapping, Aerogrid, IGN, IGP, UPR-EGP, and the GIS User Community').add_to(map_)folium.TileLayer(tiles='https://stamen-tiles.a.ssl.fastly.net/toner/{z}/{x}/{y}.png',name='Stamen Toner',attr='Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.').add_to(map_)folium.TileLayer(tiles='https://stamen-tiles.a.ssl.fastly.net/watercolor/{z}/{x}/{y}.jpg',name='Stamen Watercolor',attr='Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.').add_to(map_)# 添加起点和终点的标记folium.Marker(location=[start_point[1], start_point[0]], popup="Start Point", icon=folium.Icon(color='red')).add_to(map_)folium.Marker(location=[end_point[1], end_point[0]], popup="End Point", icon=folium.Icon(color='red')).add_to(map_)# 添加起点到终点的路径线folium.PolyLine(locations=[(start_point[1], start_point[0]), (end_point[1], end_point[0])], color="red", weight=2.5).add_to(map_)# 添加图层控制器folium.LayerControl().add_to(map_)# 保存地图到HTML文件map_.save(map_output_file)# 设置参数dem_file = 'C:/Users/LXX/Downloads/panjin.tif'start_point = (122.00159416410716, 40.722738546800684) # 起点坐标 (x, y)end_point = (121.97962150785716, 41.31952738477779) # 终点坐标 (x, y)profile_output_file = 'profile.png'map_output_file = 'map.html'# 生成地形剖面图elevations, row_coords, col_coords, transform = get_profile(dem_file, start_point, end_point)lat_lon_coords = [(y, x) for x, y in [xy(transform, row, col, offset='center') for row, col in zip(row_coords, col_coords)]]distances = calculate_distance(lat_lon_coords)plot_profile(elevations, distances, profile_output_file)# 生成路径地图create_map(start_point, end_point, map_output_file)print(f'Profile plot saved as {profile_output_file}')print(f'Map saved as {map_output_file}')
