ValueError with mpcalc.parcel_profile

[dopplershift]

Discussed in https://github.com/Unidata/MetPy/discussions/2308

^{Originally posted by **ejt4m8** January 25, 2022} Hello! After the release of MetPy 1.2 I wanted to go back and troubleshoot an old headache with mpcalc.parcel_path. I was having the same issues that others were having with repeated values of pressure which I saw was fixed in this update. I'm working with data from a well-known sounding, May 22, 2011 at 12z from SGF, the Springfield, MO sounding on the day of the Joplin tornado. I pull the data using siphon and assign to variables with the following code: ```python from datetime import datetime from siphon.simplewebservice.wyoming import WyomingUpperAir import metpy.calc as mpcalc from metpy.units import units date = datetime(2011, 5, 22, 12) station = 'SGF' df = WyomingUpperAir.request_data(date, station) p = df['pressure'].values * units(df.units['pressure']) T = df['temperature'].values * units(df.units['temperature']) Td = df['dewpoint'].values * units(df.units['dewpoint']) u = df['u_wind'].values * units(df.units['u_wind']) v = df['v_wind'].values * units(df.units['v_wind']) path = mpcalc.parcel_profile(p, T[0], Td[0]) ``` That code yields the following error message: ```pytb ValueError Traceback (most recent call last) Input In [15], in ----> 1 parcel_path = mpcalc.parcel_profile(p, T[0], Td[0]) File ~\anaconda3\envs\metpy\lib\site-packages\metpy\xarray.py:1230, in preprocess_and_wrap..decorator..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 ~\anaconda3\envs\metpy\lib\site-packages\metpy\units.py:288, in check_units..dec..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 ~\anaconda3\envs\metpy\lib\site-packages\metpy\calc\thermo.py:761, in parcel_profile(pressure, temperature, dewpoint) 722 @exporter.export 723 @preprocess_and_wrap(wrap_like='pressure') 724 @check_units('[pressure]', '[temperature]', '[temperature]') 725 def parcel_profile(pressure, temperature, dewpoint): 726 r"""Calculate the profile a parcel takes through the atmosphere. 727 728 The parcel starts at `temperature`, and `dewpoint`, lifted up (...) 759 760 """ --> 761 _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpoint) 762 return concatenate((t_l, t_u)) File ~\anaconda3\envs\metpy\lib\site-packages\metpy\calc\thermo.py:939, in _parcel_profile_helper(pressure, temperature, dewpoint) 937 # Find moist pseudo-adiabatic profile starting at the LCL 938 press_upper = concatenate((press_lcl, pressure[pressure < press_lcl])) --> 939 temp_upper = moist_lapse(press_upper, temp_lower[-1]).to(temp_lower.units) 941 # Return profile pieces 942 return (press_lower[:-1], press_lcl, press_upper[1:], 943 temp_lower[:-1], temp_lcl, temp_upper[1:]) File ~\anaconda3\envs\metpy\lib\site-packages\metpy\xarray.py:1230, in preprocess_and_wrap..decorator..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 ~\anaconda3\envs\metpy\lib\site-packages\metpy\units.py:344, in process_units..dec..wrapper(*args, **kwargs) 341 _mutate_arguments(bound_args, units.Quantity, lambda val, _: val.to_base_units().m) 343 # Evaluate inner calculation --> 344 result = func(*bound_args.args, **bound_args.kwargs) 346 # Wrap output 347 if multiple_output: File ~\anaconda3\envs\metpy\lib\site-packages\metpy\calc\thermo.py:345, in moist_lapse(pressure, temperature, reference_pressure) 342 press_side = pressure[points_above][::-1] 344 # Flip on exit so t values correspond to increasing pressure --> 345 trace = si.solve_ivp(t_span=(reference_pressure, press_side[-1]), 346 t_eval=press_side, **solver_args).y[..., ::-1] 347 ret = np.concatenate((trace, ret), axis=-1) 349 # Do we have any points below the reference pressure File ~\anaconda3\envs\metpy\lib\site-packages\scipy\integrate\_ivp\ivp.py:529, in solve_ivp(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, args, **options) 527 d = np.diff(t_eval) 528 if tf > t0 and np.any(d <= 0) or tf < t0 and np.any(d >= 0): --> 529 raise ValueError("Values in `t_eval` are not properly sorted.") 531 if tf > t0: 532 t_eval_i = 0 ValueError: Values in `t_eval` are not properly sorted. ``` This error only occurs with some dates, times and stations. For other times and stations the parcel path calculates like normal.

[kgoebber]

Okay, I just ran into this as well for the following date:

date = datetime(2007, 5, 4, 18)

In inspecting the data, it appears that there were two pressure levels that were the same. In this instance it was all the way up at 24.4 hPa. If I subset to use only levels where the pressure > 90, it runs just fine.

Checking the example provided above - there the culprit is 27.7 hPa (and would also error at 19.0 hPa).

So to overcome this error my data parsing would look like:

from datetime import datetime

import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.plots import SkewT
from metpy.units import units, pandas_dataframe_to_unit_arrays
import numpy as np
from siphon.simplewebservice.wyoming import WyomingUpperAir

# Set date that you want
# Data goes back to the 1970's
date = datetime(2007, 5, 4, 18)

# Set station ID, there are different stations back in the day
# Current station IDs found at http://weather.rap.ucar.edu/upper
station = 'DDC'

# Use Siphon module to grab data from remote server
df = WyomingUpperAir.request_data(date, station)
# Create dictionary of unit arrays
data = pandas_dataframe_to_unit_arrays(df)

# Isolate united arrays from dictionary to individual variables
p = data['pressure']
idx = p.m > 90
p = p[idx]
T = data['temperature'][idx]
Td = data['dewpoint'][idx]
u = data['u_wind'][idx]
v = data['v_wind'][idx]

So at the very least we should check the pressures to see if there are any repeats and error on that to give the suggestion of subsetting the data.

The question is whether we want to bake in any subsetting by default to limit this from happening for users?

In general, this is likely to happen only at very high altitudes (low pressure) since the change of pressure is less rapid at those levels. But I also wonder if this will be more frequent (and happen at higher pressure) with the increased amount of data coming from modern data files.

[ktyle]

@kgoebber Another very recent example is the KLIX 3/23/2022 0000 UTC sounding (just prior to the EF-3 in New Orleans). Similar to your example (and others I'm aware of) it's also an instance of adjacent levels with the same pressure, well up in the stratosphere (in this case, 17.4 hPa).

The workaround I've employed in an example notebooki is to call pandas' drop_duplicates before calling parcel_profile as follows:

df.drop_duplicates(inplace=True,subset='pressure',ignore_index=True)

Perhaps this could be integrated into the parcel_profile function?