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

背景与痛点

OSGB（OpenSceneGraph Binary）是三维 GIS 中广泛使用的模型格式，其采用分块 LOD 结构存储高精度地形和建筑数据。但在实际项目中，我们常需从中提取传统 2D 等高线用于规划分析，这面临两大挑战：

1. 数据转换复杂性：OSGB 的层级结构不同于常规 DEM 数据，直接提取需处理坐标系转换和精度损失问题
2. 计算效率瓶颈：城市级模型常包含千万级三角面片，传统算法易内存溢出

方法对比

1. 高程点插值法

原理：
– 从 OSGB 中采样离散高程点
– 通过反距离权重（IDW）或克里金插值生成规则网格 DEM
– 基于 DEM 生成等高线

优势场景：
– 数据存在局部缺失时（如建筑遮挡区域）
– 需要平滑的等高线输出（适合小比例尺制图）

2. 三角网剖分法(TIN)

原理：
– 直接读取 OSGB 的三角面片顶点
– 构建不规则三角网(TIN)
– 在 TIN 上追踪等高线

独特价值：
– 保留原始地形细节（适合大比例尺工程应用）
– 天然避免插值导致的精度损失

核心实现

方法一：高程点插值法代码示例

import numpy as np
from osgeo import gdal, ogr

def generate_contours_from_osgb(osgb_path, output_shp, interval=5):
    """
    参数说明：osgb_path: 输入 OSGB 模型路径
    output_shp: 输出等高线 Shapefile 路径
    interval: 等高距（米）"""
    # Step1: 使用 Cass 3D API 提取高程点
    points = cass3d.extract_vertices(osgb_path)  # 返回 Nx3 的 numpy 数组

    # Step2: 创建规则网格 DEM
    xmin, xmax = points[:,0].min(), points[:,0].max()
    ymin, ymax = points[:,1].min(), points[:,1].max()
    grid_res = 1.0  # 网格分辨率

    dem = idw_interpolation(points, 
                          bbox=(xmin, ymin, xmax, ymax),
                          resolution=grid_res)

    # Step3: 生成等高线
    driver = gdal.GetDriverByName('ESRI Shapefile')
    ds = driver.Create(output_shp, 0, 0, 0, gdal.GDT_Unknown)
    layer = ds.CreateLayer('contours', geom_type=ogr.wkbLineString)

    gdal.ContourGenerate(dem.GetRasterBand(1), 
                        interval, 0, [], 0, 0, 
                        layer, -1, -1)

方法二：三角网剖分法关键代码

def tin_contouring(osgb_path, output_shp):
    """直接基于 TIN 的等高线追踪"""
    # 获取三角网数据
    vertices, faces = cass3d.extract_mesh(osgb_path)

    # 构建拓扑结构
    tin = TinModel(vertices, faces)

    # 等高线追踪算法
    contours = []
    for elevation in np.arange(z_min, z_max, interval):
        contours += tin.trace_contour(elevation)

    # 保存结果
    save_to_shapefile(contours, output_shp)

性能考量

通过测试某城市 30km² OSGB 模型（约 1200 万个三角面片）：

内存优化技巧：
– 对于插值法：使用分块处理策略，每块大小建议 50 万 -100 万点
– 对于 TIN 法：启用 Cass 3D 的 LOD 简化功能，可减少 30%-50% 面片数

避坑指南

常见拓扑错误修复

1. 等高线自相交：
2. 在 GDAL 生成后执行 ogr.Geometry.Buffer(0) 操作

3. 或使用 shapely.ops.unary_union 进行几何修复

4. 高程异常值：

5. 预处理时过滤 Z 值超出合理范围的点

   valid_points = points[(points[:,2] > 50) & (points[:,2] < 300)]  # 假设合理高程范围 50-300 米

总结与延伸

选型决策树：
1. 项目精度要求 >0.5 米 → 选择三角网剖分法
2. 存在大量数据缺失 → 选择高程点插值法
3. 硬件内存 <8GB → 优先考虑分块插值方案

建议读者使用 Cass 3D 自带的 SampleData/terrain.osgb 进行首次尝试，通过调整等高距参数观察输出效果差异。对于超大规模数据，可结合 Cass 3D 的分布式计算模块实现并行处理。