Error in parcel_profile

[sgdecker]

My spidey sense says I'm overlooking something, but maybe this really is a bug?

Here's my code:

from datetime import datetime, timedelta

import xarray as xr
from xarray.backends import NetCDF4DataStore
from siphon.catalog import TDSCatalog

import metpy.calc as mpcalc


def get_nam212(init_time, valid_time):
    """Obtain NAM data on the 212 grid."""
    ymd = init_time.strftime('%Y%m%d')
    hr = init_time.strftime('%H')
    filename = f'{ymd}_{hr}00.grib2'
    ds_name = 'NAM_CONUS_40km_conduit_' + filename
    cat_name = ('https://thredds.ucar.edu/thredds/catalog/grib/NCEP/NAM/'
                'CONUS_40km/conduit/' + ds_name + '/catalog.xml')

    cat = TDSCatalog(cat_name)
    ds = cat.datasets[ds_name]
    ncss = ds.subset()
    query = ncss.query()
    query.time(valid_time).variables('all')
    nc = ncss.get_data(query)
    data = xr.open_dataset(NetCDF4DataStore(nc)).metpy.parse_cf()
    return data


init_time = datetime(2023, 2, 7, 12)
plot_time = init_time + timedelta(hours=6)

nam = get_nam212(init_time, plot_time)

tmp = nam['Temperature_isobaric']
dwp = mpcalc.dewpoint_from_relative_humidity(tmp, nam['Relative_humidity_isobaric'])

p = tmp['isobaric2']

p_env = p[::-1]
t_env = tmp[0, ::-1, 0, 0]
d_env = dwp[0, ::-1, 0, 0]

print(p_env.shape, t_env.shape, d_env.shape)

prof = mpcalc.parcel_profile(p_env, t_env, d_env)

The result is:

/home/decker/local/miniconda3/envs/meghan/lib/python3.11/site-packages/metpy/calc/thermo.py:1384: RuntimeWarning: divide by zero encountered in log
  val = np.log(vapor_pressure / mpconsts.nounit.sat_pressure_0c)
/home/decker/local/miniconda3/envs/meghan/lib/python3.11/site-packages/metpy/calc/thermo.py:1385: RuntimeWarning: invalid value encountered in divide
  return mpconsts.nounit.zero_degc + 243.5 * val / (17.67 - val)
(39,) (39,) (39,)
Traceback (most recent call last):
  File "/chariton/decker/test/metpy/xarraytest.py", line 45, in <module>
    prof = mpcalc.parcel_profile(p_env, t_env, d_env)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/decker/local/miniconda3/envs/meghan/lib/python3.11/site-packages/metpy/xarray.py", line 1329, in wrapper
    result = func(*bound_args.args, **bound_args.kwargs)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/decker/local/miniconda3/envs/meghan/lib/python3.11/site-packages/metpy/units.py", line 320, in wrapper
    return func(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^
  File "/home/decker/local/miniconda3/envs/meghan/lib/python3.11/site-packages/metpy/calc/thermo.py", line 976, in parcel_profile
    _, _, _, t_l, _, t_u = _parcel_profile_helper(pressure, temperature, dewpoint)
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/decker/local/miniconda3/envs/meghan/lib/python3.11/site-packages/metpy/calc/thermo.py", line 1179, in _parcel_profile_helper
    temp_lower = dry_lapse(press_lower, temperature)
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/decker/local/miniconda3/envs/meghan/lib/python3.11/site-packages/metpy/xarray.py", line 1329, in wrapper
    result = func(*bound_args.args, **bound_args.kwargs)
             ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/decker/local/miniconda3/envs/meghan/lib/python3.11/site-packages/metpy/units.py", line 320, in wrapper
    return func(*args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^
  File "/home/decker/local/miniconda3/envs/meghan/lib/python3.11/site-packages/metpy/calc/thermo.py", line 255, in dry_lapse
    return temperature * (pressure / reference_pressure)**mpconsts.kappa
           ~~~~~~~~~~~~^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  File "/home/decker/local/miniconda3/envs/meghan/lib/python3.11/site-packages/pint/facets/plain/quantity.py", line 1131, in __mul__
    return self._mul_div(other, operator.mul)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/decker/local/miniconda3/envs/meghan/lib/python3.11/site-packages/pint/facets/plain/quantity.py", line 102, in wrapped
    return f(self, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/decker/local/miniconda3/envs/meghan/lib/python3.11/site-packages/pint/facets/plain/quantity.py", line 82, in wrapped
    result = f(self, *args, **kwargs)
             ^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/decker/local/miniconda3/envs/meghan/lib/python3.11/site-packages/pint/facets/plain/quantity.py", line 1119, in _mul_div
    magnitude = magnitude_op(new_self._magnitude, other._magnitude)
                ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: operands could not be broadcast together with shapes (39,) (40,) 

I'm guessing the extra level has something to do with the LCL, and indeed there is no error using parcel_profile_with_lcl, but shouldn't this variation work also?

[dopplershift]

parcel_profile doesn't take arrays for temperature and dewpoint. It wants a sequence of pressures for which the profile should be calculated, and the starting thermodynamic state for the parcel (i.e. temperature and dewpoint at pressure[0]).

[sgdecker]

OK, I see I am violating the "contract" of the parcel_profile function. Expanding parcel_profile's abilities to do the right thing in this case (and be consistent with parcel_profile_with_lcl, whose docstring does not explicitly state that the parcel at index 0 is the one that is lifted, but that's what it does) would be nice, and maybe even easier than generating a better error message, so I think I'll cook up a feature request.

This would make the calling side consistent, i.e.,

prof = mpcalc.parcel_profile(p_env, t_env, d_env)
prof2 = mpcalc.parcel_profile_with_lcl(p_env, t_env, d_env)

vs.

prof = mpcalc.parcel_profile(p_env, t_env[0], d_env[0])
prof2 = mpcalc.parcel_profile_with_lcl(p_env, t_env, d_env)

(although the second form would still work also)

[dopplershift]

I would be happy to see the feature request that clearly outlines what your expected behavior would be. For a bit of history, parcel_profile(), requests the bare minimum it needs to do its job: the pressures you want and what the starting conditions of the parcel are (T0, Td0).

It later became clear that we needed another function that made it easy to get the LCL and include it in existing profiles--and allow users to avoid messy code manually interpolating the profiles of T/Td to the LCL. This of course necessitated passing those full profiles into parcel_profile_with_lcl. So while I see how the interfaces being different can create confusion, it's because we are expressly trying to avoid requiring full profiles in parcel_profile since we don't need them. This permits calculating a theoretical parcel profile based only on surface obs, for example.