地形图读取与可视化：从GeoTIFF到3D渲染

做风电项目，第一件事就是看地形。我见过太多团队，风数据算得天花乱坠，结果地形数据读错了，整个项目白干。今天咱们就聊聊，怎么用Python把地形图读进来、画出来、看明白。

1. 为什么是GeoTIFF？

说白了，普通图片不带地理坐标。你拍张照片，像素(100,200)是啥经纬度？不知道。GeoTIFF不一样，每个像素都绑定了真实的地理位置。我在内蒙古做项目时，拿到的地形数据就是GeoTIFF格式，一个文件几百兆，但信息量巨大。

GeoTIFF里存的是啥？是数字高程模型（DEM）。每个像素的值，代表那个点的海拔高度。嗯，就这么简单。

2. 用Rasterio读地形数据

我个人习惯用Rasterio，它比GDAL的Python绑定好用多了。来看看怎么读：

import rasterio
import numpy as np

# 打开GeoTIFF文件
with rasterio.open('terrain.tif') as src:
    # 读取第一波段（DEM数据）
    elevation = src.read(1)
    # 获取地理变换参数
    transform = src.transform
    # 获取坐标系
    crs = src.crs
    
print(f"数据形状: {elevation.shape}")
print(f"坐标系: {crs}")
print(f"变换参数: {transform}")

你想想看，就这么几行代码，地形数据就进内存了。我在云南一个山地风场项目里，就是用这个方式读入了30米分辨率的DEM数据，整个场区的地形一目了然。

3. 绘制等高线

等高线是风电选址的基础。风机放在哪，得看等高线判断山脊走向。我以前犯过傻，直接用原始数据画等高线，结果图又密又乱。后来学乖了，先做平滑处理。

import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter

# 平滑地形数据
elevation_smooth = gaussian_filter(elevation, sigma=2)

# 创建图形
fig, ax = plt.subplots(figsize=(12, 8))

# 绘制等高线
contour = ax.contour(elevation_smooth, 
                     levels=20,  # 等高线条数
                     colors='brown',
                     linewidths=0.8)

# 标注等高线数值
ax.clabel(contour, inline=True, fontsize=8)

# 添加颜色填充
contour_fill = ax.contourf(elevation_smooth, 
                           levels=20, 
                           cmap='terrain')

plt.colorbar(contour_fill, label='海拔 (m)')
plt.title('风电场区等高线图')
plt.show()

为什么会用gaussian_filter？因为原始DEM数据可能有噪点，直接画等高线会看到很多小圈圈，那其实是数据噪声，不是真实地形。我建议sigma值设在1-3之间，具体看数据分辨率。

4. 3D地形渲染

2D等高线看多了，换个3D视角更直观。特别是给甲方汇报时，3D地形图一出来，效果拉满。我一般用Matplotlib的mplot3d工具包：

from mpl_toolkits.mplot3d import Axes3D

# 创建网格坐标
rows, cols = elevation.shape
x = np.arange(0, cols, 1)
y = np.arange(0, rows, 1)
X, Y = np.meshgrid(x, y)

# 创建3D图形
fig = plt.figure(figsize=(14, 10))
ax = fig.add_subplot(111, projection='3d')

# 绘制地形表面
surf = ax.plot_surface(X, Y, elevation_smooth,
                       cmap='terrain',
                       rstride=4,  # 步长，控制渲染密度
                       cstride=4,
                       alpha=0.9)

# 设置视角
ax.view_init(elev=30, azim=135)

# 添加颜色条
fig.colorbar(surf, ax=ax, shrink=0.5, label='海拔 (m)')

ax.set_xlabel('X (像素)')
ax.set_ylabel('Y (像素)')
ax.set_zlabel('海拔 (m)')
plt.title('风电场区3D地形渲染')
plt.show()

这里有个坑：rstride和cstride参数。数据量大的时候，步长设1会卡死。我建议设4或5，既能看到地形细节，又不会让电脑崩溃。嗯，我在一台8GB内存的笔记本上跑过10000x10000的DEM，步长设8才勉强跑动。

5. 知识体系梳理

来，我画个图帮你理清思路：

6. 避坑指南

做地形图处理这些年，我踩过的坑不少。挑几个典型的说说：

• 坐标系不统一：我曾经拿到一个GeoTIFF，坐标系是WGS84，但另一个数据是北京54。直接叠加，位置差了500米。一定要检查src.crs，用rasterio.warp.reproject统一坐标系。
• 数据范围异常：有次读DEM，发现最高点海拔9999米。一看，原来是空值用9999填充了。用np.where(elevation > 9000, np.nan, elevation)处理一下就好。
• 内存爆炸：高分辨率DEM动不动就上GB。我建议用rasterio.windows分块读取，或者先降采样再处理。

7. 完整代码示例

最后，给你一个完整的处理流程。这是我做项目时常用的模板：

import rasterio
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter
from mpl_toolkits.mplot3d import Axes3D

def process_dem(filepath, smooth_sigma=2, contour_levels=20):
    """
    完整的地形图处理流程
    """
    # 1. 读取数据
    with rasterio.open(filepath) as src:
        elevation = src.read(1)
        transform = src.transform
        crs = src.crs
    
    # 2. 平滑处理
    elevation_smooth = gaussian_filter(elevation, sigma=smooth_sigma)
    
    # 3. 绘制等高线
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(16, 6))
    
    # 等高线图
    contour = ax1.contour(elevation_smooth, levels=contour_levels, 
                          colors='brown', linewidths=0.8)
    ax1.clabel(contour, inline=True, fontsize=8)
    contour_fill = ax1.contourf(elevation_smooth, levels=contour_levels, 
                                cmap='terrain')
    plt.colorbar(contour_fill, ax=ax1, label='海拔 (m)')
    ax1.set_title('等高线图')
    ax1.set_xlabel('X (像素)')
    ax1.set_ylabel('Y (像素)')
    
    # 4. 3D渲染
    rows, cols = elevation_smooth.shape
    x = np.arange(0, cols, 1)
    y = np.arange(0, rows, 1)
    X, Y = np.meshgrid(x, y)
    
    ax2 = fig.add_subplot(122, projection='3d')
    surf = ax2.plot_surface(X, Y, elevation_smooth, 
                            cmap='terrain', 
                            rstride=4, cstride=4, 
                            alpha=0.9)
    ax2.view_init(elev=30, azim=135)
    fig.colorbar(surf, ax=ax2, shrink=0.5, label='海拔 (m)')
    ax2.set_title('3D地形渲染')
    ax2.set_xlabel('X (像素)')
    ax2.set_ylabel('Y (像素)')
    ax2.set_zlabel('海拔 (m)')
    
    plt.tight_layout()
    plt.show()
    
    return elevation_smooth, transform, crs

# 使用示例
# dem_data, transform, crs = process_dem('terrain.tif')

这个函数把读取、平滑、可视化全包了。你拿到新数据，直接调这个函数，两分钟就能看到地形全貌。我个人习惯把平滑参数sigma和等高线层数levels做成可调参数，不同地形用不同配置。

好了，地形图读取和可视化就聊到这儿。记住，数据读对了，后面分析才有意义。下次咱们聊聊风数据清洗，那又是另一番天地。

---