做中尺度气象模拟的人多多少少都有过这样的冲动把一个真实台风或者暴雨过程装进自己的电脑然后慢慢调整地形、土地利用、物理方案看老天爷会不会“配合”。WRF就是满足这种冲动的核心工具。这篇文章不是理论课而是按我实际跑通一套完整模拟的路径来写的从编译WRF开始到用GFS与ERA5准备驱动场运行台风暴雨个例再到修改土地利用和地形数据设计敏感性试验最后用Python完成专业分析和绘图。整个过程覆盖了搭建本地“天气实验室”的主要技术点适合刚开始接触WRF、或者已经能跑通但想做过程研究和敏感性试验的气象相关专业学生、科研人员和业余爱好者。1. 把“天气实验室”搭起来硬件选型、软件依赖与WRF编译实录1.1 机器配置别只看核心数内存和磁盘才是隐藏瓶颈很多新手会先问到底需要多高配置的机器。我的观点很直接一开始不必追求服务器级别但要清楚瓶颈在哪里。模拟计算确实靠CPUWRF用MPI并行核心越多跑得越快可是大部分卡脖子的地方反而不在CPU。以一次三层嵌套台风模拟为例外层d01是27 km中间d02是9 km内层d03是3 km网格点加起来几百万如果每6分钟输出一次三维场单是wrfout就能堆出几十个GB。我见过不少人在d03跑到一半的时候磁盘写满导致结果前面全浪费。所以第一优先级是存储至少预留500 GB空闲空间再考虑CPU核心数和内存。内存方面常见的三层嵌套配置建议32 GB起步如果你想把两层都设成1 km以内的高分辨率那64 GB会更从容。另外要注意温度控制很多工作站跑长时间模拟时因为散热压不住出现节点直接卡死。我自己的习惯是正式跑之前先用一个60小时的小案例连续运行测试观察CPU温度和内存占用确认稳定再上正式任务。这样做看起来慢其实省时间。1.2 Linux环境与依赖库最省事的组合方式WRF原生支持Linux和macOSWindows需要借助WSL或者虚拟机但我不推荐在Windows上折腾问题太多。选一个长期支持版Linux发行版例如Ubuntu 20.04/22.04可以少踩很多坑。编译WRF之前需要准备编译器、MPI库和NetCDF相关库。最常见的一套组合是gcc、gfortran、g 编译器mpich 或 openmpi 并行库zlib、libpng、jasper 用于GRIB数据解码HDF5、netCDF-C、netCDF-Fortran这里特别提醒一下WRF对netCDF版本比较敏感尤其是netCDF-Fortran和netCDF-C版本之间要保持兼容。我在Ubuntu上直接apt安装时遇到过版本错配计算过程中metgrid和real.exe能跑但wrf.exe在初始化阶段就报错。后来全部手动编译或者直接用conda管理的环境问题才消失。如果不想在系统库里折腾可以用conda创建独立环境来安装编译器以外的东西。实际操作是这样的wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh conda create -n wrf_env python3.9 conda activate wrf_env conda install -c conda-forge netcdf-fortran netcdf-c hdf5 mpich jasper libpng zlib然后把环境路径写入环境的变量。编译WRF前把这些变量导出export NETCDF$CONDA_PREFIX export HDF5$CONDA_PREFIX export WRF_DIR$HOME/WRF export PATH$CONDA_PREFIX/bin:$PATH export LD_LIBRARY_PATH$CONDA_PREFIX/lib:$LD_LIBRARY_PATH用conda环境的好处是库之间版本配套基本不会出现依赖缺漏的问题。缺点是增加一层文件系统开销对大规模并行来说效率略低但对学习和单机工作完全够用。1.3 WRF 4.x的编译流程与常见坑从GitHub或者官网下载WRF源码后进入目录执行./configureconfigure会列出很多选项一般x86_64 Linux配gfortran和mpich或openmpi选带dmpar的那一项也就是分布式内存并行。如果找不到合适的选项说明依赖库环境变量没设好回头检查NETCDF路径。配置完成后执行./compile em_real这个编译过程比较久几十分钟到一两个小时都正常。编译结束后检查main目录下是否生成了wrf.exe、real.exe和ndown.exe。只编译成功但没有这三个文件等于没成功。我遇到过两个比较典型的编译错误。第一个是cpp: error: unrecognized command line option这通常是编译器版本和WRF版本不匹配导致换老版本编译器或者升级WRF。第二个是netCDF路径没找到在configure时提示NetCDF not found此时确认$NETCDF是否指向conda环境根目录。还有一个容易被忽略的是若系统里同时存在多个netcdfconfigure找到的路径和编译器实际链接的路径不一致编译时就会报各种看不懂的符号错误。编译通过之后我建议先跑一下官方自带的理想化案例确认wrf.exe能正常启动再去处理真实资料的流程。这一步跳过后面出了问题会很难区分是数据库问题还是模式本身问题。2. 驱动场从官网到本地GFS与ERA5的下载策略和WPS预处理2.1 选GFS还是ERA5实时预报 vs 再分析资料WRF模拟不能用模型自己凭空起报必须由外部的大尺度数据提供初始场和边界场。目前最常用的就是GFS和ERA5。GFS是美国国家环境预报中心发布的全球预报数据时效强、更新快每天多次发布水平分辨率大约0.25度时间间隔有1小时、3小时和6小时。如果你关心的是当前正在发生的天气过程或者想做一个“准业务化”的预报试验GFS最合适。ERA5是欧洲中期天气预报中心的第五代再分析资料把历史观测和模式输出融合起来生成了一套时间连续、变量完整的数据集水平分辨率约0.25度时间间隔为1小时。ERA5的优点是稳定、完整、适合科研统计不会因为预报时效导致模拟初期就引入较大偏差。缺点是发布时间有滞后通常需要几天到几个月不适合实时个例。做典型台风暴雨过程复盘和敏感性试验ERA5是更稳妥的选择。两者下载方式有很大差别。GFS可以走科研数据服务网站或云存储用wget直接批量拉取ERA5需要注册账号拿到API密钥后通过Python的cdsapi下载。2.2 GFS首发数据的批量下载脚本GFS文件按时间路径存储以2023年某台风过程为例假设选择2023年7月21日00时作为起报时刻需要下载f000到f120的多个时次。下载命令可以写成#!/bin/bash yy2023; mm07; dd21; hh00 fhr0 while [ $fhr -le 120 ]; do fhr3$(printf %03d $fhr) URLhttps://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl?filegfs.t${hh}z.pgrb2.0p25.f${fhr3}lev_400_mbonlev_500_mbonlev_850_mbonvar_TMPonvar_UGRDonvar_VGRDonvar_HGTonvar_RHonsubregiontoplat40leftlon110rightlon130bottomlat20dir%2Fgfs.${yy}${mm}${dd}%2F${hh}%2Fatmos wget -O gfs_${fhr3}.grb2 $URL fhr$(($fhr 6)) done需要注意的是下载时最好只选择WRF运行区域范围内的子区域和必要高度层既能降低下载时间也能减少后续WPS处理的数据量。因为WRF只需要外围区域的边界强迫不需要全球完整数据。当然如果你希望稳妥一点也可以直接下载全球数据代价是磁盘占用和处理时间都会明显增加。2.3 ERA5通过CDS API下载ERA5下载最关键的是变量列表要完整。对WRF驱动来说核心变量包括气压层上的位势、温度、风场U/V、相对湿度或比湿地面气压、2米温度、2米露点温度、10米风场海表温度、海平面气压土壤温度和土壤湿度看陆面方案是否需要如果缺了任何一个关键变量ungrib阶段会报错或者生成的中间文件变量不完整后面metgrid会出现数据空洞。在CDS官网提交请求或用Python脚本批量下载时建议把时间段按小时拆开每次请求尽量控制在合理范围以免请求排队时间过长。import cdsapi c cdsapi.Client(urlhttps://cds.climate.copernicus.eu/api, key你的KEY) c.retrieve( reanalysis-era5-pressure-levels, { variable: [geopotential, temperature, u_component_of_wind, v_component_of_wind, relative_humidity], pressure_level: [1000, 975, 950, 925, 900, 850, 800, 700, 600, 500, 400, 300, 250, 200], product_type: reanalysis, data_format: grib, year: 2023, month: 07, day: [20, 21, 22, 23], time: 00:00, }, era5_pressure_202307.grib)下载完成后务必检查文件大小和时间步数避免某个时次数据缺失。一个常见问题是ERA5的netCDF格式下载到本地后变量名和GRIB2不一样处理方式也不同。我建议统一下载GRIB格式与GFS一起走 ungrib 流程省去变量名转换的麻烦。2.4 WPS链条geogrid、ungrib、metgridWPS负责把外部数据“翻译”成WRF认识的格式。完整流程是geogrid定义模拟区域和静态地理数据ungrib把GRIB/GRIB2数据解码成中间格式metgrid把中间格式数据插值到geogrid定义的网格上运行ungrib前需要把下载的GRIB文件链接到WPS目录并把相应的Vtable复制过来。GFS和ERA5使用的Vtable不同GFS一般用Vtable.GFSERA5因为也以GRIB形式提供用Vtable.ERA5或Vtable.ECMWF都行关键是变量匹配。别把Vtable搞错否则ungrib出来的中间文件会缺变量。常见报错是“Unknown record”或者“Level not found”这往往说明数据文件缺少特定层次或者Vtable指向了错误的变量。还有一次我遇到metgrid输出场为零排查半天发现是ungrib阶段中间文件的日期标记出了问题我当时手动修改了link脚本里的小时计数补上后发现是脚本时区理解错了。从实际经验看WPS链条里最能避免坑的做法是一次处理一个时次的数据先跑通再放开批量并且每个步骤都看log文件尾部确认成功。跑批处理看似方便一旦中间某个时次数据有问题整个链条会卡在相同位置排查起来反而更慢。3. 台风暴雨个例从namelist到输出判读的完整闭环3.1 选个例不是越“极端”越好让WRF跑出漂亮结果的前提是选一个典型案例资料完整、天气系统相对清晰、模式能正确重现冷暖空气交汇或涡旋结构。很多初学者一上来就选最强台风结果模式积分几小时就出现CFL报错原因可能是驱动场太复杂、地形剧烈、积云参数化不合适。建议先选择熟悉的、过程相对平稳但降水量明显的暴雨案例等模式跑顺再挑战台风过程。我选择台风个例时通常会看一下路径是否经过资料密集区域比如沿海探空站和雷达覆盖比较好的区域方便后面用观测数据做对比验证。时间窗口建议覆盖登陆前24小时到登陆后24小时既能看到台风的涡旋结构演变也能分析降水分布的时空变化。3.2 namelist.wps里的区域设计区域设计对手是模拟成败的关键。以三层嵌套为例geogrid parent_grid_ratio 1,3,3 i_parent_start 1, 40, 60 j_parent_start 1, 40, 60 e_we 200, 301, 301 e_sn 180, 301, 301 dx 27000, 9000, 3000 dy 27000, 9000, 3000 /这里需要解释一下外层d01要覆盖台风外围大尺度环境范围不能太小否则边界强迫和内部系统发展会不协调。d02负责覆盖台风主体云系d03则聚焦重点降水区或地形关键区。如果d03范围太小台风中心稍微偏离预定路径系统就会跑到模拟区域外结果就废了。所以d03中心我的经验是放到台风登陆点附近略靠海洋一侧而不是路径中心。namelist.wps的ungrib和metgrid两个部分比较固定核心是把interval_seconds设置成与驱动数据的时间间隔一致。GFS六小时一次就设21600秒ERA5一小时一次就设3600秒。3.3 namelist.input物理方案配置思路namelist.input里最让人纠结的是物理方案选择因为它决定模拟的云和降水特征。对台风案例我的常用组合是微物理方案WSM6或Thompson两者对台风暖云和冷云过程都有较好表现。d01也可以考虑新Thompson高分辨率内层更细腻。积云参数化d01必须开d02视分辨率3km以内可以不开。台风模拟d01常用Kain-Fritsch或Grell-Freitas。边界层方案YSU适合中尺度模拟MYNN在高分辨率下的边界层发展细节更好。陆面过程Noah LSM经典稳定Noah-MP可选但计算量大。辐射方案RRTMG在长波和短波上都比较好且支持嵌套域并行。时间步长不能拍脑袋选。经验公式大约是6*dx(km)27km对应约162秒9km对应约54秒3km对应约18秒。考虑到台风强对流环境我通常会再乘0.5即9km用30秒3km用10秒。时间步长过大容易CFL报错其实调小之后很多问题迎刃而解。3.4 运行流程从real.exe到wrf.exe的监控要点real.exe是连接metgrid输出到初始场和边界场的过程。运行real之前需要确保metgrid输出时间从模拟起始时刻开始一直延续到终止时刻。real.exe会读取wrfinput和wrfbdy如果时间窗口不匹配会提示找不到数据。我习惯在运行real前用ncdump快速检查一下metgrid的时间维度。wrf.exe运行过程中日志文件是rsl.out.0000和rsl.error.0000。很多人觉得看日志很麻烦但实际上大多数错误在最后几十行里就能定位。常见错误包括CFL violation通常出现在强对流回波与地形相互作用的地方。先调小时步长再考虑降低内层分辨率或调整微物理方案。土壤温度初始化报错大多是era5或gfs驱动场缺少深层土壤温度层次需要检查ungrib阶段是否包含土壤变量。内存不足直接卡死或报segmentation fault多半是嵌套域网格数太大超出了机器内存。模拟运行后先不要急着深入分析。看一下2米温度、海平面气压的基本时空分布对比前几个积分的合理性。比如台风中心气压有没有随时间明显下降云系有没有沿路径移动。如果模拟三小时就出现气压剧烈震荡多半是初始场和模式不够协调先重查namelist。判断拟结果是否合理我会用两个快速方法一是画海平面气压场和850 hPa风场看有没有清晰的涡旋结构和闭合低压中心二是画3小时累积降水看雨带是否和台风螺旋雨带位置接近。确认这两点之后再进入后面更精细的分析。4. 修改土地利用与地形给模式“动手术”的具体操作4.1 为什么要修改土地利用和地形经典数值模拟通常直接使用自带静态地理数据但真实世界的下垫面一直在变化尤其是城市化扩张、农田灌溉、水库建设、林草地退化。这些变化会影响地表能量平衡和低层风场进而改变降水分布。你如果想研究“如果这个区域变成城市台风降水会不会增强”那就必须在模拟中人为修改土地利用类型再对比控制试验和敏感试验的差异。地形修改也有类似逻辑比如你想评估风电场对局地风场的影响或者填海造陆对台风登陆降水的影响就需要用新的地形高程替换原来的地形数据。这类研究在数值模拟里非常常见也是WRF灵敏度试验的核心操作之一。4.2 理解geo_em.d01.nc的关键变量进行修改之前你需要知道WRF的静态地理数据如何组织。geogrid步骤生成geo_em.d01.nc里面包含一系列二维和三维变量。最重要的是下面几个LU_INDEX每个格点的主导土地利用类型索引直接参与陆面过程计算。LANDUSEF12个月的土地利用占比维度是(月份, 类别, 纬度, 经度)每个格点各个月份不同类别的占比相加为1。TERRAINDATA模式用的地形海拔高度默认来自USGS/SRTM等数据。SOILTEXTURE、TOPOFRAC等次要变量一般不用改。注意LU_INDEX和LANDUSEF并不是独立的。当你修改LANDUSEF时LU_INDEX必须同步更新为占比最大的类别否则陆面过程的土地利用属性会和实际不匹配。4.3 用Python直接修改geo_em文件最常见的方法是修改某个区域的LANDUSEF把原本的农田或草地类型改为城市类型。下面以USGS 24类分类为例假设城市类别编号为1农田类别为3现在要把某个矩形区域改成城市。import numpy as np from netCDF4 import Dataset src Dataset(geo_em.d01.nc, r) lons src.variables[XLONG_M][0] lats src.variables[XLAT_M][0] mask (lats 32) (lats 34) (lons 118) (lons 120) landusef src.variables[LANDUSEF][:] # (time1, month12, category, lat, lon) ncat landusef.shape[2] # 把区域范围内所有月份的土地利用改为城市类别 for m in range(12): # 先将所有类别占比清零 landusef[0, m, :, :, :] 0.0 # 城市类别设为1这里假设城市类别索引是1 landusef[0, m, 1, :, :][mask] 1.0 src.variables[LANDUSEF][:] landusef # 重新计算LU_INDEX lu_index src.variables[LU_INDEX][0, :, :] for m in range(12): lu_index[mask] 1 # 城市类别编号 src.variables[LU_INDEX][0, :, :] lu_index[mask] src.close()这段代码的关键点在于修改完LANDUSEF后一定要重算LU_INDEX而且矩形区域的经纬度范围要和geo_em里的XLONG_M、XLAT_M一致否则修改会落在区域外。如果你想按某个GIS矢量边界来做不规则区域可以用GDAL把矢量栅格化成与geo_em相同分辨率的掩膜数组再对掩膜区域赋值。4.4 地形修改与一致性检查地形修改类似先用外部的高分辨率DEM处理成模拟网格坐标插值后替换TERRAINDATA。下面是一个简化思路用GDAL读取外部DEM用Warp重投影到与geo_em相同的坐标系Resample到相同分辨率然后写入geo_em这个过程最需要注意的是坐标匹配。geo_em里的LAT/LON是经纬度网格中心点外部DEM可能使用UTM或其他投影。必须先把DEM转换到WGS84经纬度再进行双线性插值。插值完的地形要平滑处理不然模式积分时容易出现气压场剧烈变化的数值不稳定。我个人的建议是如果只是简单研究直接修改geo_em文件就够如果研究区域地形复杂或者涉及城市精细尺度建议修改之后跑一个短时模拟对比修改前后的海平面气压场和风场确认没有出现不合理噪声。5. 敏感性试验设计让WRF帮你回答“如果……会怎样”5.1 常见敏感性试验类型敏感性试验的核心是控制变量。想研究土地利用变化的影响保持驱动场和物理方案完全不变只替换新土地利用数据想研究物理过程参数化的影响保持数据不变只改变某一个物理方案。还有地形敏感性、海温敏感性、初始场扰动敏感性等。这里一定要分清楚“敏感性试验”和“普通预报试验”的区别前者是通过对照试验来诊断模式响应的因果链条后者是尽可能让预报接近真实。所以敏感性试验里人为改动之后模拟结果与真实观测偏差变大这本身不一定是坏事反而说明该因子确实在起作用。5.2 控制变量设计从CTL到多组试验举一个土地利用敏感性实验的例子。假设我想研究“城市扩张对台风暴雨的影响”可以这样设计试验名驱动场物理方案土地利用地形CTLERA5统一方案A原始土地利用原始地形URBERA5统一方案A城市扩张方案原始地形TERRERA5统一方案A原始土地利用修改后地形URBTERRERA5统一方案A城市扩张方案修改后地形这样设计的好处是CTL和URB对比能单独分离出土地利用变化的影响CTL和TERR对比能单独分离出地形变化的影响URBTERR与URB对比还能看出两个因素叠加后有没有非线性效应。实际操作时我通常为每个试验建立单独目录目录里只存放链接过来的wrfout和namelist所有试验共用同一份驱动场和初始场数据。这样能最大程度避免误操作导致试验组之间数据不一致。命名也要清晰比如CTL、URB、TERR不要用test1、test2这种名字否则跑完一周你自己都分不清哪组是哪组。5.3 如何判断敏感信号是真实的做完敏感性试验后不能只画两张图说“有差异”。要判断差异是不是真实的模式响应需要做统计检验。常用方法是对多时次输出做配对t检验比较CTL和URB在感兴趣区域的平均降水、温度或风场差异是否显著。更具体的做法是提取每个试验的逐小时降水字段选定一个固定区域比如台风暴雨落区或城市下游区域对每个时次计算区域平均得到一组时间序列用Python的scipy.stats.ttest_rel计算配对t检验如果p值小于0.05可以认为两个试验的差异在统计上显著这里有个容易踩的坑如果样本数太少比如只有几个时次差异再多也通不过显著性检验。所以设计试验时最好把模拟时段拉长或者用多个个例重复试验。多个个例的做法更符合真实科研逻辑但计算量会成倍增加。我还想强调一个职业习惯敏感性试验里最怕无意中改变了不该变的因素。比如你在修改土地利用时不小心把海表温度覆盖了或者运行real.exe时没有复制同一份驱动数据都会让灵敏度信号失真。每次跑完一组试验建议先用diff对比namelist和相关输入文件的checksum确认除了目标变量外没有其他差异。6. 用Python做专业级后处理分析从读取wrfout到论文级绘图6.1 Python环境别手动装包直接上condaWRF后处理最常用的库是wrf-python、netCDF4、xarray、matplotlib、cartopy。新手最容易在装包阶段卡住最主要原因是用pip直接装wrf-python在某些平台上很难成功。虽然现在pip装wrf-python也可以但依赖库冲突还是很多。最省心的方法是先用conda创建环境然后一条命令装上全部依赖conda create -n wrf_post python3.9 conda activate wrf_post conda install -c conda-forge wrf-python netcdf4 xarray dask matplotlib cartopy pandas scipy如果conda下载慢换成国内镜像源。这条命令装完之后matplotlib、cartopy和wrf-python都能用基本覆盖了日常绘图需求。这里特别提醒一下wrf-python依赖的numpy版本不能太新有些版本组合会出现无法导入的问题建议用conda默认解析出来的版本组合别手动升级。6.2 读取wrfout并理解网格结构wrfout文件本质是netCDF格式可以用netCDF4直接读取但直接用wrf-python的getvar更方便。下面这段代码可以快速读取台风模拟的海平面气压和10米风场import xarray as xr import wrf import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature ds xr.open_dataset(wrfout_d03_2023-07-21_00:00:00) slp wrf.getvar(ds, slp) ua wrf.getvar(ds, ua) va wrf.getvar(ds, va) lon wrf.getvar(ds, lon2d) lat wrf.getvar(ds, lat2d) proj wrf.get_cartopy(slp) fig plt.figure(figsize(10, 8)) ax plt.axes(projectionproj) ax.add_feature(cfeature.COASTLINE, linewidth1) ax.add_feature(cfeature.BORDERS, linewidth0.5) levels list(range(960, 1020, 4)) cont ax.contourf(lon, lat, slp, levelslevels, cmapviridis, transformccrs.PlateCarree()) ax.barbs(lon[::5, ::5], lat[::5, ::5], ua[::5, ::5], va[::5, ::5], transformccrs.PlateCarree()) plt.colorbar(cont, axax, shrink0.8) plt.title(Sea Level Pressure and 10m Wind) plt.savefig(slp_wind_d03.png, dpi300, bbox_inchestight)这段代码用到了wrf-python的getvar它内部已经做了坐标变换和单位换算比直接读netCDF再手动处理方便得多。要画降水用wrf.getvar(ds, RAINC) wrf.getvar(ds, RAINNC)这是累积降水差分即可得到单位时间降水量。这里的RAINC是积云参数化降水RAINNC是显式微物理降水千万别漏掉任何一个。6.3 横坐标太密集、变量名记不住怎么办这个坑来自真实经历很多新手用matplotlib画图时横坐标显示时间序列或者地图坐标轴刻度非常密集甚至叠在一起图像很难看。解决办法是设置合适的时间分辨率或者用MaxNLocator控制刻度数量。地图坐标的经纬度标注也一样手动设置间隔即可。还有一个常见问题是wrf-python的getvar变量名不是netCDF原变量名比如温度是theta画图时需要转成摄氏温度可以用tk获取开尔文温度位势高度在WRF里是z而不是GRIB文件里的HGT降水累积量是RAINC、RAINNC对应命名要记牢想快速查看所有可用变量可以执行wrf.getvar(ds, ALL)看变量清单或者直接打印ds.data_vars看项目。最好在项目里维护一个自己的变量映射表方便写脚本时查。6.4 批量绘制和统计分析跑了一组敏感性试验通常要生成几十上百张图。我习惯用Python写一个循环脚本遍历多个wrfout文件把每个时次的降水、风场、温度图批量输出。批量脚本要特别注意变量提取和保存的循环逻辑别把时次顺序搞混。在统计分析部分我还常用xarray直接读取多个试验的累计降水做差值场和区域平均ctl xr.open_mfdataset(CTL/wrfout_d03_*, combineby_coords) urb xr.open_mfdataset(URB/wrfout_d03_*, combineby_coords) # 计算每个试验在某个时间段内的累计降水 def total_precip(ds): return (ds[RAINC].sum(dimTime) ds[RAINNC].sum(dimTime)).values precip_ctl total_precip(ctl) precip_urb total_precip(urb) diff precip_urb - precip_ctl # 区域平均 region (lat 30) (lat 34) (lon 118) (lon 122) print(区域平均降水量差值(URB-CTL):, diff[region].mean())这段代码里需要注意的坑是xr.open_mfdataset的时间维度合并方式如果文件名前缀不统一或者存在重叠时次会导致计算结果错误。稳妥起见可以用glob.glob排序后传入或者先合并再检查维度长度。6.5 绘图细节从能出图到“论文级”能把图画出来和画出图文并茂的图之间差距通常在细节。具体控制点包括地图投影选择中纬度区域用Lambert Conformal小范围区域用横轴墨卡托WRF自带get_cartopy函数可以直接获得正确投影。我自己画单层嵌套时经常直接用get_cartopy输出省去手动设置投影。色标选择降水通常用precipitation色标风场用矢量和填色结合温度用红蓝渐变。用matplotlib的ColormapBuilder调整色标起点和终点避免自动配色掩盖真实值。字体设置插入图名、标签时统一设置rcParams[font.sans-serif]为支持中文的字体。如果不想为中文乱码头疼直接用英文标签也行。保存分辨率论文投稿通常要求300 dpi以上出图时指定dpi300文件格式优先PNG或PDFPDF体积小且缩放不变形。绘图还能做箱线图、误差条形图、时间序列等结合你的研究问题自由发挥。关键是形成标准化的输入-函数-输出流程让每个试验的图风格一致后面写论文和报告时能直接复用。从搭建环境到修改地形再到敏感性实验设计和Python分析这是一条完整的WRF研究路径。我实际跑通一轮之后最大的感受是WRF的难点不是运行本身而是你对数据流和物理过程的理解。每一步的取舍都会影响结果而最好的调试工具不是教程是你对日志和数据变量的熟悉程度。当你第四次遇到CFL报错、第五次发现GEODATA缺文件时这些坑就会变成一个可复用的经验库。希望这份记录能让你在搭建自己的天气实验室时少走一些弯路多一些真正理解模式的时间。 WRF中尺度气象模拟实战：从编译到敏感性试验与Python分析
做中尺度气象模拟的人多多少少都有过这样的冲动把一个真实台风或者暴雨过程装进自己的电脑然后慢慢调整地形、土地利用、物理方案看老天爷会不会“配合”。WRF就是满足这种冲动的核心工具。这篇文章不是理论课而是按我实际跑通一套完整模拟的路径来写的从编译WRF开始到用GFS与ERA5准备驱动场运行台风暴雨个例再到修改土地利用和地形数据设计敏感性试验最后用Python完成专业分析和绘图。整个过程覆盖了搭建本地“天气实验室”的主要技术点适合刚开始接触WRF、或者已经能跑通但想做过程研究和敏感性试验的气象相关专业学生、科研人员和业余爱好者。1. 把“天气实验室”搭起来硬件选型、软件依赖与WRF编译实录1.1 机器配置别只看核心数内存和磁盘才是隐藏瓶颈很多新手会先问到底需要多高配置的机器。我的观点很直接一开始不必追求服务器级别但要清楚瓶颈在哪里。模拟计算确实靠CPUWRF用MPI并行核心越多跑得越快可是大部分卡脖子的地方反而不在CPU。以一次三层嵌套台风模拟为例外层d01是27 km中间d02是9 km内层d03是3 km网格点加起来几百万如果每6分钟输出一次三维场单是wrfout就能堆出几十个GB。我见过不少人在d03跑到一半的时候磁盘写满导致结果前面全浪费。所以第一优先级是存储至少预留500 GB空闲空间再考虑CPU核心数和内存。内存方面常见的三层嵌套配置建议32 GB起步如果你想把两层都设成1 km以内的高分辨率那64 GB会更从容。另外要注意温度控制很多工作站跑长时间模拟时因为散热压不住出现节点直接卡死。我自己的习惯是正式跑之前先用一个60小时的小案例连续运行测试观察CPU温度和内存占用确认稳定再上正式任务。这样做看起来慢其实省时间。1.2 Linux环境与依赖库最省事的组合方式WRF原生支持Linux和macOSWindows需要借助WSL或者虚拟机但我不推荐在Windows上折腾问题太多。选一个长期支持版Linux发行版例如Ubuntu 20.04/22.04可以少踩很多坑。编译WRF之前需要准备编译器、MPI库和NetCDF相关库。最常见的一套组合是gcc、gfortran、g 编译器mpich 或 openmpi 并行库zlib、libpng、jasper 用于GRIB数据解码HDF5、netCDF-C、netCDF-Fortran这里特别提醒一下WRF对netCDF版本比较敏感尤其是netCDF-Fortran和netCDF-C版本之间要保持兼容。我在Ubuntu上直接apt安装时遇到过版本错配计算过程中metgrid和real.exe能跑但wrf.exe在初始化阶段就报错。后来全部手动编译或者直接用conda管理的环境问题才消失。如果不想在系统库里折腾可以用conda创建独立环境来安装编译器以外的东西。实际操作是这样的wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh conda create -n wrf_env python3.9 conda activate wrf_env conda install -c conda-forge netcdf-fortran netcdf-c hdf5 mpich jasper libpng zlib然后把环境路径写入环境的变量。编译WRF前把这些变量导出export NETCDF$CONDA_PREFIX export HDF5$CONDA_PREFIX export WRF_DIR$HOME/WRF export PATH$CONDA_PREFIX/bin:$PATH export LD_LIBRARY_PATH$CONDA_PREFIX/lib:$LD_LIBRARY_PATH用conda环境的好处是库之间版本配套基本不会出现依赖缺漏的问题。缺点是增加一层文件系统开销对大规模并行来说效率略低但对学习和单机工作完全够用。1.3 WRF 4.x的编译流程与常见坑从GitHub或者官网下载WRF源码后进入目录执行./configureconfigure会列出很多选项一般x86_64 Linux配gfortran和mpich或openmpi选带dmpar的那一项也就是分布式内存并行。如果找不到合适的选项说明依赖库环境变量没设好回头检查NETCDF路径。配置完成后执行./compile em_real这个编译过程比较久几十分钟到一两个小时都正常。编译结束后检查main目录下是否生成了wrf.exe、real.exe和ndown.exe。只编译成功但没有这三个文件等于没成功。我遇到过两个比较典型的编译错误。第一个是cpp: error: unrecognized command line option这通常是编译器版本和WRF版本不匹配导致换老版本编译器或者升级WRF。第二个是netCDF路径没找到在configure时提示NetCDF not found此时确认$NETCDF是否指向conda环境根目录。还有一个容易被忽略的是若系统里同时存在多个netcdfconfigure找到的路径和编译器实际链接的路径不一致编译时就会报各种看不懂的符号错误。编译通过之后我建议先跑一下官方自带的理想化案例确认wrf.exe能正常启动再去处理真实资料的流程。这一步跳过后面出了问题会很难区分是数据库问题还是模式本身问题。2. 驱动场从官网到本地GFS与ERA5的下载策略和WPS预处理2.1 选GFS还是ERA5实时预报 vs 再分析资料WRF模拟不能用模型自己凭空起报必须由外部的大尺度数据提供初始场和边界场。目前最常用的就是GFS和ERA5。GFS是美国国家环境预报中心发布的全球预报数据时效强、更新快每天多次发布水平分辨率大约0.25度时间间隔有1小时、3小时和6小时。如果你关心的是当前正在发生的天气过程或者想做一个“准业务化”的预报试验GFS最合适。ERA5是欧洲中期天气预报中心的第五代再分析资料把历史观测和模式输出融合起来生成了一套时间连续、变量完整的数据集水平分辨率约0.25度时间间隔为1小时。ERA5的优点是稳定、完整、适合科研统计不会因为预报时效导致模拟初期就引入较大偏差。缺点是发布时间有滞后通常需要几天到几个月不适合实时个例。做典型台风暴雨过程复盘和敏感性试验ERA5是更稳妥的选择。两者下载方式有很大差别。GFS可以走科研数据服务网站或云存储用wget直接批量拉取ERA5需要注册账号拿到API密钥后通过Python的cdsapi下载。2.2 GFS首发数据的批量下载脚本GFS文件按时间路径存储以2023年某台风过程为例假设选择2023年7月21日00时作为起报时刻需要下载f000到f120的多个时次。下载命令可以写成#!/bin/bash yy2023; mm07; dd21; hh00 fhr0 while [ $fhr -le 120 ]; do fhr3$(printf %03d $fhr) URLhttps://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl?filegfs.t${hh}z.pgrb2.0p25.f${fhr3}lev_400_mbonlev_500_mbonlev_850_mbonvar_TMPonvar_UGRDonvar_VGRDonvar_HGTonvar_RHonsubregiontoplat40leftlon110rightlon130bottomlat20dir%2Fgfs.${yy}${mm}${dd}%2F${hh}%2Fatmos wget -O gfs_${fhr3}.grb2 $URL fhr$(($fhr 6)) done需要注意的是下载时最好只选择WRF运行区域范围内的子区域和必要高度层既能降低下载时间也能减少后续WPS处理的数据量。因为WRF只需要外围区域的边界强迫不需要全球完整数据。当然如果你希望稳妥一点也可以直接下载全球数据代价是磁盘占用和处理时间都会明显增加。2.3 ERA5通过CDS API下载ERA5下载最关键的是变量列表要完整。对WRF驱动来说核心变量包括气压层上的位势、温度、风场U/V、相对湿度或比湿地面气压、2米温度、2米露点温度、10米风场海表温度、海平面气压土壤温度和土壤湿度看陆面方案是否需要如果缺了任何一个关键变量ungrib阶段会报错或者生成的中间文件变量不完整后面metgrid会出现数据空洞。在CDS官网提交请求或用Python脚本批量下载时建议把时间段按小时拆开每次请求尽量控制在合理范围以免请求排队时间过长。import cdsapi c cdsapi.Client(urlhttps://cds.climate.copernicus.eu/api, key你的KEY) c.retrieve( reanalysis-era5-pressure-levels, { variable: [geopotential, temperature, u_component_of_wind, v_component_of_wind, relative_humidity], pressure_level: [1000, 975, 950, 925, 900, 850, 800, 700, 600, 500, 400, 300, 250, 200], product_type: reanalysis, data_format: grib, year: 2023, month: 07, day: [20, 21, 22, 23], time: 00:00, }, era5_pressure_202307.grib)下载完成后务必检查文件大小和时间步数避免某个时次数据缺失。一个常见问题是ERA5的netCDF格式下载到本地后变量名和GRIB2不一样处理方式也不同。我建议统一下载GRIB格式与GFS一起走 ungrib 流程省去变量名转换的麻烦。2.4 WPS链条geogrid、ungrib、metgridWPS负责把外部数据“翻译”成WRF认识的格式。完整流程是geogrid定义模拟区域和静态地理数据ungrib把GRIB/GRIB2数据解码成中间格式metgrid把中间格式数据插值到geogrid定义的网格上运行ungrib前需要把下载的GRIB文件链接到WPS目录并把相应的Vtable复制过来。GFS和ERA5使用的Vtable不同GFS一般用Vtable.GFSERA5因为也以GRIB形式提供用Vtable.ERA5或Vtable.ECMWF都行关键是变量匹配。别把Vtable搞错否则ungrib出来的中间文件会缺变量。常见报错是“Unknown record”或者“Level not found”这往往说明数据文件缺少特定层次或者Vtable指向了错误的变量。还有一次我遇到metgrid输出场为零排查半天发现是ungrib阶段中间文件的日期标记出了问题我当时手动修改了link脚本里的小时计数补上后发现是脚本时区理解错了。从实际经验看WPS链条里最能避免坑的做法是一次处理一个时次的数据先跑通再放开批量并且每个步骤都看log文件尾部确认成功。跑批处理看似方便一旦中间某个时次数据有问题整个链条会卡在相同位置排查起来反而更慢。3. 台风暴雨个例从namelist到输出判读的完整闭环3.1 选个例不是越“极端”越好让WRF跑出漂亮结果的前提是选一个典型案例资料完整、天气系统相对清晰、模式能正确重现冷暖空气交汇或涡旋结构。很多初学者一上来就选最强台风结果模式积分几小时就出现CFL报错原因可能是驱动场太复杂、地形剧烈、积云参数化不合适。建议先选择熟悉的、过程相对平稳但降水量明显的暴雨案例等模式跑顺再挑战台风过程。我选择台风个例时通常会看一下路径是否经过资料密集区域比如沿海探空站和雷达覆盖比较好的区域方便后面用观测数据做对比验证。时间窗口建议覆盖登陆前24小时到登陆后24小时既能看到台风的涡旋结构演变也能分析降水分布的时空变化。3.2 namelist.wps里的区域设计区域设计对手是模拟成败的关键。以三层嵌套为例geogrid parent_grid_ratio 1,3,3 i_parent_start 1, 40, 60 j_parent_start 1, 40, 60 e_we 200, 301, 301 e_sn 180, 301, 301 dx 27000, 9000, 3000 dy 27000, 9000, 3000 /这里需要解释一下外层d01要覆盖台风外围大尺度环境范围不能太小否则边界强迫和内部系统发展会不协调。d02负责覆盖台风主体云系d03则聚焦重点降水区或地形关键区。如果d03范围太小台风中心稍微偏离预定路径系统就会跑到模拟区域外结果就废了。所以d03中心我的经验是放到台风登陆点附近略靠海洋一侧而不是路径中心。namelist.wps的ungrib和metgrid两个部分比较固定核心是把interval_seconds设置成与驱动数据的时间间隔一致。GFS六小时一次就设21600秒ERA5一小时一次就设3600秒。3.3 namelist.input物理方案配置思路namelist.input里最让人纠结的是物理方案选择因为它决定模拟的云和降水特征。对台风案例我的常用组合是微物理方案WSM6或Thompson两者对台风暖云和冷云过程都有较好表现。d01也可以考虑新Thompson高分辨率内层更细腻。积云参数化d01必须开d02视分辨率3km以内可以不开。台风模拟d01常用Kain-Fritsch或Grell-Freitas。边界层方案YSU适合中尺度模拟MYNN在高分辨率下的边界层发展细节更好。陆面过程Noah LSM经典稳定Noah-MP可选但计算量大。辐射方案RRTMG在长波和短波上都比较好且支持嵌套域并行。时间步长不能拍脑袋选。经验公式大约是6*dx(km)27km对应约162秒9km对应约54秒3km对应约18秒。考虑到台风强对流环境我通常会再乘0.5即9km用30秒3km用10秒。时间步长过大容易CFL报错其实调小之后很多问题迎刃而解。3.4 运行流程从real.exe到wrf.exe的监控要点real.exe是连接metgrid输出到初始场和边界场的过程。运行real之前需要确保metgrid输出时间从模拟起始时刻开始一直延续到终止时刻。real.exe会读取wrfinput和wrfbdy如果时间窗口不匹配会提示找不到数据。我习惯在运行real前用ncdump快速检查一下metgrid的时间维度。wrf.exe运行过程中日志文件是rsl.out.0000和rsl.error.0000。很多人觉得看日志很麻烦但实际上大多数错误在最后几十行里就能定位。常见错误包括CFL violation通常出现在强对流回波与地形相互作用的地方。先调小时步长再考虑降低内层分辨率或调整微物理方案。土壤温度初始化报错大多是era5或gfs驱动场缺少深层土壤温度层次需要检查ungrib阶段是否包含土壤变量。内存不足直接卡死或报segmentation fault多半是嵌套域网格数太大超出了机器内存。模拟运行后先不要急着深入分析。看一下2米温度、海平面气压的基本时空分布对比前几个积分的合理性。比如台风中心气压有没有随时间明显下降云系有没有沿路径移动。如果模拟三小时就出现气压剧烈震荡多半是初始场和模式不够协调先重查namelist。判断拟结果是否合理我会用两个快速方法一是画海平面气压场和850 hPa风场看有没有清晰的涡旋结构和闭合低压中心二是画3小时累积降水看雨带是否和台风螺旋雨带位置接近。确认这两点之后再进入后面更精细的分析。4. 修改土地利用与地形给模式“动手术”的具体操作4.1 为什么要修改土地利用和地形经典数值模拟通常直接使用自带静态地理数据但真实世界的下垫面一直在变化尤其是城市化扩张、农田灌溉、水库建设、林草地退化。这些变化会影响地表能量平衡和低层风场进而改变降水分布。你如果想研究“如果这个区域变成城市台风降水会不会增强”那就必须在模拟中人为修改土地利用类型再对比控制试验和敏感试验的差异。地形修改也有类似逻辑比如你想评估风电场对局地风场的影响或者填海造陆对台风登陆降水的影响就需要用新的地形高程替换原来的地形数据。这类研究在数值模拟里非常常见也是WRF灵敏度试验的核心操作之一。4.2 理解geo_em.d01.nc的关键变量进行修改之前你需要知道WRF的静态地理数据如何组织。geogrid步骤生成geo_em.d01.nc里面包含一系列二维和三维变量。最重要的是下面几个LU_INDEX每个格点的主导土地利用类型索引直接参与陆面过程计算。LANDUSEF12个月的土地利用占比维度是(月份, 类别, 纬度, 经度)每个格点各个月份不同类别的占比相加为1。TERRAINDATA模式用的地形海拔高度默认来自USGS/SRTM等数据。SOILTEXTURE、TOPOFRAC等次要变量一般不用改。注意LU_INDEX和LANDUSEF并不是独立的。当你修改LANDUSEF时LU_INDEX必须同步更新为占比最大的类别否则陆面过程的土地利用属性会和实际不匹配。4.3 用Python直接修改geo_em文件最常见的方法是修改某个区域的LANDUSEF把原本的农田或草地类型改为城市类型。下面以USGS 24类分类为例假设城市类别编号为1农田类别为3现在要把某个矩形区域改成城市。import numpy as np from netCDF4 import Dataset src Dataset(geo_em.d01.nc, r) lons src.variables[XLONG_M][0] lats src.variables[XLAT_M][0] mask (lats 32) (lats 34) (lons 118) (lons 120) landusef src.variables[LANDUSEF][:] # (time1, month12, category, lat, lon) ncat landusef.shape[2] # 把区域范围内所有月份的土地利用改为城市类别 for m in range(12): # 先将所有类别占比清零 landusef[0, m, :, :, :] 0.0 # 城市类别设为1这里假设城市类别索引是1 landusef[0, m, 1, :, :][mask] 1.0 src.variables[LANDUSEF][:] landusef # 重新计算LU_INDEX lu_index src.variables[LU_INDEX][0, :, :] for m in range(12): lu_index[mask] 1 # 城市类别编号 src.variables[LU_INDEX][0, :, :] lu_index[mask] src.close()这段代码的关键点在于修改完LANDUSEF后一定要重算LU_INDEX而且矩形区域的经纬度范围要和geo_em里的XLONG_M、XLAT_M一致否则修改会落在区域外。如果你想按某个GIS矢量边界来做不规则区域可以用GDAL把矢量栅格化成与geo_em相同分辨率的掩膜数组再对掩膜区域赋值。4.4 地形修改与一致性检查地形修改类似先用外部的高分辨率DEM处理成模拟网格坐标插值后替换TERRAINDATA。下面是一个简化思路用GDAL读取外部DEM用Warp重投影到与geo_em相同的坐标系Resample到相同分辨率然后写入geo_em这个过程最需要注意的是坐标匹配。geo_em里的LAT/LON是经纬度网格中心点外部DEM可能使用UTM或其他投影。必须先把DEM转换到WGS84经纬度再进行双线性插值。插值完的地形要平滑处理不然模式积分时容易出现气压场剧烈变化的数值不稳定。我个人的建议是如果只是简单研究直接修改geo_em文件就够如果研究区域地形复杂或者涉及城市精细尺度建议修改之后跑一个短时模拟对比修改前后的海平面气压场和风场确认没有出现不合理噪声。5. 敏感性试验设计让WRF帮你回答“如果……会怎样”5.1 常见敏感性试验类型敏感性试验的核心是控制变量。想研究土地利用变化的影响保持驱动场和物理方案完全不变只替换新土地利用数据想研究物理过程参数化的影响保持数据不变只改变某一个物理方案。还有地形敏感性、海温敏感性、初始场扰动敏感性等。这里一定要分清楚“敏感性试验”和“普通预报试验”的区别前者是通过对照试验来诊断模式响应的因果链条后者是尽可能让预报接近真实。所以敏感性试验里人为改动之后模拟结果与真实观测偏差变大这本身不一定是坏事反而说明该因子确实在起作用。5.2 控制变量设计从CTL到多组试验举一个土地利用敏感性实验的例子。假设我想研究“城市扩张对台风暴雨的影响”可以这样设计试验名驱动场物理方案土地利用地形CTLERA5统一方案A原始土地利用原始地形URBERA5统一方案A城市扩张方案原始地形TERRERA5统一方案A原始土地利用修改后地形URBTERRERA5统一方案A城市扩张方案修改后地形这样设计的好处是CTL和URB对比能单独分离出土地利用变化的影响CTL和TERR对比能单独分离出地形变化的影响URBTERR与URB对比还能看出两个因素叠加后有没有非线性效应。实际操作时我通常为每个试验建立单独目录目录里只存放链接过来的wrfout和namelist所有试验共用同一份驱动场和初始场数据。这样能最大程度避免误操作导致试验组之间数据不一致。命名也要清晰比如CTL、URB、TERR不要用test1、test2这种名字否则跑完一周你自己都分不清哪组是哪组。5.3 如何判断敏感信号是真实的做完敏感性试验后不能只画两张图说“有差异”。要判断差异是不是真实的模式响应需要做统计检验。常用方法是对多时次输出做配对t检验比较CTL和URB在感兴趣区域的平均降水、温度或风场差异是否显著。更具体的做法是提取每个试验的逐小时降水字段选定一个固定区域比如台风暴雨落区或城市下游区域对每个时次计算区域平均得到一组时间序列用Python的scipy.stats.ttest_rel计算配对t检验如果p值小于0.05可以认为两个试验的差异在统计上显著这里有个容易踩的坑如果样本数太少比如只有几个时次差异再多也通不过显著性检验。所以设计试验时最好把模拟时段拉长或者用多个个例重复试验。多个个例的做法更符合真实科研逻辑但计算量会成倍增加。我还想强调一个职业习惯敏感性试验里最怕无意中改变了不该变的因素。比如你在修改土地利用时不小心把海表温度覆盖了或者运行real.exe时没有复制同一份驱动数据都会让灵敏度信号失真。每次跑完一组试验建议先用diff对比namelist和相关输入文件的checksum确认除了目标变量外没有其他差异。6. 用Python做专业级后处理分析从读取wrfout到论文级绘图6.1 Python环境别手动装包直接上condaWRF后处理最常用的库是wrf-python、netCDF4、xarray、matplotlib、cartopy。新手最容易在装包阶段卡住最主要原因是用pip直接装wrf-python在某些平台上很难成功。虽然现在pip装wrf-python也可以但依赖库冲突还是很多。最省心的方法是先用conda创建环境然后一条命令装上全部依赖conda create -n wrf_post python3.9 conda activate wrf_post conda install -c conda-forge wrf-python netcdf4 xarray dask matplotlib cartopy pandas scipy如果conda下载慢换成国内镜像源。这条命令装完之后matplotlib、cartopy和wrf-python都能用基本覆盖了日常绘图需求。这里特别提醒一下wrf-python依赖的numpy版本不能太新有些版本组合会出现无法导入的问题建议用conda默认解析出来的版本组合别手动升级。6.2 读取wrfout并理解网格结构wrfout文件本质是netCDF格式可以用netCDF4直接读取但直接用wrf-python的getvar更方便。下面这段代码可以快速读取台风模拟的海平面气压和10米风场import xarray as xr import wrf import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature ds xr.open_dataset(wrfout_d03_2023-07-21_00:00:00) slp wrf.getvar(ds, slp) ua wrf.getvar(ds, ua) va wrf.getvar(ds, va) lon wrf.getvar(ds, lon2d) lat wrf.getvar(ds, lat2d) proj wrf.get_cartopy(slp) fig plt.figure(figsize(10, 8)) ax plt.axes(projectionproj) ax.add_feature(cfeature.COASTLINE, linewidth1) ax.add_feature(cfeature.BORDERS, linewidth0.5) levels list(range(960, 1020, 4)) cont ax.contourf(lon, lat, slp, levelslevels, cmapviridis, transformccrs.PlateCarree()) ax.barbs(lon[::5, ::5], lat[::5, ::5], ua[::5, ::5], va[::5, ::5], transformccrs.PlateCarree()) plt.colorbar(cont, axax, shrink0.8) plt.title(Sea Level Pressure and 10m Wind) plt.savefig(slp_wind_d03.png, dpi300, bbox_inchestight)这段代码用到了wrf-python的getvar它内部已经做了坐标变换和单位换算比直接读netCDF再手动处理方便得多。要画降水用wrf.getvar(ds, RAINC) wrf.getvar(ds, RAINNC)这是累积降水差分即可得到单位时间降水量。这里的RAINC是积云参数化降水RAINNC是显式微物理降水千万别漏掉任何一个。6.3 横坐标太密集、变量名记不住怎么办这个坑来自真实经历很多新手用matplotlib画图时横坐标显示时间序列或者地图坐标轴刻度非常密集甚至叠在一起图像很难看。解决办法是设置合适的时间分辨率或者用MaxNLocator控制刻度数量。地图坐标的经纬度标注也一样手动设置间隔即可。还有一个常见问题是wrf-python的getvar变量名不是netCDF原变量名比如温度是theta画图时需要转成摄氏温度可以用tk获取开尔文温度位势高度在WRF里是z而不是GRIB文件里的HGT降水累积量是RAINC、RAINNC对应命名要记牢想快速查看所有可用变量可以执行wrf.getvar(ds, ALL)看变量清单或者直接打印ds.data_vars看项目。最好在项目里维护一个自己的变量映射表方便写脚本时查。6.4 批量绘制和统计分析跑了一组敏感性试验通常要生成几十上百张图。我习惯用Python写一个循环脚本遍历多个wrfout文件把每个时次的降水、风场、温度图批量输出。批量脚本要特别注意变量提取和保存的循环逻辑别把时次顺序搞混。在统计分析部分我还常用xarray直接读取多个试验的累计降水做差值场和区域平均ctl xr.open_mfdataset(CTL/wrfout_d03_*, combineby_coords) urb xr.open_mfdataset(URB/wrfout_d03_*, combineby_coords) # 计算每个试验在某个时间段内的累计降水 def total_precip(ds): return (ds[RAINC].sum(dimTime) ds[RAINNC].sum(dimTime)).values precip_ctl total_precip(ctl) precip_urb total_precip(urb) diff precip_urb - precip_ctl # 区域平均 region (lat 30) (lat 34) (lon 118) (lon 122) print(区域平均降水量差值(URB-CTL):, diff[region].mean())这段代码里需要注意的坑是xr.open_mfdataset的时间维度合并方式如果文件名前缀不统一或者存在重叠时次会导致计算结果错误。稳妥起见可以用glob.glob排序后传入或者先合并再检查维度长度。6.5 绘图细节从能出图到“论文级”能把图画出来和画出图文并茂的图之间差距通常在细节。具体控制点包括地图投影选择中纬度区域用Lambert Conformal小范围区域用横轴墨卡托WRF自带get_cartopy函数可以直接获得正确投影。我自己画单层嵌套时经常直接用get_cartopy输出省去手动设置投影。色标选择降水通常用precipitation色标风场用矢量和填色结合温度用红蓝渐变。用matplotlib的ColormapBuilder调整色标起点和终点避免自动配色掩盖真实值。字体设置插入图名、标签时统一设置rcParams[font.sans-serif]为支持中文的字体。如果不想为中文乱码头疼直接用英文标签也行。保存分辨率论文投稿通常要求300 dpi以上出图时指定dpi300文件格式优先PNG或PDFPDF体积小且缩放不变形。绘图还能做箱线图、误差条形图、时间序列等结合你的研究问题自由发挥。关键是形成标准化的输入-函数-输出流程让每个试验的图风格一致后面写论文和报告时能直接复用。从搭建环境到修改地形再到敏感性实验设计和Python分析这是一条完整的WRF研究路径。我实际跑通一轮之后最大的感受是WRF的难点不是运行本身而是你对数据流和物理过程的理解。每一步的取舍都会影响结果而最好的调试工具不是教程是你对日志和数据变量的熟悉程度。当你第四次遇到CFL报错、第五次发现GEODATA缺文件时这些坑就会变成一个可复用的经验库。希望这份记录能让你在搭建自己的天气实验室时少走一些弯路多一些真正理解模式的时间。