3. 三维地形建模：从DEM到可视化
好，咱们进入第三章。这一章，我打算聊聊三维地形建模。
你想想看，一个风场数字孪生系统，如果地形是平的，那还有什么意义？风是跟着地形跑的——山脊加速、山谷狭管效应、背风坡湍流……这些都得靠真实的地形数据才能模拟出来。
我个人习惯把地形建模拆成四步：拿数据 → 生成网格 → 坐标转换 → 可视化。每一步都有坑，咱们一个一个说。
核心流程概览
3.1 DEM数据获取：别拿错数据
DEM，Digital Elevation Model，数字高程模型。说白了就是一张灰度图，每个像素存的是海拔高度。
目前主流的免费数据源有三个：
| 数据源 | 分辨率 | 覆盖范围 | 我的评价 |
|---|---|---|---|
| SRTM | 30m / 90m | 全球（±60°纬度） | 最稳定，我首选 |
| ASTER GDEM | 30m | 全球 | 细节多，但空洞也多 |
| ALOS World 3D | 30m | 全球 | 精度不错，下载麻烦 |
我在项目中遇到过一个问题：用ASTER GDEM做某山地风场，结果山脊线附近全是空洞，补都补不回来。后来老老实实换回SRTM，虽然分辨率低一点，但数据完整。
我的建议：如果风场范围不大（比如10km×10km以内），优先用SRTM 30m数据。如果范围大，90m也够用。别一味追求高分辨率，数据完整性更重要。
获取方式也很简单：
- USGS EarthExplorer：最全，但操作略繁琐
- OpenTopography：界面友好，支持直接裁剪
- 地理空间数据云：国内用户首选，速度快
嗯，这里要注意：下载时一定要确认坐标系。大部分DEM数据默认是WGS84地理坐标系（经纬度），但有些会给你投影坐标系。我建议统一用WGS84，后面再转。
3.2 地形网格生成：从点云到Mesh
拿到DEM数据后，其实就是一堆像素点——每个点有(x, y, z)三个值。x和y是经纬度，z是海拔。
但数字孪生需要的是三角形网格，也就是Mesh。怎么把像素点变成三角形？
两种主流方法：
- 规则网格法：直接把DEM像素点连成四边形，再拆成两个三角形。简单粗暴，适合平坦地形。
- Delaunay三角剖分：根据地形复杂度自适应生成三角形。陡峭区域三角形密，平坦区域三角形疏。效果好，但计算量大。
我个人习惯用Delaunay。为什么？你想想看，风场地形往往有山脊、山谷这些特征，规则网格在这些地方会浪费大量三角形，而Delaunay能自动优化。
代码示例（Python + scipy）：
import numpy as np
from scipy.spatial import Delaunay
import matplotlib.pyplot as plt
# 假设dem_data是二维数组，shape=(rows, cols)
rows, cols = dem_data.shape
x = np.arange(cols)
y = np.arange(rows)
xx, yy = np.meshgrid(x, y)
# 展平为点云
points = np.c_[xx.ravel(), yy.ravel(), dem_data.ravel()]
# Delaunay三角剖分（只基于x,y）
tri = Delaunay(points[:, :2])
# 输出三角形索引
triangles = tri.simplices
print(f"生成了 {len(triangles)} 个三角形")注意：Delaunay剖分时只用了x和y坐标，z值只作为属性。千万别把z也当成坐标输入，否则地形会被压扁或拉伸。
我曾经犯过这个错——把(x, y, z)全扔进Delaunay，结果生成的地形像被揉过的纸团。排查了半天才发现是坐标维度搞错了。
3.3 GIS坐标转换：别让风跑偏
这一步，很多人会忽略。但我要说：坐标转换做不好，后面的仿真全白搭。
DEM数据是WGS84经纬度坐标，但风场仿真通常需要投影坐标系（比如UTM）。为什么？因为仿真软件要算距离、算面积，经纬度不是等距的——在赤道1度经度≈111km，到北极就缩成0了。
转换流程：
- 确定投影带：根据风场中心经度，选择对应的UTM带。中国主要跨UTM 43-53带。
- 执行转换：用pyproj或GDAL，把每个点的经纬度转成米为单位的东西向和南北向坐标。
- 验证：转换后检查一下，两个点之间的距离是否合理。
代码示例：
from pyproj import Transformer
# 定义坐标系：WGS84 → UTM 50N（北京附近）
transformer = Transformer.from_crs("EPSG:4326", "EPSG:32650")
# 假设有一个点：东经116.4°，北纬39.9°
lon, lat = 116.4, 39.9
x_m, y_m = transformer.transform(lat, lon)
print(f"转换后坐标：东向 {x_m:.2f} m，北向 {y_m:.2f} m")
# 输出：东向 452345.67 m，北向 4417890.12 m一个小技巧：转换后记得把坐标原点平移到(0,0)，方便后续处理。比如所有x减去最小值，所有y减去最小值。这样地形网格就在第一象限了，看着舒服。
3.4 地形可视化：让数据看得见
最后一步，把Mesh渲染出来。我常用的工具有三个：
| 工具 | 适用场景 | 优点 | 缺点 |
|---|---|---|---|
| Three.js | Web端展示 | 轻量、跨平台 | 大规模数据卡顿 |
| Unity | 高保真仿真 | 效果好、交互强 | 开发成本高 |
| Matplotlib | 快速验证 | 代码少、出图快 | 不能交互 |
我个人习惯：开发阶段用Matplotlib快速看效果，交付阶段用Three.js或Unity。
Matplotlib画地形图，几行代码搞定：
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# 假设x_mesh, y_mesh, z_mesh是转换后的网格数据
ax.plot_surface(x_mesh, y_mesh, z_mesh,
cmap='terrain', alpha=0.8)
ax.set_xlabel('东向 (m)')
ax.set_ylabel('北向 (m)')
ax.set_zlabel('海拔 (m)')
plt.show()嗯，这里要注意：颜色映射很重要。用'tetrain'色图，低海拔绿色、中海拔棕色、高海拔白色，一眼就能看出地形起伏。别用'jet'或'hot'，那些是给温度场用的。
避坑指南：
- 数据空洞：DEM数据偶尔会有空洞（NoData值）。我曾经遇到一个风场，山脊上有个大洞，仿真出来的风速直接爆表。解决方案：用周围像素插值填补，或者换数据源。
- 分辨率不匹配：如果DEM分辨率是30m，但风场网格是100m，直接降采样会导致地形细节丢失。我建议先做平滑滤波，再降采样。
- 坐标系混淆：这是最常见的坑。我见过有人把经纬度当成米来算，结果地形被压成了一张纸。每次转换后，务必打印几个点的坐标，肉眼检查一下。
好了，三维地形建模就聊到这儿。记住四个关键词：DEM、网格、坐标、可视化。每一步都不难，但每一步都有细节。做风场数字孪生，地形是地基——地基没打好，上面的仿真再漂亮也是空中楼阁。