WRF中使用SRTM高分辨率的地形资料

目前WRF自带的地形数据分辨率最高为30弧秒(约900m)，对于大部分研究和应用，此分辨率的地形数据基本能够满足。但是对于一些特殊应用，如超高分辨率的模拟或者大涡模拟，WRF 静态数据目前提供的地形数据分辨率稍显不足，因此此类应用需要解决如何添加更高分辨率的地形数据进入 WRF 的静态数据中。

SRTM(Shuttle Radar Topography Mission)，是由美国 NASA 主导的全球遥感探测计划，全球的采样间隔为3弧秒(约90m)，美国本土的采样间隔为1弧秒(约30m)，覆盖全球80%的范围，其数据的具体介绍与下载可参考：https://srtm.csi.cgiar.org/。

SRTM 的数据为 Geotiff 或 Esri ASCII 格式，瓦片形式(5°×5°或30°×30°)，而 geogrid 所需的静态数据也为瓦片形式存放，不过使用的是二进制格式，因此需要将 Geotiff 转为 geogrid 的二进制格式，本文将详细介绍其中的步骤。

1 SRTM数据下载

下载 Geotiff 格式，瓦片大小为5°×5°，将中国及其周边区域的数据都下载下来，下载地址：https://srtm.csi.cgiar.org/

选取的SRTM瓦片

由于 geogrid 的二进制数据维度(x或y轴)最大为99999，而分辨率太高所以一次制作的区域不能太大，此次选取下图中的16个瓦片数据，基本覆盖了中国东部地区。可以多制作几个区域，分开使用。

中国东部的SRTM瓦片数据

2 将 SRTM 数据合并

由于下载的 SRTM 数据为多个分散瓦片，因此需要将它们合成一个文件，使用 GDAL 库完成合并(如果只下载一个瓦片数据则不需要合并，可跳过此步骤)。GDAL(Geospatial Data Abstraction Library) 是一个在 X/MIT 许可协议下的开源栅格空间数据转换库。

（1）GDAL 的安装和使用说明

# conda安装gdal，时间较久
conda install -c conda-forge gdal -y # gdal_merge.py是gdal自带的命令（不是自己写的脚本）# gdal_merge.py: https://gdal.org/programs/gdal_merge.htmlgdal_merge.py [-o out_filename] [-of out_format] [-co NAME=VALUE]*              [-ps pixelsize_x pixelsize_y] [-tap] [-separate] [-q] [-v] [-pct]              [-ul_lr ulx uly lrx lry] [-init "value [value...]"]              [-n nodata_value] [-a_nodata output_nodata_value]              [-ot datatype] [-createonly] input_files

（2）GDAL 合并数据

#将下载的srtm合并成一个# srtm/srtm_china/*tif 下载的srtm的路径# -o 输出文件，得到srtm_china.tif# -a_nodata 缺省值，srtm缺省值为-32678gdal_merge.py srtm/srtm_china/*tif -o srtm_china.tif -a_nodata -32768

3 将 Geotiff 转为 geogrid 的二进制格式

（1）安装依赖库 GeoTIFF and LibTIFF

sudo apt-get install libgeotiff-dev

（2）安装 convert_geotiff

下载源码：

https://github.com/openwfm/convert_geotiff/releases

export CONVERT_GEOTIFF=$APPS/convert_geotiff_v0.1

export PATH=$CONVERT_GEOTIFF/bin:$PATH./configure --prefix=$CONVERT_GEOTIFF

make

make install

（3）将 Geotiff 格式转成 geogrid 所需的二进制格式

新建一个 topo_3s 文件夹，使用 convert_geotiff 进行格式转换。

mkdir topo_3s

cd topo_3s

convert_geotiff -w 4 -t 1500 -u "meters MSL" -d "3s topography" -b 0 -m -32768 srtm_china.tif# -w word length# -t 瓦片文件大小# -u 单位# -d 描述# -b # m 缺省值

在 topo_3s 文件夹下，生成 geogrid 所需二进制格式的数据以及一个 index 文件，将 topo_3s 文件夹复制或者移动到WRF静态数据文件夹 GEOG。

得到 topo_3s，geogrid 能够读取的二进制数据以及一个 index 文件

4 WPS 中的设置
（1）修改 GEOGRID.TBL.ARW

在 GEOGRID.TBL.ARW 文件中找到 name = HGT_M 部分，加入以下语句：

# GEOGRID.TBL.ARW 

interp_option=gtopo_3s:average_gcell(4.0)+four_pt+average_4pt

rel_path  = gtopo_3s:topo_3s

（2）namelist.wps中的设置

geog_fata_res 设置中，对应的 domain 的列加入 gtopo_3s。​​​​​​​

# namelist.wps

geog_data_res = 'default','default','gtopo_3s+default'

5 结果展示

下图可以看到，相比于使用30弧秒的数据(左图)，使用3弧秒(右图)的数据对模拟区域的地形刻画更加细致。

上图为200m分辨率的网格，共600×600格点，左图和右图所用的SRTM数据分辨率分别为为30s(约900m)，右图为3s(约90m)

6 发现的问题

SRTM 数据在海洋区域，取缺省值(-32768)，按照以上步骤操作，最终生成的 geo_em.d0x 文件中，海洋区域的地形高度值为32768，说明这个错误出现和缺省值有关。

解决方案如下：因为 geogrid 的地形数据为距海平面高度的值，海面上应该为0。在使用 convert_geotiff 之前，先将 tif 数据中的缺省值改为0，也就是将海洋部分的值修改为0，然后在运行 convert_geotiff 时将缺省值 -m 设置为0(虽然缺省值默认为0，但是运行命令 -m 0 也不能省略，原因未知，不过不设置亲测会发生错误)。

如果不涉及海洋区域，可以忽略此问题，以下为更改缺省值为0的脚本：

（代码可以观测气海同途公众号进行获取）

扫描下方二维码关注气海同途，大气海洋数值模式，可视化等内容。

参考：

https://wrfexplorer.wordpress.com/2015/03/24/high-resolution-topography/

http://www.meteoboy.com/90m-topo-in-WRF.html

https://github.com/openwfm/convert_geotiff

https://cloud.tencent.com/developer/article/1618268