企业官网建设流程全解析

简介：本资源面向气象建模初学者与WRF用户，提供一套轻量级、即用型的wrfout文件变量提取与可视化实践方案，解决WRF模拟输出数据难以直接读取和分析的常见痛点。压缩包仅2KB，含2个核心脚本：t2m.py（Python实现，基于netCDF4/xarray提取地表温度T2M并绘图）与wrfout.ncl（NCL实现，调用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 git+https://github.com/NCAR/wrf-python.git@main

提示：conda install会自动拉取netcdf4=1.6.4,pyproj=3.9.0,numpy=1.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.1）

2.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", timeidx=wrf.ALL_TIMES) # shape: (time, south_north, west_east) lons = wrf.getvar(ncfile, "lon", timeidx=wrf.ALL_TIMES) # 3. 提取地表气压（单位：Pa，自动去单位转换） psfc = wrf.getvar(ncfile, "PSFC", timeidx=wrf.ALL_TIMES) # shape: (time, south_north, west_east)

关键参数说明：

• timeidx=wrf.ALL_TIMES：提取所有时间步（若只取第 0 个时刻，用timeidx=0）
• wrf.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 层权重、边界外推、多时间步广播：

# 提取三维位温 T（shape: time, bottom_top, south_north, west_east） t = wrf.getvar(ncfile, "T", timeidx=wrf.ALL_TIMES) # 单位：K # 提取三维气压 P（单位：Pa） p = wrf.getvar(ncfile, "P", timeidx=wrf.ALL_TIMES) # 修正气压（perturbation） pb = wrf.getvar(ncfile, "PB", timeidx=wrf.ALL_TIMES) # 基准气压（base state） pressure = p + pb # 总气压，单位 Pa # 插值得到 850 hPa 层上的位温（注意：850 hPa = 85000 Pa） t_850 = wrf.interplevel(t, pressure, 85000, missing=np.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一致），missing=np.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()别名是 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], transform=ccrs.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_d<domain>_<date>_<time>'且io_form_history = 2（netCDF）。

4.3 现象：interplevel()报错ValueError: vertical coordinate must be monotonic

原因：pressure场在某格点上非单调（如地形起伏剧烈处，低层气压反而高于高层），scipy.interpolate.interp1d拒绝插值。
解决：用wrf.interplevel(..., missing=np.nan, keep_coords=True)并捕获异常，或改用wrf.vinterp()（支持多层插值）：

# 更鲁棒的写法 try: t_850 = wrf.interplevel(t, pressure, 85000, missing=np.nan) except ValueError: # 回退：用 vinterp 对每个格点单独插值 t_850 = wrf.vinterp(ncfile, field=t, vert_coord="pressure", interp_levels=[85000], extrapolate=True, missing=np.nan)

4.4 现象："PSFC"单位是 Pa，但业务系统要求 hPa，手动/100后发现小数点后全为 0

原因：PSFC是整型压缩存储（short类型），scale_factor=0.01，add_offset=0，NetCDF 库自动解压，但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", timeidx=0) 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, engine="wrf", combine="by_coords", concat_dim="Time", data_vars="minimal", coords="minimal") # 提取 2m 温度时间序列（自动拼接） t2_series = wrf.getvar(ds, "T2") # shape: (Time, south_north, west_east)

注意：engine="wrf"是 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", dtype="float32", compress="lzw")

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", timeidx=0) # (south_north, west_east) rh = wrf.getvar(ncfile, "rh", timeidx=0) # (south_north, west_east) psfc = wrf.getvar(ncfile, "PSFC", timeidx=0) # (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) ], dim=0) # 输入模型 # model(tensor_input.unsqueeze(0)) # 添加 batch 维度

我坚持把 wrf-python 当作“WRF 数据的 ORM 层”——就像 Django ORM 屏蔽 SQL 差异一样，它让我专注物理量逻辑，而不是ncfile.variables['QVAPOR'].scale_factor。过去三年，我用这套流程处理了超 200 TB 的 WRF 输出，零次因变量名或坐标系翻车。如果你还在用 NCL 脚本硬编码do i=1,100; do j=1,100; ...，现在就是切换的临界点。希望帮到你。

本文还有配套的精品资源，点击获取