2026/10/20 20:33:56

WRF后处理神器wrf-python 1.1.3:诊断量获取、插值与坐标转换实战

WRF后处理神器wrf-python 1.1.3:诊断量获取、插值与坐标转换实战 简介wrf-python 是用于 WRF 天气研究与预报模式输出的 Python 处理与可视化接口库面向气象科研人员、大气科学专业学生及数值模式应用开发者可在 Python 环境中直接读取 WRF 网格数据完成变量诊断与地图投影绘图。该压缩包为 1.3MB共收录 328 个文件其中 51 个核心 Python 脚本构成库的主体功能18 个 Fortran 源文件提供底层重插值计算支持222 个 rst 文档覆盖完整的 API 使用说明与示例另有 Windows 平台构建脚本、png 图表与配置文件辅助用户理解编译安装流程。通过阅读 Python 脚本与 Fortran 内核可深入了解变量标准化、坐标投影与站点插值的实现细节便于在科研项目中定制诊断计算官方 rst 文档则能帮助使用者快速对接 MetPy、Cartopy 等常用气象绘图生态。目前已有 2408 人学习下载适合需要嵌入 WRF 后处理流程或自主查阅源码开展二次开发的 Python 使用者。1. wrf-python 1.1.3 是什么WRF 输出诊断量的“一站式解算器”wrf-python 1.1.3 是 WRF 模式后处理绕不开的一个 Python 库。跑完 WRF 你会得到一堆 NetCDF 格式的 wrfout 文件最头疼的不是画图而是把模式内部的扰动位温、扰动气压还原成有物理意义的温度、风场、海平面气压。wrf-python 做的事情就是把 WRF/NCL 里繁琐的一整套诊断公式封装成 getvar、vinterp、xy_to_ll 这些函数你只需要传入 wrfout 文件路径和变量名它自动处理坐标投影、单位换算、基态叠加。它解决的是“模式输出到科研分析”中间那段最枯燥的脏活。这个库适合两类人一类是刚跑完 WRF 正在想办法从 nc 文件里抠变量的新手另一类是天天做个例分析、画了无数张剖面图的熟手。前者省去翻公式的时间后者省去重复造轮子的时间。往下我从安装编译、核心 API、常见坑位和批量验证四个层面把这个包拆给你看。2. 从 tar.gz 装起源码编译需要什么、卡在哪2.1 环境准备先解决 gdal 和 netCDF4 的系统依赖wrf-python 的底层不是纯 PythonCython 写的核心算法要调 GDAL 做投影变换、调 netCDF4 读模式文件。所以拿到 tar.gz 直接pip install大概率会翻车报错多半是找不到 gdal-config 或者cannot find -lgdal。这里有个关键认知wrf-python 是依赖系统级地理库的不是普通纯 Python 包装它的顺序应该是先把 GDAL、netCDF4、numpy 这些底座装好再编译 wrf-python 本体。我一般用 conda 建独立环境把版本隔离干净。你可以按下面这套来准备命令里的 python3.6 不是随便选的wrf-python 1.1.3 发布时主流环境是 Python 3.5/3.6在 Python 3.8 上编译 Cython 生成的代码偶尔会报函数签名错误用 3.6 最稳。conda create -n wrfpy python3.6 -y conda activate wrfpy conda install -c conda-forge gdal netcdf4 numpy wrapt -y这段命令做了三件事创建名为 wrfpy 的独立环境避免污染机器上的其他 Python 项目把 gdal、netcdf4、numpy、wrapt 装进这个环境gdal 从 conda-forge 通道装因为 conda 会连带把 proj、geos 这些底层库一起装好pip 装 gdal 经常会缺头文件或者链接不上。wrapt 是 wrf-python 运行时用到的装饰器库少了它 import 阶段就会出错。依赖表大致是下面这个样子装之前心里有个数依赖包用途推荐来源numpy数组运算基础conda-forgenetCDF4读取 wrfout 文件conda-forgegdal投影变换和地理坐标计算conda-forgewrapt运行时方法包装conda-forgematplotlib绘图可选conda-forgematplotlib 只有在用 WRF 自带绘图函数时才需要如果只做数值计算不画图可以不装。basemap 在 1.1.3 时代是默认绘图后端但现在基本被 cartopy 替代我建议不要为了画图去折腾 basemap单独计算物理量用不到它。2.2 真实的安装流程从解包到 import 成功把 wrf-python-1.1.3.tar.gz 解压到工作目录然后进入目录执行编译。很多人在这一步直接python setup.py install在 conda 环境下大概率还是找不到 gdal 的头文件路径因为 conda 的 include 和 lib 目录不在系统默认路径里。你需要手动告诉编译器 gdal 装在哪里。tar zxvf wrf-python-1.1.3.tar.gz cd wrf-python-1.1.3 python setup.py build_ext \ --include-dirs$CONDA_PREFIX/include \ --library-dirs$CONDA_PREFIX/lib python setup.py install这里解释一下两个关键参数--include-dirs指向 conda 环境的 include 目录编译 gdal 相关的 C 扩展时需要找到 gdal_priv.h 头文件--library-dirs指向 lib 目录链接阶段要在这里找 libgdal.so 和 libnetcdf.so。如果你的 conda 环境路径不一样把$CONDA_PREFIX换成实际路径即可。这一步成功之后setup.py install 会把 wrf 包写进 site-packages。还有一种更省事的做法如果 conda-forge 通道已经编好了对应版本直接装二进制包conda install -c conda-forge wrf-python1.1.3 -y我个人推荐优先试 conda 方式装完就能用不用碰编译器。但如果你需要看源码或者改底层tar.gz 方式是绕不开的。装完后验证一下python -c import wrf; print(wrf.__version__)能看到1.1.3输出说明编译链接全部正常环境算是立住了。2.3 依赖版本边界1.1.3 与 numpy、Python 3.8 的兼容wrf-python 1.1.3 的适配边界比新版本窄很多翻车点集中在两个方面。第一是 Python 版本前面说过 3.6 是稳妥选择3.8 以上编译时如果碰到_wrf扩展模块的报错不要怀疑自己的操作直接换成 Python 3.6 环境重来。第二是 numpy 版本我在 numpy 1.15 和 1.19 上都跑通过再往上的版本偶尔会出现二进制接口不匹配症状是 import 时直接段错误。这里有个容易忽略的点如果 conda 环境里的 netCDF4 和 numpy 是从不同通道混装的二进制 ABI 对不上wrf 在读取数据时也可能崩溃。解决的办法不是去逐个调试而是建新环境时统一指定 conda-forge 通道让这些包使用同一套编译工具链编出来的二进制问题基本就消失了。如果编译过程中报gdal-config not found说明 gdal 虽然装了但环境变量 PATH 里没有 gdal-config检查一下which gdal-config。conda-forge 的 gdal 包正常会把 gdal-config 放进$CONDA_PREFIX/bin如果找不到重新执行conda install -c conda-forge gdal会把它补全。3. 三个核心函数getvar、vinterp 和坐标转换3.1 getvar 取诊断量变量名不是乱给的wrf-python 最大的价值体现在 getvar 函数。WRF 的 wrfout 文件里存储的是扰动量和基态量分离的原始变量比如气压拆成 P 和 PB位温是 T 和 P 的组合如果直接从 NetCDF 里读p再自己加pb也能算但要记住每个变量的单位、叠加方式、是否要做地形追随坐标换算非常琐碎。getvar 把这一层全部封装掉你只需要告诉它“我要什么物理量”它自己知道去读哪些原始变量、怎么组合。先看一段最基础的调用import wrf from netCDF4 import Dataset wrfout Dataset(wrfout_d01_2019-08-10_12:00:00.nc) ua wrf.getvar(wrfout, ua, timeidx0) slp wrf.getvar(wrfout, slp, timeidx0) print(type(ua)) print(ua.units) print(ua.shape)timeidx0表示取文件里第一个时间步wrfout 多时间步文件里这个参数控制取哪一时刻。返回的 ua 不是普通 numpy 数组而是带投影信息的 Variable 对象它继承了数组的全部行为可以直接做切片和算术运算同时额外携带了坐标投影元数据后面做 xy_to_ll 或者画图时这些元数据会被自动使用。print 出的 units 是m s-1shape 是三维网格维度单位已经帮你换算好了。常用变量名对照表如下变量名物理含义单位维度ua / va纬向/经向风m/s3Dwa垂直速度m/s3Dtk温度K3Dtc摄氏温度degC3Dtd露点温度degC3Dslp海平面气压hPa2Dz位势高度m3Dpressure全气压hPa3D注意没有temperature这种命名wrf-python 沿用 WRF 后处理的历史命名习惯。如果你传了不存在的变量名它会抛 KeyError 并列出可用变量列表跟着提示改就行。3.2 vinterp 把模式面插到气压面参数决定成败WRF 输出是地形追随坐标垂直方向不是等压面做天气图分析时通常需要把变量插值到 1000 hPa、850 hPa、500 hPa 这些标准气压层。vinterp 就是干这个的。它沿着每个水平格点做垂直插值输入一个三维诊断量场和三维气压场输出指定气压层上的变量值。import wrf from netCDF4 import Dataset wrfout Dataset(wrfout_d01_2019-08-10_12:00:00.nc) tk wrf.getvar(wrfout, tk) pressure wrf.getvar(wrfout, pressure) levels [1000., 925., 850., 700., 500., 300., 200., 100.] tk_pl wrf.vinterp(tk, pressure, levels, field_typetk, extrapolateTrue, missing1e20) print(tk_pl.shape)这段代码把温度从模式面插到 8 个标准等压面上。field_typetk告诉插值器这是温度场内部会按温度场的特点选择更合适的插值方式extrapolateTrue表示允许在模式层顶以上外推这个参数很关键后面避坑章节会详细说missing1e20是缺测值标记插值过程中如果某层数据不足这个值会被写入结果。vinterp 实现的原理是对每个格点沿垂直方向逐层扫描找到目标气压层所在的高度区间然后做线性或对数线性插值。返回数组的第一个维度对应你传入的 levels 列表长度后面两维是水平格点。使用时注意pressure 必须是 hPa 单位别把 wrfout 里的 Pa 原始变量直接喂进去wrf.getvar(wrfout, pressure) 返回的已经是 hPa。3.3 用 xy_to_ll 和 ll_to_xy 定位站点网格和经纬度互转做站点对比分析时手里常常是一个站点经纬度想知道它落在 WRF 网格的哪个位置或者反过来想知道某个格点对应的经纬度。坐标转换函数在这时候发挥作用。import wrf from netCDF4 import Dataset wrfout Dataset(wrfout_d01_2019-08-10_12:00:00.nc) lat0, lon0 35.2, 117.4 x, y wrf.ll_to_xy(wrfout, lat0, lon0) print(站点对应的网格浮点索引:, x, y) lat, lon wrf.xy_to_ll(wrfout, x, y) print(索引反算经纬度:, lat, lon)ll_to_xy返回 x 和 y 两个浮点值分别对应网格的东西向和南北向索引。注意 WRF 网格的 x 是列方向y 是行方向二维数组切片时用的是[y, x]的顺序写反了会气死人。返回的是浮点数是因为站点不一定正好落在格点上你可以四舍五入取最近格点也可以用它周边的四个格点做双线性插值。xy_to_ll是它的逆运算传入格点索引返回经纬度。坐标转换的底层依赖 wrfout 文件全局属性里的投影参数也就是 2.1 节为什么必须装 gdal 的原因。这套投影信息在 getvar 返回的变量里也被保存了所以从变量对象出发做转换通常更稳。4. 常见问题排查安装失败与计算异常的 4 个高发点4.1 编译时报 gdal-config not found现象执行python setup.py install或者pip install wrf-python1.1.3时终端输出gdal-config not found随后报错退出整个安装流程中断。原因wrf-python 的扩展模块编译时需要调用 gdal-config 程序来获取 gdal 的编译参数和链接参数。这个程序是 gdal 开发包的一部分只装了运行库或者根本没装 gdal自然找不到它。解决先通过 conda-forge 把 gdal 完整装好再确认 gdal-config 是否在环境变量 PATH 里。用which gdal-config查看如果没输出说明 conda 安装不完整或者环境激活有问题。重新执行conda activate wrfpy后再检查。如果 gdal-config 存在但仍然报这个错多半是权限或环境变量问题改用带绝对路径的方式编译python setup.py build_ext \ --include-dirs$CONDA_PREFIX/include \ --library-dirs$CONDA_PREFIX/lib4.2 import wrf 直接段错误现象辛辛苦苦装完import wrf一执行Python 解释器直接崩溃退出什么错误信息都没有只剩一行 Segmentation Fault。原因这是二进制接口不匹配的典型症状最常见的是 numpy 或 netCDF4 不是同一套工具链编译的。比如 conda 环境里装了 conda-forge 的 numpy又通过 pip 装了 pre-built 的 netCDF4两者底层 ABI 对不上C 扩展在加载时访问越界直接崩溃。解决不要在这个环境里盲目重装某单个包那是在赌运气。直接建一个全新环境所有包统一从 conda-forge 安装conda create -n wrfpy_fresh python3.6 -c conda-forge -y conda activate wrfpy_fresh conda install -c conda-forge gdal netcdf4 numpy wrapt -y conda install -c conda-forge wrf-python1.1.3 -y我遇到过的多起段错误全是混装引起统一用 conda-forge 后未再复发。4.3 getvar 返回结果全是 masked现象调用wrf.getvar(wrfout, ua)成功后打印结果发现所有数值都是 masked 状态或者全是 nan根本没法用。原因最常见的是直接用 netCDF4 的变量名去套 getvar比如用了U或P这种 wrfout 里的原始变量名wrf-python 内部按自己的元数据表去匹配匹配不到内容就会返回空值。另一个原因是 timeidx 超出文件时间步范围返回的空切片被整体标记为无效。解决用 wrf-python 文档里的标准变量名不要自己猜。写一个快速排查脚本先打印出文件里有几个时间步再确认变量名是否在支持列表里from netCDF4 import Dataset import wrf wrfout Dataset(wrfout_d01_2019-08-10_12:00:00.nc) print(wrfout.variables.keys()) print(时间步数:, wrfout.dimensions[Time].size) try: v wrf.getvar(wrfout, ua) print(读取成功shape:, v.shape) except KeyError as e: print(变量名不支持:, e)4.4 vinterp 顶层出现异常尖峰现象vinterp 插值到 200 hPa 或 100 hPa 后数值出现明显不合理的大尖峰比如温度一下子比周围高几十度或者在某个局部区域出现负的绝对温度。原因模式层顶高度有限如果目标气压层高于模式最顶层的实际气压插值程序只能在外推模式下用最外层的垂直梯度继续延拓。强对流区域垂直梯度本来就大梯度外推很容易放大误差形成尖峰。另外某些物理量在模式层顶附近本身就有不稳定性直接外推就是玄学。解决先检查 wrfout 的模式层顶气压值再决定 levels 是否合理。稳妥做法是关闭外推让超出范围的值变成缺测tk_pl wrf.vinterp(tk, pressure, levels, field_typetk, extrapolateFalse, missing-9999.0)这样超过模式顶的插值结果会被 fill 成 -9999不参与后续绘图和统计。如果你必须在高层用数据那就接受外推带来的不确定性至少做一次和探空或再分析资料的对比验证确认偏差在可控范围。5. 进阶用法批处理、站点对比与计算验证批量处理是实际项目中避免重复劳动的必经之路。WRF 一次跑几十个小时wrfout 文件多到按小时切分逐个手工操作根本不现实。用 glob 配合循环把整个时次序列的诊断量一次算完import glob from netCDF4 import Dataset import wrf files sorted(glob.glob(wrfout_d01_2019-08-10_??:00:00.nc)) for path in files: nc Dataset(path) slp wrf.getvar(nc, slp).values print(path, 海平面气压最小值:, slp.min()) nc.close()变量返回的对象在循环里要及时 close 文件避免文件句柄堆积。这里.values把带投影的 Variable 转成纯 numpy 数组如果后面不需要投影信息这一步能省不少内存。站点对比是验证 wrf-python 计算结果是否靠谱的硬指标。我习惯写一个函数把任意三维变量插到指定站点经纬度上得到该站的垂直廓线def extract_point_profile(ncfile, lat0, lon0, varname): var3d wrf.getvar(ncfile, varname) x, y wrf.ll_to_xy(ncfile, lat0, lon0) i int(round(x)) j int(round(y)) return var3d[:, j, i]注意切片顺序是先 j 后 i对应网格的行列这个顺序搞反了取到的就是错位格点。用这个函数把插值后的温度廓线和探测资料画在同一张图上偏差不超过 2 到 3 度基本可以信任。换了一台机器重新搭环境后不要急着出图先跑一遍这个验证流程。我最初带着别人的脚本换环境装完直接 batch 跑画出来的高空温度场比探空偏了 5 度排查了两天才发现是 vinterp 在顶层外推被顶穿了从那以后每次新环境我都会用同一个站点、同一个时次的廓线做回归校验确认没问题再批处理。这套小习惯看起来笨省下来的排查时间却是实打实的希望帮到你。本文还有配套的精品资源点击获取