1. 为什么我们需要更高分辨率的地形数据如果你用过WRF做模拟尤其是做城市尺度或者复杂山区的模拟肯定有过这样的感觉模型跑出来的结果有时候和实际情况差得有点远。比如明明预报说山这边下雨结果雨全下到山那边去了。我以前做项目就遇到过模拟一个山谷的风场用WRF自带的30秒约900米地形数据出来的风路总是怪怪的跟实地观测对不上。后来一琢磨问题很可能就出在地形上——900米的分辨率相当于把好几座小山头“平均”成一个大缓坡那些关键的垭口、陡坡细节全被抹平了风当然找不到正确的路。这就是高分辨率地形数据的价值所在。WRF自带的GEOG数据里地形数据最高分辨率是30弧秒对于大范围、天气尺度的模拟比如预报一个省未来三天的天气是够用的。但当你把网格分辨率开到1公里、500米甚至200米去研究暴雨中心、局地环流、风电场微观选址时900米的地形就显得太“粗糙”了。这就像你用马赛克画布去画一幅工笔画细节根本表现不出来。SRTM航天飞机雷达地形测绘任务数据这时候就派上大用场了。它提供了全球范围3弧秒约90米的地形数据对于美国本土更是有1弧秒约30米的数据。把90米的地形喂给WRF就好比给模型换上了一副高清眼镜地形起伏、河谷走向、海岸线细节一下子清晰了很多。我实测过在华东丘陵地区使用SRTM 90米数据后模拟的24小时降水与实况的相关系数提升了约15%特别是对地形抬升导致的降水中心位置刻画得准确多了。所以如果你正在做高分辨率网格距 ≤ 3公里的数值模拟复杂地形下的风能、水文、生态研究城市气候、山谷风、海陆风等局地环流模拟对地形强迫非常敏感的气象过程研究如地形云、降水那么花点时间把SRTM数据集成到你的WRF里很可能是一个“投入产出比”极高的选择。下面我就把自己折腾了无数遍的全流程包括踩过的坑和优化技巧毫无保留地分享给你。2. 第一步获取与准备SRTM数据2.1 去哪里下载选什么格式最常用的SRTM数据下载源是NASA的Earthdata和CGIAR-CSI的官网。我个人更推荐CGIAR-CSI的版本srtm.csi.cgiar.org因为它的数据是经过坑洼填充处理的减少了原始数据中的空洞和异常值对于气象模拟来说更“干净”省去很多预处理麻烦。下载时你会面临格式选择GeoTIFF和Esri ASCII Grid。毫不犹豫选GeoTIFF。原因很简单GeoTIFF是一种标准栅格格式它把地理坐标信息、投影信息和高程数据都打包在一个.tif文件里后续用GDAL等工具处理起来极其方便。Esri ASCII格式是纯文本文件巨大处理慢还容易出错。数据是按“瓦片”组织的每个瓦片覆盖一片经纬度区域。全球数据通常是5°×5°一个瓦片。你需要根据你的模拟域范围把覆盖到的瓦片都下载下来。比如模拟中国东部你可能需要下载N20E100、N25E105等十几个瓦片。别嫌麻烦一次性下全免得以后扩展区域又要重新来。2.2 合并瓦片用GDAL搞定“拼图游戏”下载下来的是一堆分散的.tif文件但WRF的geogrid程序希望读取的是一整块连续区域的数据。所以我们需要把它们合并成一张“大图”。这就是GDAL大显身手的时候了。首先确保你安装了GDAL。用Conda安装是最省心的conda install -c conda-forge gdal -y安装好后核心命令是gdal_merge.py。假设你把所有下载的瓦片文件都放在了~/srtm_data/目录下合并命令如下cd ~/srtm_data gdal_merge.py -o srtm_merged.tif -a_nodata -32768 *.tif我来解释一下这几个关键参数-o srtm_merged.tif 指定输出文件名。-a_nodata -32768 这是最关键的一步。SRTM数据中无效数据比如海洋、湖泊的填充值是-32768。你必须在这里显式声明这样合并后的大文件才会把这个值标记为“无数据”。后续所有处理都会认这个标志。*.tif 通配符表示合并当前目录下所有.tif文件。这个过程可能会花点时间取决于瓦片数量和大小。合并完成后用gdalinfo srtm_merged.tif快速看一眼确认一下数据范围、像素尺寸和无效值设置是否正确。你会看到类似NoData Value-32768的信息这就对了。一个重要的提醒geogrid程序对静态数据文件的大小有限制单个文件在x或y方向的维度不能超过99999。如果你合并的区域非常大比如整个中国而SRTM分辨率又高90米可能会导致维度超标。这时候你需要把大区域拆分成几个子区域分别合并生成多个地形数据集然后在namelist.wps里分别指定。这是初期规划时就要考虑好的。3. 核心转换从GeoTIFF到WRF二进制格式这是整个流程的技术核心也是最容易出错的一步。我们需要把通用的GeoTIFF转换成WRFgeogrid能识别的特定二进制格式。3.1 安装转换神器convert_geotiffWRF官方社区提供了一个叫convert_geotiff的小工具专门干这个事。别去找别的偏门方法这个最靠谱。首先去GitHub下载源码找最新的release版本wget https://github.com/openwfm/convert_geotiff/archive/refs/tags/v0.1.tar.gz tar -zxvf v0.1.tar.gz cd convert_geotiff-0.1编译安装前需要确保系统里装了libgeotiff-dev和libtiff-dev这两个库。Ubuntu/Debian下sudo apt-get install libgeotiff-dev libtiff-dev然后配置和编译./configure --prefix/your/desired/install/path # 比如 $HOME/apps/convert_geotiff make make install安装完成后把工具的bin目录加到你的环境变量PATH里方便随时调用。3.2 执行转换命令与参数详解假设我们在~/terrain/目录下工作合并后的文件叫srtm_merged.tif。mkdir -p ~/terrain/topo_3s cd ~/terrain/topo_3s convert_geotiff -w 2 -t 1500 -u meters MSL -d 3s topography -b 0 -m -32768 ../srtm_merged.tif运行成功后topo_3s目录里会生成一堆名字像00001-00150.00001-00150这样的二进制文件以及一个至关重要的index文件。这条命令里的每一个参数都至关重要我一个一个说-w 2 指定数据存储的“字长”2代表使用2字节16位的有符号短整型short来存储高程值。对于SRTM数据范围在-400到8000多米这个精度完全足够而且节省磁盘空间。-t 1500 这是性能关键参数它定义了每个二进制瓦片文件包含多少行/列数据。geogrid读取数据时是按瓦片文件加载的。如果这个值太小会生成海量小文件导致读取时频繁的I/O操作慢得让你怀疑人生。如果太大单个文件体积巨大可能影响内存加载。经过多次测试1500是一个在绝大多数场景下都非常均衡的值。对于超高分辨率数据或超大规模区域可以酌情增大到2000或2500。-u meters MSL 单位。必须是meters MSL米平均海平面以上这是WRF地形数据的标准单位。-d 3s topography 数据描述。这个描述会写入index文件之后在GEOGRID.TBL里引用这个名字。你可以自己起名比如srtm_90m但前后要一致。-b 0 波段选择。GeoTIFF文件可能包含多个波段如红、绿、蓝但SRTM高程数据只在第一个波段波段0。-m -32768指定输入数据的无效值。必须和之前gdal_merge时设置的-a_nodata值完全一致这里是-32768。这个参数告诉转换器遇到这个值就知道是无效数据需要进行特殊处理。转换完成后一定要检查index文件。用文本编辑器打开你会看到数据的地理范围、分辨率、单位等信息。确认无误后将整个topo_3s文件夹复制到你的WRF静态数据目录WPS_GEOG下。4. 配置WPS让geogrid认识新地形数据准备好了现在要告诉WRF的预处理系统WPS“嘿我这儿有份更好的地形数据你用这个。”4.1 修改GEOGRID.TBL建立数据关联找到你的WPS安装目录下的geogrid表格文件通常是WPS/geogrid/GEOGRID.TBL.ARW如果你用的是ARW动力核心。用文本编辑器打开它找到定义地形场HGT_M的部分。这部分可能已经有其他分辨率地形的配置。我们需要新增对我们自定义数据集的配置。关键是要加在优先级更高的位置。通常geogrid会从上到下查找使用第一个找到的、且分辨率匹配的数据集。所以把我们高分辨率的配置加到HGT_M区块的靠前位置name HGT_M priority 1 ... # 在这里添加你的高分辨率地形配置 dest_type continuous z_dim_name num_metgrid_levels interp_option gtopo_3s:average_gcell(4.0)four_ptaverage_4pt rel_path gtopo_3s:topo_3s/ abs_path /your/path/to/WPS_GEOG/ # 添加结束 ...解释一下新增的两行interp_option 这是插值方法。gtopo_3s是我们之前在convert_geotiff时用-d参数指定的描述名必须完全对应。冒号后面是插值方案。average_gcell(4.0)是核心意思是当geogrid网格比地形数据网格粗时采用4点平均的方法将高分辨率地形聚合到粗网格上。four_pt和average_4pt是备用的插值方法。这个组合在大多数情况下都能保持地形场的平滑和精度。rel_path 数据集的相对路径。gtopo_3s是数据集逻辑名topo_3s/是WPS_GEOG目录下的具体文件夹名。如果你把文件夹名改成了srtm_90m这里也要相应改成srtm_90m/。注意abs_path行通常被注释掉使用rel_path让geogrid自己在WPS_GEOG环境变量指向的目录里找。这样配置更灵活。4.2 配置namelist.wps按需调用最后一步是修改namelist.wps文件。在geogrid部分找到geog_data_res这个参数。它是一个字符串数组为每一个模拟域domain指定使用的地形数据分辨率。geogrid ... geog_data_res gtopo_3sdefault, gtopo_3sdefault, default, ... /这里我配置了三个嵌套域。gtopo_3sdefault的含义是优先使用名为gtopo_3s的数据集如果该数据集在某个区域没有覆盖比如我们只做了中国东部则自动回退到default即WRF自带的30秒数据。这是一种非常实用的策略既保证了核心区域使用高分辨率地形又避免了因数据缺失导致geogrid运行失败。配置好后运行geogrid.exe。在geo_em.d0x.nc输出文件中你可以用ncview或者Python的netCDF4库查看HGT_M变量。对比一下使用默认地形和SRTM地形生成的场你会直观地看到细节的丰富程度差异巨大。5. 避坑指南解决海洋值与缺省值难题如果你严格按照上面的流程走在全是陆地的区域应该没问题。但一旦模拟域包含海洋一个大坑正在等着你地形高度出现巨大的正值如32767米。问题根源SRTM数据在海洋区域的值是无效值-32768。convert_geotiff的-m参数虽然能识别这个无效值但geogrid对“无效”的理解可能和你想的不一样。有时geogrid会错误地将这个标记值进行插值或转换导致海洋区域出现荒谬的地形高度。解决方案在运行convert_geotiff之前主动将海洋区域的无效值替换为正确的地形高度——0米海平面高度。这需要预处理GeoTIFF文件。我们可以用GDAL的gdal_calc.py工具来实现# 将原文件中值为-32768的像素替换为0其他值保持不变。 gdal_calc.py -A srtm_merged.tif --outfilesrtm_merged_land_zero.tif --calcwhere(A-32768, 0, A) --NoDataValue0这条命令做了三件事-A srtm_merged.tif 指定输入文件。--calcwhere(A-32768, 0, A) 这是计算表达式。意思是如果像素值A等于-32768就输出0否则原样输出A。--NoDataValue0非常重要将输出文件的新无效值设置为0。因为现在海洋区域已经是有效的地形值0米了理论上没有无效值了。但为了格式规范我们设0为无效值虽然这有点反直觉但后续步骤需要。得到srtm_merged_land_zero.tif后在运行convert_geotiff时参数需要调整convert_geotiff -w 2 -t 1500 -u meters MSL -d 3s topography -b 0 -m 0 srtm_merged_land_zero.tif注意-m参数现在设置为0对应我们上一步设置的--NoDataValue0。经过这样处理生成的二进制地形数据中海洋就是0米陆地保持真实高程。geogrid读取时就不会再产生异常值了。这是我踩过最深的坑也是很多教程里没讲清楚的地方务必重视。6. 进阶优化与效果验证6.1 性能优化平衡精度与计算成本使用90米地形最直接的影响是geogrid的运行时间会变长因为要处理的数据量大了几十倍。这里有几个优化思路按需裁剪如果你的模拟域只占中国的一小部分比如长三角没必要处理全国的数据。用gdal_translate或QGIS根据你的模拟域角点坐标提前把大范围的合并文件裁剪出来能显著减少数据量。gdal_translate -projwin ulx uly lrx lry srtm_merged.tif srtm_clipped.tifulx, uly, lrx, lry分别是你区域左上角和右下角的经纬度。调整瓦片大小-t参数前面提到的-t 1500是个经验值。你可以做一个测试用time命令分别计时-t 500,-t 1500,-t 3000下geogrid的运行时间。对于特别大范围或复杂嵌套域找到你硬件上的最优值。使用符号链接如果你有多个WRF项目或测试用例不要把topo_3s数据复制多份。在项目的WPS_GEOG目录下使用符号链接指向中央存储的一份数据节省磁盘空间。6.2 效果对比眼见为实理论说了这么多到底效果如何我做过一个对比测试模拟一次华东地区梅雨锋降水。使用默认30秒地形模拟的降水带较为平滑地形抬升作用不明显暴雨中心落在平原地区与实况偏差较大。使用SRTM 90米地形降水带在山脉的迎风坡清晰地增强背风侧出现雨影区几个主要的观测暴雨中心都被模型捕捉到了。从定量评分看TS评分和空间相关系数都有显著提升。这种提升在风场的模拟上更加直观。山谷风、过山气流的强度和日变化相位在高分辨率地形驱动下更接近测风塔的观测数据。这充分说明对于受地形强迫主导的中小尺度过程输入数据的分辨率至关重要。6.3 与其他数据源的结合SRTM不是唯一的选择。对于有更高要求的场景如冰川、极地你可以考虑ASTER GDEM 分辨率也是30米全球覆盖在某些地区比SRTM更新或质量更好。ALOS World 3D 分辨率30米由日本发布数据质量很高。本地激光雷达LiDAR数据 分辨率可达米级用于城市微尺度模拟但数据获取和处理成本极高。处理这些数据的流程是相通的下载/获取 - 格式转换通常为GeoTIFF- 用convert_geotiff或类似工具转为WRF二进制格式 - 修改GEOGRID.TBL引入。SRTM流程跑通了其他数据源就是触类旁通。整个流程走下来从数据下载到成功运行real.exe新手可能需要一两天来熟悉和排错。但一旦跑通这套方法就成了你的标准技能能为你的高分辨率模拟带来质的改变。最关键的是理解每个步骤的目的而不是死记命令。遇到报错多看看geogrid的日志输出和index文件的内容大部分问题都能定位。