使用WRF数据绘制Metpy SkewT时出现量纲错误的解决方法

问题场景 #

我正在使用WRF输出数据绘制SkewT，编写的Python代码如下：

import wrf
from netCDF4 import Dataset
import matplotlib.pyplot as plt
import numpy as np
import metpy.calc as mpcalc
from metpy.plots import SkewT
from metpy.units import units

wrfin = Dataset(r'wrfout_d02_2022-06-19_00_00_00')
lat_lon = [25.0803, 121.2183]
x_y = wrf.ll_to_xy(wrfin, lat_lon[0], lat_lon[1])

p1 = wrf.getvar(wrfin,"pressure",timeidx=0)
T1 = wrf.getvar(wrfin,"tc",timeidx=0)
Td1 = wrf.getvar(wrfin,"td",timeidx=0)
u1 = wrf.getvar(wrfin,"ua",timeidx=0)
v1 = wrf.getvar(wrfin,"va",timeidx=0)

p = p1[:,x_y[0],x_y[1]] * units.hPa
T = T1[:,x_y[0],x_y[1]] * units.degC
Td = Td1[:,x_y[0],x_y[1]] * units.degC
u = v1[:,x_y[0],x_y[1]] * units('m/s')
v = u1[:,x_y[0],x_y[1]] * units('m/s')

skew = SkewT()
skew.plot(p, T, 'r')
skew.plot(p, Td, 'g')

my_interval =  np.arange(100, 1000, 50) * units('mbar')
ix = mpcalc.resample_nn_1d(p, my_interval)
skew.plot_barbs(p[ix], u[ix], v[ix])

skew.plot_dry_adiabats()
skew.plot_moist_adiabats()
skew.plot_mixing_lines()
skew.ax.set_ylim(1000, 100)
skew.ax.set_xlim(-60, 40)
skew.ax.set_xlabel('Temperature ($^\circ$C)')
skew.ax.set_ylabel('Pressure (hPa)')

plt.savefig('SkewT.png', bbox_inches='tight')

运行代码时出现错误：

raise DimensionalityError(
pint.errors.DimensionalityError: Cannot convert from 'dimensionless' (dimensionless) to 'millibar' ([mass] / [length] / [time] ** 2)

当前环境：Python 3.9.7，Metpy 0.12.0

问题原因及解决方案 #

1. 旧版本Metpy的单位处理局限性 #

Metpy 0.12.0是较早期版本，对Pint单位的支持不完善，resample_nn_1d函数无法正确识别等价单位（如hPa与mbar），或要求输入为无单位的数值数组。

方案一：使用无单位数值数组

提取压力变量和重采样区间的纯数值（去掉单位）后传入函数：

# 修改重采样部分代码
p_vals = p.magnitude
my_interval = np.arange(100, 1000, 50)  # 数值对应hPa单位
ix = mpcalc.resample_nn_1d(p_vals, my_interval)
skew.plot_barbs(p[ix], u[ix], v[ix])

方案二：统一单位表述

将重采样区间的单位改为与p一致的hPa，避免单位识别冲突：

my_interval = np.arange(100, 1000, 50) * units.hPa
ix = mpcalc.resample_nn_1d(p, my_interval)
skew.plot_barbs(p[ix], u[ix], v[ix])

2. 修正风场变量赋值错误 #

代码中u和v分量的赋值完全搞反：WRF的ua对应东向风（u分量），va对应北向风（v分量），正确赋值应为：

u = u1[:,x_y[0],x_y[1]] * units('m/s')
v = v1[:,x_y[0],x_y[1]] * units('m/s')

3. 可选：升级Metpy版本 #

如果环境允许，建议升级到Metpy 1.0+版本，新版本对单位处理、WRF数据兼容性有大幅提升，可从根源避免这类旧版本bug：

pip install --upgrade metpy

内容的提问来源于stack exchange，提问作者rock0789