Python自动化处理断面高程数据

原创已于 2025-07-08 13:28:30 修改 · 991 阅读

CC 4.0 BY-SA版权

文章标签：

于 2025-07-07 15:07:51 首次发布

本文摘要：提出一个基于Python的自动化流程，用于沿断面线每1米间隔生成断面点并提取高程数据。首先读取断面线矢量数据（shp）和DEM栅格数据（tif），进行坐标转换和几何校验。核心算法包含：1）沿断面线按整数长度采样；2）坐标投影转换；3）DEM高程值提取。程序输出包含坐标、里程和高程的shp文件及结构化Excel表格。关键处理包括：线要素方向标准化、坐标边界检查、异常值处理等。该方案实现了从空间数据到断面特征的自动化提取，适用于水文分析、地形测量等领域，代码具备完整错误处理机制。

一、ai交流过程：现在给你一个断面线数据，“95-102E-断面线.shp”，如图1所示。

属性如图2所示。

再给你一个“105E-0.2m-dem.tif”，如图3所示。

现在想让你沿着断面线从左向右，每1m间隔生成一个断面点（注：左端点为起点，如果该断面线的长度不是整数，假设是n+a（其中n是断面长度整数部分，a是小数部分）那么最后一小段（即小数部分a）不要了。即终点为长度为n位置的点）。然后每条断面线都进行这样的操作。生成所有的断面点后，需要提取所有点位的x，y，里程值，高程值（值提取至点操作）。生成的excel文件格式如图4所示（所有生成的点需在一个excel文件里面）：

能用python实现这个过程吗？

二、最终代码：

import geopandas as gpd
import rasterio
from shapely.geometry import Point, LineString
import numpy as np
import os
from pyproj import Transformer
import pandas as pd

# 1. 读取断面线数据
line_shp_path = r"D:\0617test\0620\96-102E-唐家村断面线\96-102E-断面线.shp"

# 验证文件路径是否存在
if not os.path.exists(line_shp_path):
    raise FileNotFoundError(f"错误: 未找到断面线文件 '{line_shp_path}'")

# 读取断面线数据并添加异常处理
try:
    line_gdf = gpd.read_file(line_shp_path)
    print(f"成功读取断面线数据: {len(line_gdf)} 条记录")
    print(f"数据坐标系: {line_gdf.crs}")

    # 检查是否有几何列
    if 'geometry' not in line_gdf.columns:
        raise ValueError("错误: 数据中不包含几何列 'geometry'")

    # 验证是否为线要素
    geom_types = line_gdf.geometry.geom_type.unique()
    if not all(geom in ['LineString', 'MultiLineString'] for geom in geom_types):
        print(f"警告: 数据中包含非线要素类型: {geom_types}")

except Exception as e:
    print(f"读取断面线数据失败: {e}")
    raise

# 2. 读取 DEM 数据
dem_path = r"D:\0617test\0620\唐家村.tif"

# 验证 DEM 文件路径
if not os.path.exists(dem_path):
    raise FileNotFoundError(f"错误: 未找到 DEM 文件 '{dem_path}'")

try:
    with rasterio.open(dem_path) as dem_ds:
        dem_transform = dem_ds.transform
        dem_arr = dem_ds.read(1)
        dem_crs = dem_ds.crs
        print(f"成功读取 DEM 数据，尺寸: {dem_arr.shape}")
        print(f"DEM 坐标系: {dem_crs}")

except Exception as e:
    print(f"读取 DEM 数据失败: {e}")
    raise

# 创建坐标转换器：从断面线坐标系(EPSG:4543)转换到DEM坐标系
transformer = Transformer.from_crs("EPSG:4543", dem_crs, always_xy=True)


# 3. 定义函数：沿断面线生成指定间隔的点并提取高程
def generate_section_points(line_geom, interval=1):
    """
    沿断面线生成指定间隔的点，提取坐标、里程和高程
    :param line_geom: 断面线的 shapely LineString 几何对象（已确保左→右方向）
    :param interval: 采样间隔，单位米（需保证线要素坐标是投影坐标，单位为米）
    :return: 包含点几何、里程、x、y、高程的列表
    """
    # 计算线长度
    line_length = line_geom.length
    # 确定最大整数里程（去掉小数部分，如长度 10.9 → 取 10）
    max_integer_length = int(line_length)
    # 生成采样里程列表（从 0 到 max_integer_length，步长 interval）
    mileages = np.arange(0, max_integer_length+1, interval)

    points_data = []
    for mile in mileages:
        # 在断面上按里程插值点（从左端点开始）
        point = line_geom.interpolate(mile)
        x, y = point.x, point.y

        # 将点坐标从断面线坐标系(EPSG:4543)转换到DEM坐标系
        x_proj, y_proj = transformer.transform(x, y)

        # 从 DEM 提取高程
        # 转换为行列号（使用投影后的坐标）
        row, col = rasterio.transform.rowcol(dem_transform, x_proj, y_proj)

        # 检查行列号是否在DEM范围内
        if 0 <= row < dem_arr.shape[0] and 0 <= col < dem_arr.shape[1]:
            elevation = dem_arr[row, col]
        else:
            elevation = None  # 点在 DEM 范围外时设为 None

        points_data.append({
            "geometry": Point(x_proj, y_proj),  # 使用投影后的坐标创建点
            "section_na": "",  # 用于存储 DMName 字段值
            "section_id": "",  # 用于存储 Id 字段值
            "POINT_X": x_proj,  # 存储投影后的X坐标
            "POINT_Y": y_proj,  # 存储投影后的Y坐标
            "distance": mile,
            "elevation": elevation
        })
    return points_data


# 4. 遍历所有断面线，统一方向（左→右）后生成断面点，并赋值对应字段
all_points = []
for idx, row in line_gdf.iterrows():
    line_geom = row["geometry"]

    # 获取断面线的 Id 和 DMName 字段值
    line_id = row.get("Id")
    line_dm_name = row.get("DMName")

    # 处理 LineString：反转方向，确保左→右
    if isinstance(line_geom, LineString):
        # 反转坐标顺序，让起点变为左端点
        # reversed_coords = list(reversed(line_geom.coords))
        reversed_coords = line_geom.coords
        reversed_line = LineString(reversed_coords)
        section_points = generate_section_points(reversed_line, interval=1)

    # 处理 MultiLineString：逐个 LineString 反转方向
    elif isinstance(line_geom, MultiLineString):
        section_points = []
        for line_part in line_geom.geoms:
            if isinstance(line_part, LineString):
                reversed_coords = list(reversed(line_part.coords))
                reversed_line = LineString(reversed_coords)
                section_points.extend(generate_section_points(reversed_line, interval=1))

    else:
        print(f"警告: 第 {idx} 行不是线要素，跳过")
        continue

    # 根据线数据属性为断面名称（section_na）、断面序号（section_id）赋值
    for p in section_points:
        p["section_na"] = line_dm_name
        p["section_id"] = line_id
    all_points.extend(section_points)

print(f"共生成 {len(all_points)} 个断面点")

# 5. 将所有点数据转换为 GeoDataFrame
if not all_points:
    raise ValueError("错误: 未生成任何断面点，请检查断面线数据")

# 使用DEM的坐标系
points_gdf = gpd.GeoDataFrame(all_points, crs=dem_crs)

# 6. 保存为 shp 文件
output_shp_path = "D:/0617test/0620/section_points.shp"
try:
    points_gdf.to_file(output_shp_path)
    print(f"成功保存断面点数据到: {output_shp_path}")
except Exception as e:
    print(f"保存断面点数据失败: {e}")
    raise
shp_path = output_shp_path
# 替换为实际生成的 shp 文件路径
# 替换为要保存的 excel 文件路径
excel_output_path = "D:/0617test/0620/唐家村.xlsx"

# 读取 shp 文件
gdf = gpd.read_file(shp_path)

# 构造目标 DataFrame 的字段列表
target_columns = [
    "断面名称", "断面序号", "POINT_X", "POINT_Y", "距离", "高程",
    "点类型标记", "断面类型", "left_roughness", "river_roughness",
    "right_roughness", "sigma_s", "slope"
]

# 准备要写入 Excel 的数据列表
data = []
for idx, row in gdf.iterrows():
    row_data = {
        "断面名称": row["section_na"],
        "断面序号": row["section_id"],
        "POINT_X": row["POINT_X"],
        "POINT_Y": row["POINT_Y"],
        "距离": row["distance"],
        "高程": row["elevation"],
        # 其他字段先填充为空字符串，可根据实际需求调整
        "点类型标记": "",
        "断面类型": "",
        "left_roughness": "",
        "river_roughness": "",
        "right_roughness": "",
        "sigma_s": "",
        "slope": ""
    }
    data.append(row_data)

# 创建 DataFrame
result_df = pd.DataFrame(data, columns=target_columns)

# 保存为 Excel 文件
result_df.to_excel(excel_output_path, index=False)
print(f"成功生成 Excel 文件：{excel_output_path}")