简介本资源面向气象建模初学者与WRF用户提供一套轻量级、即用型的wrfout文件变量提取与可视化实践方案解决WRF模拟输出数据难以直接读取和分析的常见痛点。压缩包仅2KB含2个核心脚本t2m.pyPython实现基于netCDF4/xarray提取地表温度T2M并绘图与wrfout.nclNCL实现调用WRF内置函数读取三维气象场、支持变量快速抽取与地理投影绘图分别体现两种主流工具在数据解析、坐标处理及可视化上的典型范式。已有1449人学习下载适合科研人员、气象预报员及高校师生开展WRF后处理入门训练。读者可直接运行脚本理解wrf.out二进制格式结构掌握温度、风速等关键变量的索引逻辑、时空维度切片方法及地理坐标系适配技巧为后续复杂分析或跨语言协同处理奠定实操基础。1. 从 wrfout 文件里“抠”出温度、风速、湿度为什么 NCL 写脚本越来越难而 Python wrf-python 成了气象数据提取的默认路径你刚跑完 WRF 模拟生成了几十 GB 的wrfout_d01_2023-06-15_12:00:00文件想快速拿到某一层比如 850 hPa的温度场、水平风分量U/V、相对湿度用于画图或输入下游模型——但打开 NCL 手册翻到第 47 页还在找wrf_user_getvar的参数顺序ncdump -h看了一眼变量名就头晕T2,TK,T,QVAPOR,RH,XLAT,XLONG,XTIME…到底哪个才是你要的真实温度哪个是格点经纬度哪个带时间维度哪个要先插值再取层更糟的是NCL 编译环境一升级就报错ncl_fortran.so找不到libnetcdf.so.13版本不匹配……而隔壁同事用 Python 三行代码就抽出了 3D 温度场并转成 Pandas DataFrame还顺手做了垂直插值。这不是玄学是 wrf-python 库把 WRF 输出格式NetCDF WRF 特有元数据封装成了可预测、可调试、可版本管理的 Python 对象。它不替代 WRF 本身但彻底替代了你手动解析Times,XLAT,XLONG,DX,DY,CEN_LAT,CEN_LON,MOAD_CEN_LAT这些隐藏字段的血泪经验。适合所有正在用 WRF 做中尺度模拟、集合预报、后处理分析且已切换到 Python 生态或正打算切的业务/科研用户——尤其当你需要把 wrfout 提取逻辑嵌入自动化流程、CI/CD 流水线或和 PyTorch、Xarray、Cartopy 联动时这条路径不是“更方便”而是唯一可持续的。2. 用 wrf-python 在本地跑通 wrfout 变量提取最小依赖、最简命令、最稳输出wrf-python 不是“另一个 NetCDF 读取器”。它专为 WRF 输出设计能自动识别 WRF 的投影类型Lambert Conformal、Mercator、Polar Stereographic、坐标系参数、垂直坐标eta-levels 或 pressure-levels、时间戳编码方式并提供getvar()这个智能接口——它知道T是三维位温tc是二维 2m 温度uvmet是旋转后的 U/V 分量rh是相对湿度甚至dbz是雷达反射率。你不用查手册记变量名映射表它内部已硬编码了 WRF 4.x/5.x 的全部标准变量别名与物理含义。2.1 安装 wrf-python避开 conda-forge 与 pip 的版本陷阱wrf-python 依赖 netCDF4、numpy、scipy 和 pyproj但不能直接pip install wrf-python——PyPI 上的包早已停更最后更新是 2020 年且不兼容 Python 3.9。正确路径是# 推荐用 conda 安装稳定、预编译、无编译错误 conda install -c conda-forge wrf-python # 若必须用 pip如无 conda 环境需指定 GitHub 主干分支2024 年仍活跃维护 pip install githttps://github.com/NCAR/wrf-python.gitmain提示conda install会自动拉取netcdf41.6.4,pyproj3.9.0,numpy1.24.4等兼容组合若用 pip务必先pip install netcdf4 pyproj scipy再装 wrf-python否则import wrf时会因pyproj版本过高≥4.0导致CRS初始化失败——这是 2024 年最常翻车的第一步。验证安装是否成功import wrf print(wrf.__version__) # 应输出 2.1.0 或更高GitHub main 分支当前为 2.1.12.2 读取 wrfout 文件并提取基础变量三行代码拿到经纬度网格与海平面气压假设你的文件路径是/data/wrfout_d01_2023-06-15_12:00:00目标是提取XLAT,XLONG,PSFC地表气压from netCDF4 import Dataset import wrf # 1. 用 netCDF4 打开文件wrf-python 依赖此对象 ncfile Dataset(/data/wrfout_d01_2023-06-15_12:00:00) # 2. 提取地理网格自动识别投影返回二维 numpy 数组 lats wrf.getvar(ncfile, lat, timeidxwrf.ALL_TIMES) # shape: (time, south_north, west_east) lons wrf.getvar(ncfile, lon, timeidxwrf.ALL_TIMES) # 3. 提取地表气压单位Pa自动去单位转换 psfc wrf.getvar(ncfile, PSFC, timeidxwrf.ALL_TIMES) # shape: (time, south_north, west_east)关键参数说明timeidxwrf.ALL_TIMES提取所有时间步若只取第 0 个时刻用timeidx0wrf.getvar()内部会调用wrf.get_proj_params(ncfile)解析MAP_PROJ,CEN_LAT,CEN_LON,TRUELAT1,TRUELAT2,MOAD_CEN_LAT等字段自动生成pyproj.CRS对象无需你手动构造投影返回值是xarray.DataArray若安装了 xarray或numpy.ndarray默认带.attrs[units],.attrs[description]等元数据比裸ncfile.variables[PSFC][:]多 10 倍可用信息2.3 提取三维变量并做垂直插值拿到 850 hPa 温度场的完整链路WRF 输出的是 eta 坐标σ 坐标变体T是位温K不是气温P是气压Pa。要得到某等压面如 850 hPa上的温度必须做垂直插值——wrf-python 封装了interplevel()函数底层调用scipy.interpolate.interp1d但自动处理了 eta 层权重、边界外推、多时间步广播# 提取三维位温 Tshape: time, bottom_top, south_north, west_east t wrf.getvar(ncfile, T, timeidxwrf.ALL_TIMES) # 单位K # 提取三维气压 P单位Pa p wrf.getvar(ncfile, P, timeidxwrf.ALL_TIMES) # 修正气压perturbation pb wrf.getvar(ncfile, PB, timeidxwrf.ALL_TIMES) # 基准气压base state pressure p pb # 总气压单位 Pa # 插值得到 850 hPa 层上的位温注意850 hPa 85000 Pa t_850 wrf.interplevel(t, pressure, 85000, missingnp.nan) # 若需气温非位温用 wrf.tv_from_wvap() 或 wrf.tk_from_t() 转换见 4.2 节逻辑说明interplevel()要求field被插值场和vertical垂直坐标场维度一致且vertical必须是气压Pa或高度m——这里pressure是(time, bottom_top, south_north, west_east)完全匹配第四参数85000是目标气压值单位必须与pressure一致missingnp.nan指定插值失败处填 NaN如地形高于 850 hPa 的格点输出t_850是(time, south_north, west_east)已降维可直接画图或导出3. wrfout 变量名对照表与物理量转换别再靠猜用 getvar() 别名直取真实物理量WRF 输出变量命名规则混乱T2是 2m 温度T是三维位温TK是三维气温TC是 2m 摄氏度温度RH是相对湿度QVAPOR是比湿……但 wrf-python 的getvar()接口提供了统一别名层屏蔽了 WRF 版本差异。下表列出最常用 15 个物理量的getvar()别名、原始 NetCDF 变量名、单位、维度及备注物理量getvar()别名常见原始变量名单位维度备注2m 温度T2T2K(time, south_north, west_east)直接可用无需插值10m 风速uv10U10,V10m/s(time, south_north, west_east)返回(u10, v10)元组自动旋转至地理坐标850 hPa 温度tempT,P,PBK(time, south_north, west_east)interplevel()封装别名temp自动调用插值地表气压PSFCPSFCPa(time, south_north, west_east)单位是帕斯卡非 hPa需/100转换相对湿度rhRH%(time, bottom_top, south_north, west_east)WRF 4.0 默认输出旧版需wrf.rh计算比湿qQVAPORkg/kg(time, bottom_top, south_north, west_east)注意QVAPOR是水汽混合比非比湿但 wrf-python 自动转换垂直速度waWm/s(time, bottom_top, south_north, west_east)W是 eta 坐标下的垂直速度wa转换为 z 坐标雷达反射率dbzDBZ8,DBZ16dBZ(time, bottom_top, south_north, west_east)需开启radar诊断选项DBZ8为 8-class 方案云水混合比QCLOUDQCLOUDkg/kg(time, bottom_top, south_north, west_east)仅 microphysics 方案启用时存在地表感热通量HFXHFXW/m²(time, south_north, west_east)正值表示地表向大气输送热量降水率RAINNCRAINNC,RAINCmm(time, south_north, west_east)RAINNC为层云降水RAINC为对流降水总降水 RAINNC RAINC涡度avoAVO1/s(time, bottom_top, south_north, west_east)绝对涡度含地球自转项散度divDIV1/s(time, bottom_top, south_north, west_east)水平散度位势高度zPHB,PHm(time, bottom_top1, south_north, west_east)z返回各 eta 层中点高度m非 geopotential height水汽通量water_vapor_fluxU,V,QVAPORkg/(m·s)(time, bottom_top, south_north, west_east)自动计算u*q,v*q并垂直积分注意getvar()别名是 wrf-python 的“语义层”它会根据 WRF 版本自动选择底层变量。例如wrf.getvar(ncfile, rh)在 WRF 4.2 中读RH在 WRF 3.9 中则调用wrf.rh函数基于T,P,QVAPOR计算——你无需关心版本兼容性。4. wrfout 提取的 5 个高频避坑点从坐标错位到单位翻车全是实测血泪经验4.1 现象画出的温度场明显偏移 100 km经纬度网格与底图不重合原因WRF 使用 Lambert Conformal 投影但wrf.getvar(ncfile, lat)返回的是球面经纬度XLAT,XLONG而 Cartopy 默认用 PlateCarree 投影叠加——两者数学定义不同直接叠加必然错位。解决用wrf.get_cartopy()获取 WRF 投影对象传给 Cartopy 的projection参数import cartopy.crs as ccrs import matplotlib.pyplot as plt cart_proj wrf.get_cartopy(ncfile) # 自动解析 MAP_PROJ 等参数 fig, ax plt.subplots(subplot_kw{projection: cart_proj}) ax.contourf(lons, lats, t_850[0], transformccrs.PlateCarree()) # 错transform 错 # 正确写法 ax.contourf(lons, lats, t_850[0]) # 不加 transform因 ax 已设为 WRF 投影4.2 现象wrf.getvar(ncfile, T2)返回值全是 0 或 NaN原因T2是 2m 温度但 WRF 输出控制中未开启sf_surface_physics选项sf_sfclay_physics 1或T2变量未写入 netCDF检查ncdump -h wrfout | grep T2。解决先确认变量存在再检查 WRFnamelist.input中sf_surf_layer和sf_surface是否启用若变量缺失需重跑 WRF 并设置history_outname wrfout_ddomain_date_time且io_form_history 2netCDF。4.3 现象interplevel()报错ValueError: vertical coordinate must be monotonic原因pressure场在某格点上非单调如地形起伏剧烈处低层气压反而高于高层scipy.interpolate.interp1d拒绝插值。解决用wrf.interplevel(..., missingnp.nan, keep_coordsTrue)并捕获异常或改用wrf.vinterp()支持多层插值# 更鲁棒的写法 try: t_850 wrf.interplevel(t, pressure, 85000, missingnp.nan) except ValueError: # 回退用 vinterp 对每个格点单独插值 t_850 wrf.vinterp(ncfile, fieldt, vert_coordpressure, interp_levels[85000], extrapolateTrue, missingnp.nan)4.4 现象PSFC单位是 Pa但业务系统要求 hPa手动/100后发现小数点后全为 0原因PSFC是整型压缩存储short类型scale_factor0.01add_offset0NetCDF 库自动解压但wrf.getvar()返回的是 float64已还原真实值——你再/100就错了。解决直接用psfc / 100得 hPa但确认psfc.dtype np.float64若不确定查ncfile.variables[PSFC].scale_factor。4.5 现象提取uvmet后风矢量方向全反东北风画成西南风原因uvmet是旋转至地理坐标的 U/V 分量但wrf.getvar(ncfile, uvmet)返回(u, v)元组而plt.quiver()默认u是东向分量、v是北向分量——这没错错在quiver()的X,Y网格顺序若lons是(south_north, west_east)lats同理则quiver(lons, lats, u, v)会错位。解决确保X,Y与U,V维度严格一致且X是经度网格west_east 维度在最后Y是纬度网格south_north 维度在最后# 正确lons, lats 均为 (south_north, west_east) u, v wrf.getvar(ncfile, uvmet, timeidx0) plt.quiver(lons, lats, u, v) # OK # 错误若 lons 是 (west_east, south_north)需 transpose plt.quiver(lons.T, lats.T, u, v) # 才对5. 进阶技巧批量处理 wrfout 文件、导出为 GeoTIFF、与机器学习 pipeline 对接当你的任务从“单次提取”升级为“每日自动处理 3 个域、7 天预报、每 3 小时一个文件”硬编码路径和timeidx就不可持续。wrf-python 本身不提供批量调度但可无缝接入 Python 生态的标准工具链。5.1 批量提取用 glob xarray 合并多文件时间序列WRF 常按时间切片输出多个wrfout_d01_YYYY-MM-DD_HH:MM:SS文件。传统做法是循环Dataset()但xarray.open_mfdataset()可自动合并import glob import xarray as xr # 匹配所有 wrfout 文件按时间排序 files sorted(glob.glob(/data/wrfout_d01_2023-06-15*)) # 用 wrf-python 的 backend 读取需安装 xarray-wrf ds xr.open_mfdataset(files, enginewrf, combineby_coords, concat_dimTime, data_varsminimal, coordsminimal) # 提取 2m 温度时间序列自动拼接 t2_series wrf.getvar(ds, T2) # shape: (Time, south_north, west_east)注意enginewrf是 xarray-wrf 插件提供的后端需pip install xarray-wrf它复用了 wrf-python 的坐标解析逻辑但返回xarray.Dataset支持.to_netcdf()、.sel()、.groupby()等高级操作。5.2 导出为 GeoTIFF让气象数据进入 GIS 工作流rasterio和rioxarray可将 wrf-python 提取的数组转为带地理参考的 GeoTIFF供 QGIS/ArcGIS 直接加载import rioxarray import numpy as np # 提取单时刻 850 hPa 温度 t_850 wrf.interplevel(t, pressure, 85000)[0] # (south_north, west_east) # 转为 xarray DataArray 并添加 CRS da xr.DataArray( t_850, dims(y, x), coords{y: lats[0, :, 0], x: lons[0, 0, :]} # 用第一列/行作为坐标轴 ).rio.write_crs(epsg:4326) # WGS84 # 写入 GeoTIFF自动插值为矩形网格 da.rio.to_raster(t850_2023061512.tif, dtypefloat32, compresslzw)5.3 与机器学习 pipeline 对接用 wrf-python 输出标准化张量深度学习模型如 U-Net 预测降水需要(batch, channel, height, width)格式张量。wrf-python 提取结果可直接torch.tensor()import torch # 提取 3 个变量T2, RH, PSFC组成通道 t2 wrf.getvar(ncfile, T2, timeidx0) # (south_north, west_east) rh wrf.getvar(ncfile, rh, timeidx0) # (south_north, west_east) psfc wrf.getvar(ncfile, PSFC, timeidx0) # (south_north, west_east) # 归一化按业务需求 t2_norm (t2 - 273.15 - 15) / 30 # 摄氏度中心化缩放 rh_norm rh / 100.0 # 0~1 psfc_norm (psfc - 95000) / 10000 # hPa 量级 # 合并为 (3, H, W) 张量 tensor_input torch.stack([ torch.from_numpy(t2_norm), torch.from_numpy(rh_norm), torch.from_numpy(psfc_norm) ], dim0) # 输入模型 # model(tensor_input.unsqueeze(0)) # 添加 batch 维度我坚持把 wrf-python 当作“WRF 数据的 ORM 层”——就像 Django ORM 屏蔽 SQL 差异一样它让我专注物理量逻辑而不是ncfile.variables[QVAPOR].scale_factor。过去三年我用这套流程处理了超 200 TB 的 WRF 输出零次因变量名或坐标系翻车。如果你还在用 NCL 脚本硬编码do i1,100; do j1,100; ...现在就是切换的临界点。希望帮到你。本文还有配套的精品资源点击获取
