数值模式 | WRF台风模拟地形敏感性实验教程，以“格美”为例

用户11172986

发布于 2026-05-19 18:49:10

1090

文章涵盖初始场数据下载，地形修改、wrf全流程以及数据可视化。很适合入门

WRF台风模拟地形敏感性实验教程，以“格美”为例

WRF（Weather Research and Forecasting）数值模式是当前最广泛使用的中尺度气象模拟系统，其通过嵌套网格技术显著提升低空气象要素的模拟精度。本次实操，我们以2024年台风“格美”为例，基于WRF（Weather Research and Forecasting）模式，结合NCL工具，系统介绍如何设计并执行一套完整的地形敏感性实验。实验通过修改地形数据来研究地形对台风路径影响。

实验通过地形敏感性分析，验证了WRF模式在模拟台风与地形相互作用方面的能力，并强调了地形数据精度在台风预报中的重要性。

一、实验简介

本教程将指导您在超算互联网平台使用WRF模式对台风“格美”进行地形敏感性实验，通过修改地形数据来研究地形对台风路径影响。实验将使用NCL（NCAR Command Language）工具修改地形数据，并比较不同地形条件下的台风模拟结果。

二、实验准备

1、环境配置

首先在超算互联网平台搜索wrf4.5，点击购买，点击去使用。通过命令行快速打开软件。首先加载必要的环境模块并设置Python环境：

module use /public/software/modules/base /public/software/modules/apps
source /work/home/username/apprepo/ncl/6.6.2-gcc4.8.5/scripts/env.sh
source /work/home/username/apprepo/anaconda3/2024.02.1-none/scripts/env.sh
conda create -n wrfdeal python==3.9
conda activate wrfdeal
conda install -c conda-forge wrf-python
pip install netcdf4
pip install meteva

2、数据准备

下载FNL初始场数据。使用以下Python脚本下载FNL数据：

import sys
import os
from datetime import datetime, timedelta
from urllib.request import urlretrieve

def generate_gdas1_urls(start_date, end_date):
    base_url = "https://data.rda.ucar.edu/d083002/grib2/{year}/{year}.{month:02d}/fnl_{year}{month:02d}{day:02d}_{hour:02d}_00.grib2"
    urls = []
    current_date = start_date
    while current_date <= end_date:
        if current_date.hour in [0, 6, 12, 18]:
            url = base_url.format(
                year=current_date.year,
                month=current_date.month,
                day=current_date.day,
                hour=current_date.hour
            )
            urls.append(url)
            print(url)
        current_date += timedelta(hours=6)
    return urls

def download_files(urls):
    for url in urls:
        filename = os.path.basename(url)
        print(f"Downloading {filename} ... ", end='', flush=True)
        try:
            urlretrieve(url, filename)
            print("done")
        except Exception as e:
            print(f"Failed: {e}")

if __name__ == "__main__":
    if len(sys.argv) != 3:
        print("Usage: python download_gdas1.py <start_date> <end_date>")
        print("Date format should be YYYYMMDDHH, for example: 2008072512")
        sys.exit(1)
    start_str = sys.argv[1]
    end_str = sys.argv[2]
    try:
        start_date = datetime.strptime(start_str, "%Y%m%d%H")
        end_date = datetime.strptime(end_str, "%Y%m%d%H")
        if start_date > end_date:
            print("Error: Start date must not be later than end date.")
            sys.exit(1)
        file_urls = generate_gdas1_urls(start_date, end_date)
        download_files(file_urls)
    except ValueError as e:
        print(f"Invalid date format: {e}")
        print("Date format should be YYYYMMDDHH, for example: 2008072512")
        sys.exit(1)

运行脚本下载数据：

python download_gdas1.py 2008072512 2008073012

三、实验流程

1. 创建项目目录

mkdir -p Project/WRF Project/WPS
cd Project

2. 准备WRF和WPS

cp -rf /work/home/username/apprepo/wrf_wps/4.5-intelmpi2017/app/WPS/ WRF
cp -rf /work/home/username/apprepo/wrf_wps/4.5-intelmpi2017/app/WRF/ WPS

3. 修改地形数据

使用脚本修改地形数据：

import netCDF4 as nc
import numpy as np
import wrf

# 1. 打开 geo_em.d01 文件
file_dir = '/work/home/acin2esoa8/project/WPS/geo_em.d01.nc'
data = nc.Dataset(file_dir, 'r+')

# 2. 读取经纬度变量（XLAT_M 和 XLONG_M）
lat = data['XLAT_M'][0, :, :]   # 纬度（2D数组）
lon = data['XLONG_M'][0, :, :]  # 经度（2D数组）

# 3. 定义台湾省的经纬度范围（大致范围）
taiwan_lat_min, taiwan_lat_max = 21.5, 25.5   # 台湾纬度范围
taiwan_lon_min, taiwan_lon_max = 119.0, 123.0 # 台湾经度范围

# 4. 找到台湾省范围内的网格点（生成一个布尔掩码）
mask = ((lat >= taiwan_lat_min) & (lat <= taiwan_lat_max) &
        (lon >= taiwan_lon_min) & (lon <= taiwan_lon_max))

# 5. 修改 HGT_M 的值（示例：台湾省地形高度降低 30%）
hgt = data['HGT_M'][:]          # 读取原始地形数据
hgt[0, mask] *= 0.5             # 仅修改台湾省范围内的地形
data['HGT_M'][:] = hgt          # 写回文件

# 6. 关闭文件
data.close()
print("台湾省地形高度修改完成！")

运行python脚本：

Python ter_mod.py

4. 运行WPS

修改namelist.wps：

cd WPS
vi namelist.wps

示例配置：

&share
 wrf_core = 'ARW',
 max_dom = 1,
 start_date = '2024-07-22_00:00:00',
 end_date   = '2024-07-25_18:00:00',
 interval_seconds = 10800
/
&geogrid
 parent_id         = 1,
 parent_grid_ratio = 1,
 i_parent_start    = 1,
 j_parent_start    = 1,
 e_we              = 260,
 e_sn              = 211,
 geog_data_res     = 'default',
 dx                = 9000,
 dy                = 9000,
 map_proj          = 'lambert',
 ref_lat           = 21.00,
 ref_lon           = 128.00,
 truelat1          = 20.0,
 truelat2          = 20.0,
 stand_lon         = 129.6,
 geog_data_path    = '/public/software/meteorology/WPS_GEOG'
/
&ungrib
 out_format = 'WPS',
 prefix = 'FILE',
/
&metgrid
 fg_name = 'FILE'
/

运行WPS预处理：

./geogrid.exe
ln -sf ungrib/Variable_Tables/Vtable.GFS Vtable
./link_grib.csh ../fnl_*
./ungrib.exe
./metgrid.exe

5. 运行WRF

修改namelist.input：

cd ../WRF/run
vi namelist.input

示例配置：

&time_control
 start_year           = 2024, 2019,
 start_month          = 07, 09,
 start_day            = 22, 04,
 start_hour           = 00, 12,
 end_year             = 2024, 2019,
 end_month            = 07, 09,
 end_day              = 25, 06,
 end_hour             = 18, 00,
 interval_seconds     = 10800
 input_from_file      = .true.,.true.,
 history_interval     = 60, 60,
 frames_per_outfile   = 1, 1,
 restart              = .false.,
 restart_interval     = 7200,
 io_form_history      = 2
 io_form_restart      = 2
 io_form_input        = 2
 io_form_boundary     = 2
/
&domains
 time_step            = 15,
 time_step_fract_num  = 0,
 time_step_fract_den  = 1,
 max_dom              = 1,
 e_we                 = 260, 220,
 e_sn                 = 211, 214,
 e_vert               = 45, 45,
 dzstretch_s          = 1.1
 p_top_requested      = 5000,
 num_metgrid_levels   = 34,
 num_metgrid_soil_levels = 4,
 dx                   = 9000, 9000,
 dy                   = 9000, 9000,
 grid_id              = 1, 2,
 parent_id            = 0, 1,
 i_parent_start       = 1, 53,
 j_parent_start       = 1, 25,
 parent_grid_ratio    = 1, 3,
 parent_time_step_ratio = 1, 3,
 feedback             = 1,
 smooth_option        = 0
/
&physics
 mp_physics           = 6, -1,
 cu_physics           = 0, -1,
 ra_lw_physics        = 1, -1,
 ra_sw_physics        = 1, -1,
 bl_pbl_physics       = 1, -1,
 sf_sfclay_physics    = 1, -1,
 sf_surface_physics   = 2, -1,
 radt                 = 15, 15,
 bldt                 = 0, 0,
 cudt                 = 0, 0,
 icloud               = 1,
 num_land_cat         = 21,
 sf_urban_physics     = 0, 0,
 fractional_seaice    = 1,
/
&fdda
/
&dynamics
 hybrid_opt           = 2,
 w_damping            = 0,
 diff_opt             = 2, 2,
 km_opt               = 4, 4,
 diff_6th_opt         = 0, 0,
 diff_6th_factor      = 0.12, 0.12,
 base_temp            = 290.
 damp_opt             = 3,
 zdamp                = 5000., 5000.,
 dampcoef             = 0.2, 0.2,
 khdif                = 0, 0,
 kvdif                = 0, 0,
 non_hydrostatic      = .true., .true.,
 moist_adv_opt        = 1, 1,
 scalar_adv_opt       = 1, 1,
 gwd_opt              = 1, 0,
/
&bdy_control
 spec_bdy_width       = 5,
 specified            = .true.
/
&grib2
/
&namelist_quilt
 nio_tasks_per_group  = 0,
 nio_groups           = 1,
/

运行WRF：

ln -sf ../../WPS/met_em.d* .
./real.exe
# 提交作业
cp /work/home/username/apprepo/wrf_wps/4.5-intelmpi2017/case/wrf.slurm .
sbatch wrf.slurm

6. 控制实验

重复上述步骤，不修改地形即可。

四、结果分析

使用Python脚本分析不同地形条件下的台风路径和强度：

import numpy as np
from netCDF4 import Dataset
from wrf import (to_np, getvar, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords)
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
import glob

# 定义绘图函数
def plot_typhoon_track(wrf_files, titles):
    plt.figure(figsize=(15, 10))
    ax = plt.axes(projection=ccrs.PlateCarree())
    colors = ['red', 'green']  # Only two colors now
    for i, (files, title) in enumerate(zip(wrf_files, titles)):
        lons, lats, slps = [], [], []
        for file in sorted(files):
            ncfile = Dataset(file)
            slp = getvar(ncfile, "slp")
            smooth_slp = smooth2d(slp, 3, cenweight=4)
            # 找到最低气压点
            min_idx = np.unravel_index(np.argmin(to_np(smooth_slp)), smooth_slp.shape)
            lon = to_np(ncfile['XLONG'][0,:,:])[min_idx]
            lat = to_np(ncfile['XLAT'][0,:,:])[min_idx]
            lons.append(lon)
            lats.append(lat)
            slps.append(to_np(smooth_slp)[min_idx])
        # 绘制路径
        ax.plot(lons, lats, color=colors[i], marker='o', linestyle='-', linewidth=2, markersize=6, label=title, transform=ccrs.PlateCarree())
    # 添加地图要素
    ax.coastlines(resolution='10m')
    ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
    ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--')
    ax.set_extent([115, 130, 15, 30])
    plt.title('Typhoon Track Comparison with Different Terrain Conditions')
    plt.legend()
    plt.savefig('terrain_sensitivity_track.png', dpi=300, bbox_inches='tight')
    plt.close()

# 使用glob自动查找文件
orig_files = sorted(glob.glob('wrfout_orig_d01_*'))   # 原始地形结果文件
half_files = sorted(glob.glob('wrfout_half_d01_*'))   # 减半地形结果文件
# 绘图比较 (only two scenarios now)
plot_typhoon_track([orig_files, half_files], ['Original Terrain', 'Half-height Terrain'])

import numpy as np
from netCDF4 import Dataset
from wrf import (to_np, getvar, smooth2d)
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker
from matplotlib.colors import ListedColormap
from meteva.base.tool.plot_tools import add_china_map_2basemap

def setup_map(ax, extent):
    """Setup map with common elements using meteva's China map"""
    # Add China map layers
    add_china_map_2basemap(ax, name="river", edgecolor='k', lw=0.5, encoding='gbk')     # Rivers
    add_china_map_2basemap(ax, name="nation", edgecolor='k', lw=0.5, encoding='gbk')   # National borders
    add_china_map_2basemap(ax, name="province", edgecolor='k', lw=0.5, encoding='gbk') # Provincial borders
    # Set extent
    ax.set_extent(extent, crs=ccrs.PlateCarree())
    # Configure gridlines
    gl = ax.gridlines(draw_labels=True, linewidth=1, color='none', alpha=0.5, linestyle='--', x_inline=False, y_inline=False)
    gl.top_labels = False
    gl.right_labels = False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    gl.rotate_labels = False
    gl.xlocator = mticker.FixedLocator(np.arange(117, 128, 2))
    gl.ylocator = mticker.FixedLocator(np.arange(20, 28, 1))
    gl.xlabel_style = {'size': 12}
    gl.ylabel_style = {'size': 12}
    return ax

# Color settings
rgb = ([237, 237, 237], [209, 209, 209], [173, 173, 173], [131, 131, 131], [93, 93, 93],
       [151, 198, 223], [111, 176, 214], [49, 129, 189], [26, 104, 174], [8, 79, 153],
       [62, 168, 91], [110, 193, 115], [154, 214, 149], [192, 230, 185], [223, 242, 217],
       [255, 255, 164], [255, 243, 0], [255, 183, 0], [255, 123, 0], [255, 62, 0],
       [255, 2, 0], [196, 0, 0], [136, 0, 0])
colors = np.array(rgb)/255.
cmap = ListedColormap(colors)
rain_levels = [0.1, 1, 2, 5, 7.5, 10, 13, 16, 20, 25, 30, 35, 40, 50, 60, 70, 80, 90, 100, 125, 150, 175, 200, 250]
rain_levels = [value * 5 for value in rain_levels]

# Read data
ncfile = Dataset('/data/wrfout_d01_2008-07-27_21_00_00')
lon = np.array(ncfile['XLONG'])[0,:,:]
lat = np.array(ncfile['XLAT'])[0,:,:]
R = (to_np(getvar(ncfile, "RAINC")) + to_np(getvar(ncfile, "RAINNC")) + to_np(getvar(ncfile, "RAINSH")))
slp = getvar(ncfile, "slp")
smooth_slp = smooth2d(slp, 3, cenweight=4)

# Plot rainfall
fig1, ax1 = plt.subplots(figsize=(10, 8), subplot_kw={'projection': ccrs.PlateCarree()})
cf1 = ax1.contourf(lon, lat, R, levels=rain_levels, cmap=cmap, transform=ccrs.PlateCarree(), extend='max')
ax1 = setup_map(ax1, [117, 128, 20, 28])
cbar1 = fig1.colorbar(cf1, ax=ax1, ticks=rain_levels[::3], shrink=0.65)
cbar1.set_label('Rainfall (mm)', size=12)
plt.title('Accumulated Rainfall', fontsize=14, pad=20)
plt.savefig('rainfall_plot.png', dpi=300, bbox_inches='tight')

# Plot SLP
fig2, ax2 = plt.subplots(figsize=(10, 8), subplot_kw={'projection': ccrs.PlateCarree()})
cf2 = ax2.contourf(lon, lat, to_np(smooth_slp), levels=10, transform=ccrs.PlateCarree())
ax2 = setup_map(ax2, [117, 128, 20, 28])
cbar2 = fig2.colorbar(cf2, ax=ax2, shrink=0.65)
cbar2.set_label('Sea Level Pressure (hPa)', size=12)
plt.title('Sea Level Pressure', fontsize=14, pad=20)
plt.savefig('slp_plot.png', dpi=300, bbox_inches='tight')
plt.close('all')

降水对比地形减半降水

正常地形降水

台风路径

五、实验结论

通过比较不同地形条件下的模拟结果，可以分析：

台风路径的变化：降低地形高度后，台风的路径发生了明显改变。与原始地形条件下的路径相比，减半地形条件下的台风路径整体偏南，并且在某些时段出现了打转现象，这表明地形对台风的引导气流和移动方向具有显著影响。
台风强度的影响：海平面气压（SLP）的对比图显示，地形高度的改变对台风的强度也产生了影响。在地形减半的模拟中，台风的最低海平面气压更低，这暗示了地形对台风结构和强度的调制作用。
降水分布的差异：累积降水量的对比图揭示了地形对降水分布和强度的重要影响。在地形减半的模拟中，降水大值区集中在台湾省南部，且降水极值相对于正常地形更大。而在原始地形区域，降水分布较为均匀。

综上所述，地形在台风的路径、强度和降水分布中扮演着关键角色。本实验通过地形敏感性分析，验证了WRF模式在模拟台风与地形相互作用方面的能力，并强调了地形数据精度在台风预报中的重要性。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2026-05-16，如有侵权请联系 cloudcommunity@tencent.com 删除

数据