pythonmetpy

Getting error message while ploting metpy SkewT


I am using WRF output data to plot SkweT, here is code:

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')

but running error message especially in the block:

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

It seems that mpcalc.resample_nn_1d doesn't work.

how can I solve it?

python version : 3.9.7

metpy version : 0.12.0


Solution

  • My guess is that wrf.getvar() is returning NumPy masked arrays, which behave a little weird with the Pint Quantity instances (compared with normal arrays). I would recommend trying this syntax to "attach" units:

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

    Also, a heads up that MetPy 0.12 is quite a bit out of date, and I would recommend updating to the latest version (currently 1.3.1) as soon as you can.