基于Cass高程点生成三维地形的完整指南：从数据准备到可视化实现

1次阅读

没有评论

共计 2559 个字符，预计需要花费 7 分钟才能阅读完成。

Cass 高程数据是测绘工程中常见的格式，通常以文本文件存储，包含点号、坐标 X /Y/ Z 值（平面坐标 + 高程）。这类数据的特点在于：

离散采样：通过全站仪或 RTK 测量获得，仅记录关键地形特征点（如山脊、谷底、坡度变化点）
非均匀分布：在平坦区域采样稀疏，地形复杂区域采样密集
工程标准：符合《工程测量规范》GB50026 要求，平面精度通常达厘米级

典型应用场景包括：

土方工程计算
洪水淹没分析
三维数字孪生基础
地质稳定性评估

将离散高程点转为连续地形表面需要空间插值，常见算法特性如下：

反距离权重法（IDW）
原理：距离越近的点权重越大
优点：计算简单，适合均匀分布数据
缺点：易产生 ” 牛眼效应 ”，不适合突变地形
克里金法（Kriging）
原理：基于地理统计学，考虑空间自相关性
优点：可生成误差估计曲面
缺点：需要计算变异函数，计算量大
自然邻域法（Natural Neighbor）
原理：基于 Voronoi 图动态调整权重
优点：适应不规则采样
缺点：边缘效果较差

建议选择策略：

平坦区域：IDW（计算效率高）
复杂地形：普通克里金（精度优先）
大数据量：线性插值（性能折衷）

import pandas as pd
from osgeo import gdal, osr

# 读取 Cass 格式文本数据
def parse_cass_txt(file_path):
    """
    解析 Cass 格式高程点文本文件
    格式示例：点号, 编码,Y 坐标,X 坐标,H 高程
    1,,500000.123,3000000.456,102.34
    """
    df = pd.read_csv(file_path, header=None, 
                    names=['id', 'code', 'y', 'x', 'z'],
                    dtype={'x':float, 'y':float, 'z':float})
    return df[['x', 'y', 'z']].values  # 返回 N * 3 数组

# 坐标转换（如需要）def transform_coords(points, src_epsg=4528, dst_epsg=4326):
    """将平面坐标转为经纬度"""
    transformer = osr.CoordinateTransformation(osr.SpatialReference().ImportFromEPSG(src_epsg),
        osr.SpatialReference().ImportFromEPSG(dst_epsg))
    return [transformer.TransformPoint(x,y,z)[:2] + (z,) 
            for x,y,z in points]

from scipy.interpolate import griddata
import numpy as np

def interpolate_surface(points, method='cubic', resolution=10):
    """
    执行插值计算
    :param points: (N,3)数组
    :param method: 'linear'/'cubic'/'nearest'
    :param resolution: 输出格网间距(米)
    :return: (grid_x, grid_y, grid_z) 插值结果
    """
    # 生成规则格网
    x_min, x_max = points[:,0].min(), points[:,0].max()
    y_min, y_max = points[:,1].min(), points[:,1].max()

    grid_x = np.arange(x_min, x_max, resolution)
    grid_y = np.arange(y_min, y_max, resolution)
    grid_x, grid_y = np.meshgrid(grid_x, grid_y)

    # 执行插值
    grid_z = griddata(points[:,:2], points[:,2], 
        (grid_x, grid_y), method=method)

    return grid_x, grid_y, grid_z

import plotly.graph_objects as go

def plot_3d_surface(x, y, z):
    """使用 Plotly 生成交互式三维地形"""
    fig = go.Figure(data=[go.Surface(z=z, x=x, y=y)])
    fig.update_layout(
        scene=dict(xaxis_title='东坐标(m)',
            yaxis_title='北坐标(m)',
            zaxis_title='高程(m)')
    )
    fig.show()

处理百万级高程点时建议：

数据分块
将研究区域划分为 500x500m 的网格
对每个网格单独处理后再拼接

并行计算

from multiprocessing import Pool

def parallel_interpolation(points, n_workers=4):
    """多进程分块插值"""
    chunks = np.array_split(points, n_workers)
    with Pool(n_workers) as p:
        results = p.map(lambda x: interpolate_surface(x, 'linear'), 
            chunks)
    return combine_results(results)  # 结果合并函数