如何在Python中沿Shapefile线提取Raster的数值剖面？

我刚好做过类似的需求，分享一套实用的Python实现方案，完全能对标QGIS里的Terrain Profile工具，还能保留完整的地理信息：

结合rasterio处理栅格地理信息，用shapely处理线要素，再通过栅格采样获取剖线上的数值，最后把采样点的地理坐标和对应栅格值关联起来——既保留了地理信息，又实现了剖面提取。

步骤1：加载并对齐数据 #

先读取GeoTIFF栅格和线Shapefile，确保两者的坐标系一致：

import rasterio
import geopandas as gpd

# 加载栅格（比如DEM地形数据）
with rasterio.open("dem.tif") as src:
    raster_data = src.read(1)
    raster_transform = src.transform
    raster_crs = src.crs

# 加载线要素（如果坐标系不一致，先转换）
line_gdf = gpd.read_file("profile_line.shp")
if line_gdf.crs != raster_crs:
    line_gdf = line_gdf.to_crs(raster_crs)

# 取目标剖面线（如果有多个线要素，按需调整索引）
target_line = line_gdf.geometry.iloc[0]

步骤2：生成剖线上的采样点 #

为了得到连续的剖面，在线上生成均匀分布的采样点，点数越多剖面越精细：

from shapely.geometry import Point

def generate_profile_samples(line, sample_count=500):
    sample_points = []
    line_length = line.length
    # 按距离均匀生成采样点
    for i in range(sample_count + 1):
        distance = i * line_length / sample_count
        sample_point = line.interpolate(distance)
        sample_points.append((sample_point.x, sample_point.y))
    return sample_points

# 生成500个采样点（可根据需求调整）
sample_coords = generate_profile_samples(target_line, sample_count=500)

步骤3：提取采样点的栅格值 #

用rasterio的sample方法直接提取地理坐标对应的栅格值，自动处理投影和坐标转换：

with rasterio.open("dem.tif") as src:
    # 提取每个采样点的栅格值（单波段取第一个元素）
    raster_values = [val[0] for val in src.sample(sample_coords)]

步骤4：关联地理信息并可视化 #

把采样点的坐标、距离、栅格值整合到DataFrame，还能画出和QGIS类似的剖面图表：

import pandas as pd
import matplotlib.pyplot as plt

# 构建包含完整地理信息的数据集
profile_df = pd.DataFrame({
    "x": [coord[0] for coord in sample_coords],
    "y": [coord[1] for coord in sample_coords],
    "distance_along_line": [i * target_line.length / len(sample_coords) for i in range(len(sample_coords))],
    "raster_value": raster_values
})

# 绘制地形剖面
plt.figure(figsize=(10, 4))
plt.plot(profile_df["distance_along_line"], profile_df["raster_value"], linewidth=2, color="#2ecc71")
plt.xlabel("Distance along profile (meters)")
plt.ylabel("Elevation (meters)")
plt.title("Terrain Profile")
plt.grid(True, alpha=0.3)
plt.show()

额外优化：处理NoData值 #

如果栅格存在无数据区域，可以过滤或填充这些值：

with rasterio.open("dem.tif") as src:
    nodata_value = src.nodata

# 过滤掉NoData的采样点
clean_profile_df = profile_df[profile_df["raster_value"] != nodata_value]

这套流程完全用Python实现，既能保留栅格的地理信息，又能实现和QGIS Terrain Profile一致的剖面提取效果，还能轻松集成到自动化工作流里。

内容的提问来源于stack exchange，提问作者Simon