1

When attempting to use the find the wet-bulb temperature for a parcel returned by MetPy's most_unstable_parcel function, I am getting an error that seems to be related to the Pint attributes of the most unstable parcel pressure. It is necessary to strip away those attributes by coercing the pressure to be of type int or float.

See discussion, below.

I would greatly appreciate some help in understanding how to use the most unstable parcel pressure without resorting to the artifice of coercing it to a type that requires re-assigning units. Thank you!

Example Calculation (data from OUN 20220328 00Z sounding):

from metpy.units import units
import metpy.calc as mpcalc

#  Construct a sample sounding:
p = [972.0, 942.8, 925.0, 910.1, 878.2, 850.0, 847.3, 826.0, 817.0, 787.5]*units('hPa')
T = [25.4, 22.4, 20.6, 19.3, 16.4, 13.8, 13.6, 11.6, 11.1, 9.4]*units('degC')
Td = [6.4, 5.3, 4.6, 4.1, 2.9, 1.8, 1.7, 0.6, -0.2, -2.9]*units('degC')  
h = [345,  610,  776,  914, 1219, 1497, 1524, 1737, 1829, 2134]*units('m')

#  Find the Most Unstable Parcel
MUp = mpcalc.most_unstable_parcel(p,T,Td,h)

# MUp
# (972.0 <Unit('hectopascal')>,
#  25.4 <Unit('degree_Celsius')>,
#  6.4 <Unit('degree_Celsius')>,
#  0)

#  Next compute the wet-bulb-temperature for this parcel:

Twb = mpcalc.wet_bulb_temperature(MUp[0],MUp[1],MUp[2])

When this is executed, we get the following error traceback:

IndexError                                Traceback (most recent call last)
Input In [4], in <module>
----> 1 Twb = mpcalc.wet_bulb_temperature(MUp[0],MUp[1],MUp[2])

File ~/miniconda3/envs/MetPy/lib/python3.8/site-packages/metpy/xarray.py:1230, in preprocess_and_wrap.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
   1227     _mutate_arguments(bound_args, units.Quantity, lambda arg, _: arg.m)
   1229 # Evaluate inner calculation
-> 1230 result = func(*bound_args.args, **bound_args.kwargs)
   1232 # Wrap output based on match and match_unit
   1233 if match is None:

File ~/miniconda3/envs/MetPy/lib/python3.8/site-packages/metpy/units.py:288, in check_units.<locals>.dec.<locals>.wrapper(*args, **kwargs)
    285 @functools.wraps(func)
    286 def wrapper(*args, **kwargs):
    287     _check_units_inner_helper(func, sig, defaults, dims, *args, **kwargs)
--> 288     return func(*args, **kwargs)

File ~/miniconda3/envs/MetPy/lib/python3.8/site-packages/metpy/calc/thermo.py:3073, in wet_bulb_temperature(pressure, temperature, dewpoint)
   3071 ret = it.operands[3]
   3072 if ret.size == 1:
-> 3073     ret = ret[0]
   3074 return units.Quantity(ret, temperature.units)

IndexError: too many indices for array: array is 0-dimensional, but 1 were indexed

After running a few permutations of inputs, I discovered that it is possible to use the output pressure value if one first coerces it to, say, type float or int.

#  Reformat the pressure value as input to
#  wet_bulb_temperature()

d = float(MUp[0].magnitude)*units('hPa')
Twb1 = mpcalc.wet_bulb_temperature(d,MUp[1],MUp[2])
print(Twb1)

#  output is 14.209821596989343 degree_Celsius

# As far as Python is concerned, the two instances of pressure are
# identical:


print(MUp[0]==d)

#  Output:  True

It's just not clear to me why the pressure output from most_unstable_parcel (e.g. MUp[0]) cannot be used directly as an input to the wet-bulb temperature calculation.

gdlewen
  • 73
  • 4
  • can you describe some background of your goal and the modules you use? – Lei Yang Mar 29 '22 at 00:45
  • I am trying to implement the Lid Strength Index as described in Keller, David L. "P7. 1 A COMPARATIVE VERIFICATION OF TWO “CAP” INDICES IN FORECASTING THUNDERSTORMS." The only modules I am currently using are the imports given as part of the example in the question and the MetPy Plot SkewT function. – gdlewen Mar 29 '22 at 00:52
  • Thank you, but I'm sorry I do not understand the question, "by *units... what do you expect?" – gdlewen Mar 29 '22 at 01:00
  • I am taking the array of values and assigning units of hPa. p will now be an object with magnitude (the values in the array) and units of hPa. Most MetPy calculations require the inputs to have units associated with them. If you run the example you will see that the two temperatures returned by most_unstable_parcel() have units of ˚C assigned to them, but I can coerce those to Kelvin if I choose. – gdlewen Mar 29 '22 at 01:04
  • Let us [continue this discussion in chat](https://chat.stackoverflow.com/rooms/243398/discussion-between-gdlewen-and-lei-yang). – gdlewen Mar 29 '22 at 01:09

1 Answers1

3

This is a bug with MetPy's wet_bulb_temperature function in handling NumPy scalars, so thank you for identifying the issue! A pull request fixing this has been submitted, which should make it into the upcoming MetPy v1.3.0 release.

For now, you will need to somehow override the array shape recognition that exists in wet_bulb_temperature. Converting from a NumPy scalar to a Python scalar as you discovered is one such way to do so. Another is converting to 1D NumPy arrays, like the following:

Twb = mpcalc.wet_bulb_temperature(
    *(np.atleast_1d(arg) for arg in (MUp[0], MUp[1], MUp[2]))
)
Jon Thielen
  • 479
  • 2
  • 7