WRF中地形数据（海拔高度）分辨率最高为30s，差不多就是900 m，当模型空间分辨率较高时，比如在低于1 km的情况下，经常会考虑增加地形高度的分辨率。

数据概述

SRTM(Shuttle Radar Topography Mission)，是由美国 NASA 主导的全球遥感探测计划，全球的采样间隔为3弧秒(约90m)，美国本土的采样间隔为1弧秒(约30m)，覆盖全球80%的范围。数据的具体介绍与下载可参考-SRTM 90m DEM Digital Elevation Database。

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

数据下载

下载SRMT 30m地形数据，下载地址：SRTM Data

替代地址：SRTM Tile Grabber

选取中国地区瓦片数据，进行下载：

选中需要下载的瓦片，点击【Search】，会出现以下界面：

可采用插件（Downthemall）批量下载所有瓦片数据，如下：

注意： 由于 geogrid 的二进制数据维度(x或y轴)最大为99999，而分辨率太高所以一次制作的区域不能太大

本次研究合并的瓦片数据如下：共9个瓦片数据

数据处理

打开下载的压缩包，其中有readme.txt和后缀为 .hdr、.tfw、.tif的4个文件，根据readme.txt说明，这里下载的是4.1版本，.hdr，.tfw存储了投影信息。

合并多个SRTM 数据-GDAL库

由于下载的 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.html
gdal_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

详细安装步骤可参见另一博客-【Python库安装】Python环境安装GDAL库。

（2）GDAL 合并数据

在win10上合并.tif文件，键入命令（在cmd窗口下运行）如下：

#将下载的srtm合并成一个
# srtm/srtm_china/*tif 下载的srtm的路径
# -o 输出文件，得到srtm_china.tif
# -a_nodata 缺省值，srtm缺省值为-32678
python D:\Anaconda\Scripts\gdal_merge.py -o srtm_GBA.tif -a_nodata -32768 --optfile D:\SRTM_DEM\tif_list.txt

conda activate myenv
python D:\Anaconda\envs\myenv\Scripts\gdal_merge.py -o srtm_GBA.tif -a_nodata -32768 --optfile D:\SRTM_DEM\tif_list.txt

gdal_merge.py文件路径在执行环境键入gdal_merge.py以后就会自动出现了，照着复制一下就好。

最后的参数–optfile D:\SRTM_DEM\tif_list.txt是指定输入的合并.tif文件，tif_list.txt内容如下：

D:\SRTM_DEM\srtm_58_07\srtm_58_07.tif
D:\SRTM_DEM\srtm_58_08\srtm_58_08.tif
D:\SRTM_DEM\srtm_58_09\srtm_58_09.tif
D:\SRTM_DEM\srtm_59_07\srtm_59_07.tif
D:\SRTM_DEM\srtm_59_08\srtm_59_08.tif
D:\SRTM_DEM\srtm_59_09\srtm_59_09.tif
D:\SRTM_DEM\srtm_60_07\srtm_60_07.tif
D:\SRTM_DEM\srtm_60_08\srtm_60_08.tif
D:\SRTM_DEM\srtm_60_09\srtm_60_09.tif

执行成功屏幕会出现以下信息：

并在当前目录下生成了srtm_GBA.tif。

说明： 由于需要GDAL库的安装，我是在Anaconda环境下运行的，还需要下载numpy库，在尝试多个版本后，修改为1.24.0版本（适配Python3.8-3.11版本）后运行成功

转为geogrid二进制格式

（1）安装依赖库 GeoTIFF and LibTIFF

sudo apt-get install libgeotiff-dev

（2）安装 convert_geotiff
convert_geotiff安装步骤可参见另一博客-【WRF工具】服务器上安装convert_geotiff。

（3）将 Geotiff 格式转成 geogrid 所需的二进制格式
新建一个 topo_3s 文件夹，使用 convert_geotiff 进行格式转换。

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'

数据对比

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

海洋区域缺省值

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

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

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

from osgeo import gdal
import numpy as np

def read_tiff(inpath):
  ds=gdal.Open(inpath, gdal.GA_ReadOnly)
  row=ds.RasterXSize
  col=ds.RasterYSize
  band=ds.RasterCount
  geoTransform = ds.GetGeoTransform()
  proj = ds.GetProjection()
  data=np.zeros([row,col,band])
  for i in range(band):
   dt=ds.GetRasterBand(1)
   data[:,:,i]=dt.ReadAsArray(0,0,col,row)
  return data[:,:,0], geoTransform, proj
  
def array2raster(outpath,array,geoTransform,proj):
 cols=array.shape[1]
 rows=array.shape[0]
 driver=gdal.GetDriverByName('Gtiff')
 #outRaster=driver.Create(outpath,cols,rows,1,gdal.GDT_Byte) # 0-255
 outRaster=driver.Create(outpath,cols,rows,1, gdal.GDT_Float32)
 outRaster.SetGeoTransform(geoTransform)#参数2,6为水平垂直分辨率，参数3,5表示图片是指北的
 outband=outRaster.GetRasterBand(1)
 outband.WriteArray(array)
 outRaster.SetProjection(proj)#将几何对象的数据导出为wkt格式
 outRaster.FlushCache()
 
if __name__=="__main__":
    data, geoTransform, proj = read_tiff('srtm_china.tif')
    data[data==np.min(data)] = 0.0
    array2raster("srtm_china2.tif",data,geoTransform, proj)

参考

1、知乎-WRF中如何使用SRTM的3s高分辨率地形数据集
2、WRF中使用SRTM高分辨率的地形资料