kylejgillett / sounderpy

A python package that helps you to access and plot vertical profile data for meteorological analysis
https://kylejgillett.github.io/sounderpy/
MIT License
51 stars 13 forks source link

ValueError: Specified bound is outside pressure range. #31

Closed wreed1989 closed 5 hours ago

wreed1989 commented 4 days ago

Greetings! I'm working on a project for class where MetPy was a bit limited on what I could do with SkewTs so I was testing out SounderPy to see how if it could do what I wanted it to and I'm running into some issues with using a custom dataset. I've formatted the data accordingly based on the example from the Documents but I'm getting the "ValueError: Specified bound is outside pressure range." and I'm not real sure what is occurring to throw the error.

I will show my code followed by the error for troubleshooting. I'd appreciate the help!

# Convert 'Time [sec]' to int
Jul14_1700z['Time [sec]'] = pd.to_numeric(Jul14_1700z['Time [sec]'], errors='coerce')

# Convert 'Time (UTC)' to datetime
Jul14_1700z['Time (UTC)'] = pd.to_datetime(Jul14_1700z['Time (UTC)'], errors='coerce')

# Convert numerical columns (remove units and convert to float)
Jul14_1700z['P [hPa]'] = pd.to_numeric(Jul14_1700z['P [hPa]'], errors='coerce')
Jul14_1700z['T [°C]'] = pd.to_numeric(Jul14_1700z['T [°C]'], errors='coerce')
Jul14_1700z['Hu [%]'] = pd.to_numeric(Jul14_1700z['Hu [%]'], errors='coerce')
Jul14_1700z['Ws [m/s]'] = pd.to_numeric(Jul14_1700z['Ws [m/s]'], errors='coerce')
Jul14_1700z['Ws U [m/s]'] = pd.to_numeric(Jul14_1700z['Ws U [m/s]'], errors='coerce')
Jul14_1700z['Ws V [m/s]'] = pd.to_numeric(Jul14_1700z['Ws V [m/s]'], errors='coerce')
Jul14_1700z['Geopot [m]'] = pd.to_numeric(Jul14_1700z['Geopot [m]'], errors='coerce')
Jul14_1700z['Dew [°C]'] = pd.to_numeric(Jul14_1700z['Dew [°C]'], errors='coerce')

# Convert 'Wd [°]' to numeric
Jul14_1700z['Wd [°]'] = pd.to_numeric(Jul14_1700z['Wd [°]'], errors='coerce')

# Check the DataFrame
print(Jul14_1700z.dtypes)

Jul14_1700z = Jul14_1700z.sort_values(by = 'P [hPa]', ascending = True)

# Example from Sounderpy
# Create clean_data dictionary
clean_data = {}

# Add profile data | make sure to include p, z, T, Td, u, & v
clean_data['p']  = np.array(Jul14_1700z['P [hPa]']) * units.hPa
clean_data['z']  = np.array(Jul14_1700z['Geopot [m]']) * units.m
clean_data['T']  = np.array(Jul14_1700z['T [°C]']) * units.degC
clean_data['Td'] = np.array(Jul14_1700z['Dew [°C]']) * units.degC
clean_data['u']  = np.array(Jul14_1700z['Ws [m/s]']) * 1.94384 * units.kts
clean_data['v']  = np.array(Jul14_1700z['Ws [m/s]']) * 1.94384 * units.kts

# Declare site metadata
clean_data['site_info'] = {
            'site-id'     : 'Holland Airport',                  # Station ID, Site ID, etc.
            'site-name'   : 'Holland',                          # Location
            'site-lctn'   : 'MI',                               # State
            'site-latlon' : [42.7452, -86.0961],                # Lat/Lon
            'site-elv'    : 213,                                # Elevation (m)
            'source'      : 'MITTEN-CI',                        # Title of Project
            'model'       : 'none',                             # Model Name, if applicable
            'fcst-hour'   : f'none',                            # Forecast Hour, if applicable
            'run-time'    : ['none', 'none', 'none', 'none'],   # Model Run Time, if applicable
            'valid-time'  : ["2024", "07", "14", "17:00"]}      # Profile Valid Time

# Plot Titles
clean_data['titles'] = {
            'top_title': 'UNIVERSITY OF ILLINOIS: MITTEN-CI | OBSERVED SOUNDING',
            'left_title': 'VALID: 07/14/2024 - 17:00Z',
            'right_title': 'HOLLAND, MI [42.7452, -86.0961]    '}

spy.build_sounding(clean_data, dark_mode = False, special_parcels = 'none')

ValueError Traceback (most recent call last) Cell In[107], line 1 ----> 1 spy.build_sounding(clean_data, dark_mode = False, special_parcels = 'none')

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/sounderpy/sounderpy.py:237, in build_sounding(clean_data, style, color_blind, dark_mode, storm_motion, special_parcels, show_radar, radar_time, map_zoom, modify_sfc, save, filename) 235 __full_sounding(clean_data, color_blind, dark_mode, storm_motion, special_parcels, show_radar, radar_time, map_zoom, modify_sfc).savefig(filename, bbox_inches='tight') 236 else: --> 237 __full_sounding(clean_data, color_blind, dark_mode, storm_motion, special_parcels, show_radar, radar_time, map_zoom, modify_sfc).show()

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/sounderpy/plot.py:217, in __full_sounding(clean_data, color_blind, dark_mode, storm_motion, special_parcels, show_radar, radar_time, map_zoom, modify_sfc) 214 ws = mpcalc.wind_speed(u, v) 216 # calculate other sounding parameters using SounderPy Calc --> 217 general, thermo, kinem, intrp = sounding_params(sounding_data, storm_motion).calc() 218 ################################################################# 219 220 (...) 223 ### DECLARE PLOT TITLES FROM CLEAN_DATA ### 224 ################################################################# 226 top_title = clean_data['titles']['top_title']

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/sounderpy/calc.py:360, in sounding_params.calc(self) 356 thermo['mu6cape'] = mupcl.b6km 358 #--- MIXED LAYER PARCEL PROPERTIES ---# 359 # --------------------------------------------------------------- --> 360 mlpcl_p, mlpcl_T, mlpcl_Td = mpcalc.mixed_parcel(p, T, Td, bottom=p[0], depth=50*units.hPa, interpolate=True)[0:3] 361 mlpcl = parcelx(prof, flag=5, pres=mlpcl_p.m, tmpc=mlpcl_T.m, dwpc=mlpcl_Td.m) 362 thermo['mlT_trace'] = mlpcl.ttrace

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/metpy/xarray.py:1330, in preprocess_and_wrap..decorator..wrapper(*args, *kwargs) 1327 _mutate_arguments(boundargs, units.Quantity, lambda arg, : arg.m) 1329 # Evaluate inner calculation -> 1330 result = func(bound_args.args, **bound_args.kwargs) 1332 # Wrap output based on match and match_unit 1333 if match is None:

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/metpy/units.py:333, in check_units..dec..wrapper(*args, kwargs) 330 @functools.wraps(func) 331 def wrapper(*args, *kwargs): 332 _check_units_inner_helper(func, sig, defaults, dims, args, kwargs) --> 333 return func(*args, **kwargs)

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/metpy/calc/thermo.py:3351, in mixed_parcel(pressure, temperature, dewpoint, parcel_start_pressure, height, bottom, depth, interpolate) 3348 mixing_ratio = saturation_mixing_ratio(pressure, dewpoint) 3350 # Mix the variables over the layer -> 3351 mean_theta, mean_mixing_ratio = mixed_layer(pressure, theta, mixing_ratio, bottom=bottom, 3352 height=height, depth=depth, 3353 interpolate=interpolate) 3355 # Convert back to temperature 3356 mean_temperature = mean_theta * exner_function(parcel_start_pressure)

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/metpy/xarray.py:1330, in preprocess_and_wrap..decorator..wrapper(*args, *kwargs) 1327 _mutate_arguments(boundargs, units.Quantity, lambda arg, : arg.m) 1329 # Evaluate inner calculation -> 1330 result = func(bound_args.args, **bound_args.kwargs) 1332 # Wrap output based on match and match_unit 1333 if match is None:

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/metpy/units.py:333, in check_units..dec..wrapper(*args, kwargs) 330 @functools.wraps(func) 331 def wrapper(*args, *kwargs): 332 _check_units_inner_helper(func, sig, defaults, dims, args, kwargs) --> 333 return func(*args, **kwargs)

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/metpy/calc/thermo.py:3441, in mixed_layer(pressure, height, bottom, depth, interpolate, args) 3439 if depth is None: 3440 depth = units.Quantity(100, 'hPa') -> 3441 layer = get_layer(pressure, args, height=height, bottom=bottom, 3442 depth=depth, interpolate=interpolate) 3443 p_layer = layer[0] 3444 datavars_layer = layer[1:]

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/metpy/xarray.py:1330, in preprocess_and_wrap..decorator..wrapper(*args, *kwargs) 1327 _mutate_arguments(boundargs, units.Quantity, lambda arg, : arg.m) 1329 # Evaluate inner calculation -> 1330 result = func(bound_args.args, **bound_args.kwargs) 1332 # Wrap output based on match and match_unit 1333 if match is None:

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/metpy/units.py:333, in check_units..dec..wrapper(*args, kwargs) 330 @functools.wraps(func) 331 def wrapper(*args, *kwargs): 332 _check_units_inner_helper(func, sig, defaults, dims, args, kwargs) --> 333 return func(*args, **kwargs)

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/metpy/calc/tools.py:603, in get_layer(pressure, height, bottom, depth, interpolate, *args) 600 else: 601 raise ValueError('Depth must be specified in units of length or pressure') --> 603 toppressure, = _get_bound_pressure_height(pressure, top, height=height, 604 interpolate=interpolate) 606 ret = [] # returned data variables in layer 608 # Ensure pressures are sorted in ascending order

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/metpy/calc/tools.py:427, in _get_bound_pressure_height(pressure, bound, height, interpolate) 424 # If the bound is out of the range of the data, we shouldn't extrapolate 425 if not (_greater_or_close(bound_pressure, np.nanmin(pressure)) 426 and _less_or_close(bound_pressure, np.nanmax(pressure))): --> 427 raise ValueError('Specified bound is outside pressure range.') 428 if height is not None and not (_less_or_close(bound_height, np.nanmax(height)) 429 and _greater_or_close(bound_height, np.nanmin(height))): 430 raise ValueError('Specified bound is outside height range.')

ValueError: Specified bound is outside pressure range.

kylejgillett commented 4 days ago

Hi @wreed1989 ! Thanks for trying out SounderPy!

Cool to see some MITTEN-CI data being used! (CMU is my Alma mater😎).

It looking like the failure comes during the calculation of the mixed-level parcel (using MetPy).

What does your pressure array look like?

wreed1989 commented 3 days ago

The pressure array looks fine from what I can tell. I actually may have solved it by sorting in descending order instead of ascending. However, that has lead to the next error with 't_eval' not being properly sorted. Seems like if I sort one, it's flagging another haha 🤦🏻‍♂️

kylejgillett commented 3 days ago

Ok great! Just make sure all of your arrays are sorted and flipped correctly -- i.e., whatever you do to the pressure array must be done to the other arrays.

wreed1989 commented 3 days ago

I sent you a message on Twitter to troubleshoot a little further.

wreed1989 commented 3 days ago

Here is the t_eval error:


ValueError Traceback (most recent call last) Cell In[98], line 1 ----> 1 spy.build_sounding(clean_data, dark_mode = False, special_parcels = 'none')

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/sounderpy/sounderpy.py:237, in build_sounding(clean_data, style, color_blind, dark_mode, storm_motion, special_parcels, show_radar, radar_time, map_zoom, modify_sfc, save, filename) 235 __full_sounding(clean_data, color_blind, dark_mode, storm_motion, special_parcels, show_radar, radar_time, map_zoom, modify_sfc).savefig(filename, bbox_inches='tight') 236 else: --> 237 __full_sounding(clean_data, color_blind, dark_mode, storm_motion, special_parcels, show_radar, radar_time, map_zoom, modify_sfc).show()

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/sounderpy/plot.py:217, in __full_sounding(clean_data, color_blind, dark_mode, storm_motion, special_parcels, show_radar, radar_time, map_zoom, modify_sfc) 214 ws = mpcalc.wind_speed(u, v) 216 # calculate other sounding parameters using SounderPy Calc --> 217 general, thermo, kinem, intrp = sounding_params(sounding_data, storm_motion).calc() 218 ################################################################# 219 220 (...) 223 ### DECLARE PLOT TITLES FROM CLEAN_DATA ### 224 ################################################################# 226 top_title = clean_data['titles']['top_title']

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/sounderpy/calc.py:388, in sounding_params.calc(self) 384 #--- DOWNDRAFT CAPE ---# 385 # --------------------------------------------------------------- 387 if not np.any(np.isnan(sounding_data['Td'])): --> 388 thermo['dcape'], thermo['dcin'], thermo['dparcel_p'], thermo['dparcel_T'] = dcape_calc(sounding_data['p'], 389 sounding_data['T'], 390 sounding_data['Td']) 391 else: 392 thermo['dcape'] = ma.masked

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/sounderpy/calc.py:66, in dcape_calc(pressure, temperature, dewpoint) 64 # Descend parcel moist adiabatically to surface 65 down_pressure = pressure[pressure >= parcel_start_p].to(units.hPa) ---> 66 down_parcel_trace = mpcalc.moist_lapse(down_pressure, parcel_start_wb, 67 reference_pressure=parcel_start_p) 69 # Find virtual temperature of parcel and environment 70 parcel_virt_temp = mpcalc.virtual_temperature_from_dewpoint(down_pressure, down_parcel_trace, 71 down_parcel_trace)

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/metpy/xarray.py:1330, in preprocess_and_wrap..decorator..wrapper(*args, *kwargs) 1327 _mutate_arguments(boundargs, units.Quantity, lambda arg, : arg.m) 1329 # Evaluate inner calculation -> 1330 result = func(bound_args.args, **bound_args.kwargs) 1332 # Wrap output based on match and match_unit 1333 if match is None:

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/metpy/units.py:389, in process_units..dec..wrapper(*args, *kwargs) 386 _mutate_arguments(boundargs, units.Quantity, lambda val, : val.to_base_units().m) 388 # Evaluate inner calculation --> 389 result = func(bound_args.args, **bound_args.kwargs) 391 # Wrap output 392 if multiple_output:

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/metpy/calc/thermo.py:389, in moist_lapse(pressure, temperature, reference_pressure) 386 if np.any(points_below): 387 # Integrate downward 388 press_side = pressure[points_below] --> 389 result = si.solve_ivp(t_span=(reference_pressure, press_side[-1]), 390 t_eval=press_side, **solver_args) 391 if result.success: 392 ret = np.concatenate((ret, result.y), axis=-1)

File /opt/anaconda3/envs/devel/lib/python3.10/site-packages/scipy/integrate/_ivp/ivp.py:608, in solve_ivp(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, args, **options) 606 d = np.diff(t_eval) 607 if tf > t0 and np.any(d <= 0) or tf < t0 and np.any(d >= 0): --> 608 raise ValueError("Values in t_eval are not properly sorted.") 610 if tf > t0: 611 t_eval_i = 0

ValueError: Values in t_eval are not properly sorted.

kylejgillett commented 5 hours ago

For anyone who finds this issue in the future, we discovered that...

  1. the pressure array needed to be flipped (and thus so did the rest of the arrays)
  2. and there were repeat pressure values because the observations were so dense. This can cause disruptions with MetPy functions as they desire monotonically decreasing pressure values - always. Because of the number of observations, the pressure array (and all other arrays) could be sliced using "skips", i.e. pressure_array[0::5] (or use the zeroth value and every fifth value after that).