Repository: Unidata/MetPy

Issue:
ValueError with mpcalc.parcel_profile
### Discussed in https://github.com/Unidata/MetPy/discussions/2308

<div type='discussions-op-text'>

<sup>Originally posted by **ejt4m8** January 25, 2022</sup>
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 <module>
----> 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.<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 ~\anaconda3\envs\metpy\lib\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 ~\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.<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 ~\anaconda3\envs\metpy\lib\site-packages\metpy\units.py:344, in process_units.<locals>.dec.<locals>.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.</div>

Discussion:
dopplershift: The rule for `parcel_profile` is to give back the temperatures at the pressures requested, so we need to do that. Otherwise, there's no way to communicate programmatically what's been dropped.

We *could* let `parcel_profile_with_lcl` just drop the duplicates, but I'm not sure it's worth doing special handling. This problem really goes away once we stop using Wyoming's parsed data (which only has 0.1 hPa resolution).
dcamron: That also matches pre-1.1 behavior, now we'll just warn people that it's happening. Changes made. I can let people know of this behavior in the docstring if needed.
dcamron: I've updated the docstrings for these two. Go ahead and take a look if you'd like. I'm going to throw together a quick test for duplicates in the `_as_dataset` calculation before we get this in, though.
dcamron: _Now_ it's done; apologies!

Patch:
diff --git a/docs/_templates/autosummary/module.rst b/docs/_templates/autosummary/module.rst
index d0bc51d266..86e72828eb 100644
--- a/docs/_templates/autosummary/module.rst
+++ b/docs/_templates/autosummary/module.rst
@@ -3,54 +3,42 @@
 {% include 'overrides/' ~ fullname ~ '.rst' with context %}
 
 {% else %}
-   {% if name == 'io' %}
-      {% set nice_name = 'Reading Data' %}
-   {% else %}
-      {% set nice_name = name | title | escape %}
-   {% endif %}
 
-{{ (nice_name ~ ' ``(' ~ fullname ~ ')``')|underline }}
+{{ objname }}
+{{ underline }}
 
 .. automodule:: {{ fullname }}
 
-   {% block attributes %}
-   {% if attributes %}
-   .. rubric:: {{ _('Module Attributes') }}
-
-   .. autosummary::
-   {% for item in attributes %}
-      {{ item }}
-   {%- endfor %}
-   {% endif %}
-   {% endblock %}
-
    {% block functions %}
    {% if functions %}
-   .. rubric:: {{ _('Functions') }}
+   .. rubric:: Functions
 
    .. autosummary::
       :toctree: ./
+
    {% for item in functions %}
-      {{ item }}
+      {{ item}}
    {%- endfor %}
    {% endif %}
    {% endblock %}
 
    {% block classes %}
    {% if classes %}
-   .. rubric:: {{ _('Classes') }}
+   .. rubric:: Classes
 
    .. autosummary::
       :toctree: ./
+
    {% for item in classes %}
       {{ item }}
    {%- endfor %}
+
    {% endif %}
    {% endblock %}
 
    {% block exceptions %}
    {% if exceptions %}
-   .. rubric:: {{ _('Exceptions') }}
+   .. rubric:: Exceptions
 
    .. autosummary::
       :toctree: ./
@@ -60,4 +48,4 @@
    {% endif %}
    {% endblock %}
 
-{% endif %}
+{% endif %}
\ No newline at end of file
diff --git a/docs/_templates/overrides/metpy.calc.rst b/docs/_templates/overrides/metpy.calc.rst
index f6eeeaf372..089ea63390 100644
--- a/docs/_templates/overrides/metpy.calc.rst
+++ b/docs/_templates/overrides/metpy.calc.rst
@@ -1,5 +1,5 @@
-Calculations ``(metpy.calc)``
-=============================
+calc
+==========
 
 .. automodule:: metpy.calc
 
diff --git a/docs/_templates/overrides/metpy.plots.ctables.rst b/docs/_templates/overrides/metpy.plots.ctables.rst
index e59d608373..64fe3b3b4f 100644
--- a/docs/_templates/overrides/metpy.plots.ctables.rst
+++ b/docs/_templates/overrides/metpy.plots.ctables.rst
@@ -1,5 +1,7 @@
-Colortables ``(metpy.plots.ctables)``
-=====================================
+
+
+ctables
+===================
 
 .. automodule:: metpy.plots.ctables
 
diff --git a/docs/_templates/overrides/metpy.xarray.rst b/docs/_templates/overrides/metpy.xarray.rst
index bfc6d154fb..5a02e12a7e 100644
--- a/docs/_templates/overrides/metpy.xarray.rst
+++ b/docs/_templates/overrides/metpy.xarray.rst
@@ -1,5 +1,5 @@
-Xarray Integration ``(metpy.xarray)``
-=====================================
+xarray
+==========
 
 .. automodule:: metpy.xarray
 
diff --git a/docs/conf.py b/docs/conf.py
index da0eecb48f..0960987843 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -90,7 +90,7 @@ source_suffix = ['.rst', '.md']
 
 # Controlling automatically generating summary tables in the docs
 autosummary_generate = True
-autosummary_ignore_module_all = False
+autosummary_imported_members = True
 
 # The encoding of source files.
 # source_encoding = 'utf-8-sig'
diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py
index 4e02366335..79f635c304 100644
--- a/src/metpy/calc/thermo.py
+++ b/src/metpy/calc/thermo.py
@@ -416,6 +416,10 @@ def lcl(pressure, temperature, dewpoint, max_iters=50, eps=1e-5):
 
     The function is guaranteed to finish by virtue of the `max_iters` counter.
 
+    Only functions on 1D profiles (not higher-dimension vertical cross sections or grids).
+    Since this function returns scalar values when given a profile, this will return Pint
+    Quantities even when given xarray DataArray profiles.
+
     .. versionchanged:: 1.0
        Renamed ``dewpt`` parameter to ``dewpoint``
 
@@ -744,11 +748,17 @@ def parcel_profile(pressure, temperature, dewpoint):
 
     See Also
     --------
-    lcl, moist_lapse, dry_lapse
+    lcl, moist_lapse, dry_lapse, parcel_profile_with_lcl, parcel_profile_with_lcl_as_dataset
 
     Notes
     -----
     Only functions on 1D profiles (not higher-dimension vertical cross sections or grids).
+    Duplicate pressure levels return duplicate parcel temperatures. Consider preprocessing
+    low-precision, high frequency profiles with tools like `scipy.medfilt`,
+    `pandas.drop_duplicates`, or `numpy.unique`.
+
+    Will only return Pint Quantities, even when given xarray DataArray profiles. To
+    obtain a xarray Dataset instead, use `parcel_profile_with_lcl_as_dataset` instead.
 
     .. versionchanged:: 1.0
        Renamed ``dewpt`` parameter to ``dewpoint``
@@ -805,7 +815,11 @@ def parcel_profile_with_lcl(pressure, temperature, dewpoint):
     Notes
     -----
     Only functions on 1D profiles (not higher-dimension vertical cross sections or grids).
-    Also, will only return Pint Quantities, even when given xarray DataArray profiles. To
+    Duplicate pressure levels return duplicate parcel temperatures. Consider preprocessing
+    low-precision, high frequency profiles with tools like `scipy.medfilt`,
+    `pandas.drop_duplicates`, or `numpy.unique`.
+
+    Will only return Pint Quantities, even when given xarray DataArray profiles. To
     obtain a xarray Dataset instead, use `parcel_profile_with_lcl_as_dataset` instead.
 
     .. versionchanged:: 1.0
@@ -930,9 +944,20 @@ def _parcel_profile_helper(pressure, temperature, dewpoint):
         return (press_lower[:-1], press_lcl, units.Quantity(np.array([]), press_lower.units),
                 temp_lower[:-1], temp_lcl, units.Quantity(np.array([]), temp_lower.units))
 
-    # Find moist pseudo-adiabatic profile starting at the LCL
+    # Establish profile above LCL
     press_upper = concatenate((press_lcl, pressure[pressure < press_lcl]))
-    temp_upper = moist_lapse(press_upper, temp_lower[-1]).to(temp_lower.units)
+
+    # Remove duplicate pressure values from remaining profile. Needed for solve_ivp in
+    # moist_lapse. unique will return remaining values sorted ascending.
+    unique, indices, counts = np.unique(press_upper.m, return_inverse=True, return_counts=True)
+    unique = units.Quantity(unique, press_upper.units)
+    if np.any(counts > 1):
+        warnings.warn(f'Duplicate pressure(s) {unique[counts > 1]:~P} provided. '
+                      'Output profile includes duplicate temperatures as a result.')
+
+    # Find moist pseudo-adiabatic profile starting at the LCL, reversing above sorting
+    temp_upper = moist_lapse(unique[::-1], temp_lower[-1]).to(temp_lower.units)
+    temp_upper = temp_upper[::-1][indices]
 
     # Return profile pieces
     return (press_lower[:-1], press_lcl, press_upper[1:],
diff --git a/src/metpy/io/__init__.py b/src/metpy/io/__init__.py
index c1f0aabbc3..4632c46c1e 100644
--- a/src/metpy/io/__init__.py
+++ b/src/metpy/io/__init__.py
@@ -1,15 +1,11 @@
 # Copyright (c) 2015,2016,2018,2021 MetPy Developers.
 # Distributed under the terms of the BSD 3-Clause License.
 # SPDX-License-Identifier: BSD-3-Clause
-"""Tools for reading various file formats.
+"""Classes for reading various file formats.
 
-Classes supporting formats are written to take both file names (for local files) or file-like
-objects; this allows reading files that are already in memory
-(using :class:`python:io.StringIO`) or remote files
-(using :func:`~python:urllib.request.urlopen`).
-
-`station_info` is an instance of `StationLookup` to find information about station locations
-(e.g. latitude, longitude, altitude) from various sources.
+These classes are written to take both file names (for local files) or file-like objects;
+this allows reading files that are already in memory (using :class:`python:io.StringIO`)
+or remote files (using :func:`~python:urllib.request.urlopen`).
 """
 
 from .gempak import *  # noqa: F403
diff --git a/src/metpy/io/station_data.py b/src/metpy/io/station_data.py
index e1cdb019bb..9417e93854 100644
--- a/src/metpy/io/station_data.py
+++ b/src/metpy/io/station_data.py
@@ -3,7 +3,6 @@
 # SPDX-License-Identifier: BSD-3-Clause
 """Pull out station metadata."""
 from collections import ChainMap, namedtuple
-from collections.abc import Mapping
 from functools import cached_property
 
 import numpy as np
@@ -121,18 +120,8 @@ def _read_airports_file(input_file=None):
                          }).to_dict()
 
 
-@exporter.export
-class StationLookup(Mapping):
-    """Look up station information from multiple sources.
-
-    This class follows the `Mapping` protocol with station ID as the key. This makes it
-    possible to e.g. iterate over all locations and get all of a certain criteria:
-
-    >>> import metpy.io
-    >>> conus_stations = [s for s in metpy.io.station_info if s.startswith('K')]
-    >>> conus_stations[:3]
-    ['KEET', 'K8A0', 'KALX']
-    """
+class StationLookup:
+    """Look up station information from multiple sources."""
 
     @cached_property
     def tables(self):
@@ -142,14 +131,6 @@ class StationLookup(Mapping):
                         dict(_read_station_text_file()),
                         dict(_read_airports_file()))
 
-    def __len__(self):
-        """Get the number of stations."""
-        return len(self.tables)
-
-    def __iter__(self):
-        """Allow iteration over the stations."""
-        return iter(self.tables)
-
     def __getitem__(self, stid):
         """Lookup station information from the ID."""
         try:
diff --git a/src/metpy/units.py b/src/metpy/units.py
index 3ad136420b..bbf5c60071 100644
--- a/src/metpy/units.py
+++ b/src/metpy/units.py
@@ -6,8 +6,6 @@ r"""Module to provide unit support.
 This makes use of the :mod:`pint` library and sets up the default settings
 for good temperature support.
 
-See Also: :doc:`Working with Units </tutorials/unit_tutorial>`.
-
 Attributes
 ----------
 units : :class:`pint.UnitRegistry`
@@ -214,12 +212,12 @@ def _check_argument_units(args, defaults, dimensionality):
         # See if the value passed in is appropriate
         try:
             if val.dimensionality != parsed:
-                yield arg, val, val.units, need
+                yield arg, val.units, need
         # No dimensionality
         except AttributeError:
             # If this argument is dimensionless, don't worry
             if parsed != '':
-                yield arg, val, 'none', need
+                yield arg, 'none', need
 
 
 def _get_changed_version(docstring):
@@ -258,20 +256,12 @@ def _check_units_inner_helper(func, sig, defaults, dims, *args, **kwargs):
     if bad:
         msg = f'`{func.__name__}` given arguments with incorrect units: '
         msg += ', '.join(
-            f'`{arg}` requires "{req}" but given "{given}"' for arg, _, given, req in bad
+            f'`{arg}` requires "{req}" but given "{given}"' for arg, given, req in bad
         )
         if 'none' in msg:
-            if any(isinstance(x, np.ma.core.MaskedArray) for _, x, _, _ in bad):
-                msg += ('\nA masked array `m` can be assigned a unit as follows:\n'
-                        '    from metpy.units import units\n'
-                        '    m = units.Quantity(m, "m/s")')
-            else:
-                msg += ('\nA xarray DataArray or numpy array `x` can be assigned a unit as '
-                        'follows:\n'
-                        '    from metpy.units import units\n'
-                        '    x = x * units("m/s")')
-            msg += ('\nFor more information see the Units Tutorial: '
-                    'https://unidata.github.io/MetPy/latest/tutorials/unit_tutorial.html')
+            msg += ('\nAny variable `x` can be assigned a unit as follows:\n'
+                    '    from metpy.units import units\n'
+                    '    x = units.Quantity(x, "m/s")')
 
         # If function has changed, mention that fact
         if func.__doc__:


Test patch:
diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py
index df53327214..cc8e5827bf 100644
--- a/tests/calc/test_thermo.py
+++ b/tests/calc/test_thermo.py
@@ -2092,3 +2092,68 @@ def test_cross_totals():
 
     ct = cross_totals(pressure, temperature, dewpoint)
     assert_almost_equal(ct, 21.40 * units.delta_degC, 2)
+
+
+def test_parcel_profile_drop_duplicates():
+    """Test handling repeat pressures in moist region of profile."""
+    pressure = np.array([962., 951., 937.9, 925., 908., 905.7, 894., 875.,
+                         41.3, 40.8, 37., 36.8, 32., 30., 27.7, 27.7, 26.4]) * units.hPa
+
+    temperature = units.Quantity(19.6, 'degC')
+
+    dewpoint = units.Quantity(18.6, 'degC')
+
+    truth = np.array([292.75, 291.78965331, 291.12778784, 290.61996294,
+                      289.93681828, 289.84313902, 289.36183185, 288.5626898,
+                      135.46280886, 134.99220142, 131.27369084, 131.07055878,
+                      125.93977169, 123.63877507, 120.85291224, 120.85291224,
+                      119.20448296]) * units.kelvin
+
+    with pytest.warns(UserWarning, match='Duplicate pressure'):
+        profile = parcel_profile(pressure, temperature, dewpoint)
+        assert_almost_equal(profile, truth, 5)
+
+
+def test_parcel_profile_with_lcl_as_dataset_duplicates():
+    """Test that parcel profile dataset creation works with duplicate pressures in profile."""
+    pressure = np.array(
+        [951., 951., 937.9, 925., 908., 30., 27.7, 27.7, 26.4, 25.1]
+    ) * units.hPa
+
+    temperature = np.array(
+        [20., 20., 19.5, 19., 18.6, -58.5, -58.1, -58.1, -57.2, -56.2]
+    ) * units.degC
+
+    dewpoint = np.array(
+        [19.4, 19.4, 19., 18.6, 18.3, -73.5, -75.1, -75.1, -77., -78.8]
+    ) * units.degC
+
+    truth = xr.Dataset(
+        {
+            'ambient_temperature': (
+                ('isobaric',),
+                np.insert(temperature.m, 2, 19.679237747615478) * units.degC
+            ),
+            'ambient_dew_point': (
+                ('isobaric',),
+                np.insert(dewpoint.m, 2, 19.143390198092384) * units.degC
+            ),
+            'parcel_temperature': (
+                ('isobaric',),
+                [
+                    293.15, 293.15, 292.40749167, 292.22841462, 291.73069653, 291.06139433,
+                    125.22698955, 122.40534065, 122.40534065, 120.73573642, 119.0063293
+                ] * units.kelvin
+            )
+        },
+        coords={
+            'isobaric': (
+                'isobaric',
+                np.insert(pressure.m, 2, 942.6)
+            )
+        }
+    )
+
+    profile = parcel_profile_with_lcl_as_dataset(pressure, temperature, dewpoint)
+
+    xr.testing.assert_allclose(profile, truth, atol=1e-5)
diff --git a/tests/io/test_station_data.py b/tests/io/test_station_data.py
index 7e0e9b0d70..1d7b25b4d4 100644
--- a/tests/io/test_station_data.py
+++ b/tests/io/test_station_data.py
@@ -7,7 +7,7 @@ from numpy.testing import assert_almost_equal
 import pandas as pd
 import pytest
 
-from metpy.io import add_station_lat_lon, station_info
+from metpy.io import add_station_lat_lon
 
 
 def test_add_lat_lon_station_data():
@@ -39,19 +39,3 @@ def test_add_lat_lon_station_data_not_found():
 
     with pytest.raises(KeyError):
         add_station_lat_lon(df)
-
-
-def test_station_lookup_get_station():
-    """Test that you can get a station by ID from the lookup."""
-    assert station_info['KOUN'].id == 'KOUN'
-
-
-def test_station_lookup_len():
-    """Test that you can get the length of the station data."""
-    assert len(station_info) == 13798
-
-
-def test_station_lookup_iter():
-    """Test iterating over the station data."""
-    for stid in station_info:
-        assert stid in station_info
diff --git a/tests/units/test_units.py b/tests/units/test_units.py
index 013bafb2bb..0f5ff78bdf 100644
--- a/tests/units/test_units.py
+++ b/tests/units/test_units.py
@@ -104,16 +104,6 @@ def test_bad(func, args, kwargs, bad_parts):
         assert 'unitless_const' not in message
 
 
-@pytest.mark.parametrize('func', test_funcs, ids=['some kwargs', 'all kwargs', 'all pos'])
-def test_bad_masked_array(func):
-    """Test getting a masked array-specific message when missing units."""
-    with pytest.raises(ValueError) as exc:
-        func(np.ma.array([30]), 1000 * units.mbar, 1.0 * units('kg/m^3'), 1, 5.)
-
-    message = str(exc.value)
-    assert 'units.Quantity' in message
-
-
 def test_pandas_units_simple():
     """Simple unit attachment to two columns."""
     df = pd.DataFrame(data=[[1, 4], [2, 5], [3, 6]], columns=['cola', 'colb'])