地形预处理：裁剪、重采样、投影转换、无效值处理

各位同行，咱们直接进入正题。地形预处理，说白了就是给原始DEM数据「洗个澡、理个发、换身衣服」。我刚开始做风电项目时，拿到SRTM数据就直接扔进模型里跑，结果出来的风资源图惨不忍睹——后来才发现，地形数据没处理好，后面全是白费功夫。

今天咱们就聊聊这四个核心步骤：裁剪、重采样、投影转换、NoData处理。每一步都有坑，我踩过的坑，你就不用再踩了。

1. 裁剪：别让无关数据拖累你

原始DEM数据通常覆盖很大范围，比如一个省甚至整个国家。你只需要风场周边几十公里，剩下的全是冗余。我个人习惯的做法是：先画一个风场边界，然后向外扩展5-10公里作为缓冲区，再用这个缓冲区去裁剪DEM。

为什么会这样？因为风模拟时，地形对气流的扰动范围远超风场边界。我见过有人直接用风场红线裁剪，结果边缘区域的风速计算偏差很大——嗯，后来他重新跑了一遍，多花了三天时间。

代码示例（GDAL）：

# 用矢量边界裁剪DEM
gdalwarp -cutline wind_farm_boundary.shp \
         -crop_to_cutline \
         -of GTiff \
         input_dem.tif \
         clipped_dem.tif

# 如果边界是GeoJSON格式
gdalwarp -cutline boundary.geojson \
         -crop_to_cutline \
         -of GTiff \
         input_dem.tif \
         clipped_dem.tif

2. 重采样：统一分辨率，别让网格打架

重采样这个事，说白了就是让你的DEM分辨率跟计算网格对齐。比如你的风模型跑的是100米网格，但DEM是30米分辨率——这时候就需要重采样到100米。

我常用的重采样方法有三种：

方法	适用场景	优缺点
最邻近法	分类数据、坡度坡向	速度快，但会丢失细节
双线性插值	连续地形、一般风场	平滑过渡，计算适中
三次卷积	高精度要求、复杂山地	精度最高，但慢

代码示例：

# 重采样到100米分辨率，使用双线性插值
gdalwarp -tr 100 100 \
         -r bilinear \
         -of GTiff \
         clipped_dem.tif \
         resampled_dem.tif

# 检查重采样后的分辨率
gdalinfo resampled_dem.tif | grep "Pixel Size"

3. 投影转换：坐标系统一是底线

这个问题我吃过亏。有一次项目，DEM是WGS84经纬度坐标，风场边界是UTM投影，我没做转换直接叠在一起算——结果风场位置偏了将近2公里。你想想看，这要是建了风机，基础位置全错。

我个人建议：所有数据统一到UTM投影或高斯-克吕格投影。为什么？因为经纬度坐标下，一个像素代表的实际距离在不同纬度不一样，这会导致后续计算出现系统性偏差。

投影转换的步骤：

• 确认原始DEM的投影信息（用gdalinfo查看）
• 确定目标投影（根据风场所在区域选择UTM带号）
• 执行投影转换

# 查看原始投影
gdalinfo resampled_dem.tif | grep "Coordinate System"

# 转换到UTM Zone 50N（中国东部常用）
gdalwarp -t_srs EPSG:32650 \
         -r bilinear \
         -of GTiff \
         resampled_dem.tif \
         projected_dem.tif

# 验证转换结果
gdalinfo projected_dem.tif | grep -E "Coordinate System|Pixel Size"

4. 无效值处理：NoData不是小事

DEM数据中经常有空洞区域，比如水体、云层遮挡、雷达阴影等。这些NoData值如果不处理，风模型会把它当成0米高程——这显然不对。

我曾经接手过一个项目，原始DEM在湖泊区域全是NoData，前一个人直接填了0。结果那个区域的风速模拟值异常偏高，因为模型把湖面当成了平地，忽略了水面粗糙度的影响。

处理NoData的常用方法：

• 插值填充：用周围有效值进行插值，适合小面积空洞
• 掩膜处理：如果空洞面积太大（比如整个湖泊），建议用水体掩膜单独处理
• 替换为固定值：仅适用于你知道确切高程的场景（比如海面填0）

# 查看NoData值
gdalinfo projected_dem.tif | grep "NoData"

# 用周围像素插值填充NoData（GDAL 3.0+）
gdal_fillnodata.py projected_dem.tif \
                   filled_dem.tif \
                   -md 10  # 最大搜索距离（像素）

# 或者用gdalwarp直接替换NoData
gdalwarp -srcnodata -9999 \
         -dstnodata -9999 \
         -of GTiff \
         projected_dem.tif \
         filled_dem.tif

预处理流程总结

好了，四个步骤讲完了。我个人习惯的完整流程是这样的：

1. 先裁剪，去掉无关区域，减少数据量
2. 再重采样，统一到目标分辨率
3. 然后投影转换，确保坐标系一致
4. 最后处理NoData，填补空洞

这个顺序不能乱。你想想看，如果先做投影转换再做裁剪，投影后的数据范围会变化，裁剪边界可能对不上。我刚开始时犯过这个错，后来就长记性了。

---