Repository: Unidata/MetPy

Issue:
isentropic_interpolation_as_dataset returns unplottable dataset
### What went wrong?

Students are great at inadvertently stress-testing MetPy. Here is our latest discovery.

The error took a while to figure out (especially in the original code where the offending line is far away from the declarative plotting section). Although using the intended code:
```python
isen_lev = (295, 300, 305) * units('K')
```
avoids the error, I think it is a bug for `isentropic_interpolation_as_dataset` to return a Dataset that the declarative interface can't plot. I can see three ways to fix the bug:
1) Have `isentropic_interpolation_as_dataset` ignore duplicate levels and issue a warning to the user. In this example, that function would return a dataset with two levels rather than three.
2) Have `isentropic_interpolation_as_dataset` throw an error stating that requesting the same level twice is invalid.
3) If there is a legitimate use case for duplicate isentropic levels, the declarative interface should be able to generate plots from such datasets (especially when the requested level is not a duplicate).

### Operating System

Linux

### Version

1.5.1

### Python Version

3.11.5

### Code to Reproduce

```python
from datetime import datetime
import metpy.calc as mpcalc
from metpy.plots import ContourPlot, MapPanel, PanelContainer
from metpy.units import units
import xarray as xr
from xarray.backends import NetCDF4DataStore
from siphon.catalog import TDSCatalog


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


plot_time = datetime.utcnow().replace(hour=0, minute=0, second=0)
ds = get_nam212(plot_time, plot_time).metpy.assign_latitude_longitude()

hgt = ds['Geopotential_height_isobaric'].squeeze()
tmp = ds['Temperature_isobaric'].squeeze()

isen_lev = (295, 300, 300) * units('K')
ds_isen = mpcalc.isentropic_interpolation_as_dataset(isen_lev, tmp, hgt)

isobars = ContourPlot()
isobars.data = ds_isen
isobars.field = 'pressure'
isobars.level = 295 * units('K')
isobars.contours = range(50, 1050, 50)

mp = MapPanel()
mp.area = [-90, -60, 33, 49]
mp.layers = ['coastline', 'borders', 'states']
mp.plots = [isobars]

pc = PanelContainer()
pc.size = (10, 8)
pc.panels = [mp]
pc.show()
```


### Errors, Traceback, and Logs

```pytb
InvalidIndexError                         Traceback (most recent call last)
File /chariton/decker/test/bug/duplicate_level.py:52
     50 pc.size = (10, 8)
     51 pc.panels = [mp]
---> 52 pc.show()

File ~/local/miniconda3/envs/met433/lib/python3.11/site-packages/metpy/plots/declarative.py:179, in PanelContainer.show(self)
    177 def show(self):
    178     """Show the constructed graphic on the screen."""
--> 179     self.draw()
    180     plt.show()

File ~/local/miniconda3/envs/met433/lib/python3.11/site-packages/metpy/plots/declarative.py:166, in PanelContainer.draw(self)
    164 for panel in self.panels:
    165     with panel.hold_trait_notifications():
--> 166         panel.draw()

File ~/local/miniconda3/envs/met433/lib/python3.11/site-packages/metpy/plots/declarative.py:456, in MapPanel.draw(self)
    454     self.ax.set_global()
    455 elif self.area is not None:
--> 456     self.ax.set_extent(self.area, ccrs.PlateCarree())
    458 # Draw all of the plots.
    459 for p in self.plots:

File ~/local/miniconda3/envs/met433/lib/python3.11/site-packages/metpy/plots/declarative.py:428, in MapPanel.ax(self)
    425 # If we haven't actually made an instance yet, make one with the right size and
    426 # map projection.
    427 if getattr(self, '_ax', None) is None:
--> 428     self._ax = self.parent.figure.add_subplot(*self.layout, projection=self._proj_obj)
    430 return self._ax

File ~/local/miniconda3/envs/met433/lib/python3.11/site-packages/metpy/plots/declarative.py:364, in MapPanel._proj_obj(self)
    362 if isinstance(self.projection, str):
    363     if self.projection == 'data':
--> 364         if isinstance(self.plots[0].griddata, tuple):
    365             proj = self.plots[0].griddata[0].metpy.cartopy_crs
    366         else:

File ~/local/miniconda3/envs/met433/lib/python3.11/site-packages/metpy/plots/declarative.py:767, in PlotScalar.griddata(self)
    765     if selector is not None:
    766         subset[dim_coord] = selector
--> 767 data_subset = data.metpy.sel(**subset).squeeze()
    768 if (data_subset.ndim != 2):
    769     if data_subset.ndim == 3:

File ~/local/miniconda3/envs/met433/lib/python3.11/site-packages/metpy/xarray.py:644, in MetPyDataArrayAccessor.sel(self, indexers, method, tolerance, drop, **indexers_kwargs)
    642 indexers = either_dict_or_kwargs(indexers, indexers_kwargs, 'sel')
    643 indexers = _reassign_quantity_indexer(self._data_array, indexers)
--> 644 return self._data_array.sel(indexers, method=method, tolerance=tolerance, drop=drop)

File ~/local/miniconda3/envs/met433/lib/python3.11/site-packages/xarray/core/dataarray.py:1582, in DataArray.sel(self, indexers, method, tolerance, drop, **indexers_kwargs)
   1472 def sel(
   1473     self: T_DataArray,
   1474     indexers: Mapping[Any, Any] | None = None,
   (...)
   1478     **indexers_kwargs: Any,
   1479 ) -> T_DataArray:
   1480     """Return a new DataArray whose data is given by selecting index
   1481     labels along the specified dimension(s).
   1482 
   (...)
   1580     Dimensions without coordinates: points
   1581     """
-> 1582     ds = self._to_temp_dataset().sel(
   1583         indexers=indexers,
   1584         drop=drop,
   1585         method=method,
   1586         tolerance=tolerance,
   1587         **indexers_kwargs,
   1588     )
   1589     return self._from_temp_dataset(ds)

File ~/local/miniconda3/envs/met433/lib/python3.11/site-packages/xarray/core/dataset.py:3020, in Dataset.sel(self, indexers, method, tolerance, drop, **indexers_kwargs)
   2959 """Returns a new dataset with each array indexed by tick labels
   2960 along the specified dimension(s).
   2961 
   (...)
   3017 DataArray.sel
   3018 """
   3019 indexers = either_dict_or_kwargs(indexers, indexers_kwargs, "sel")
-> 3020 query_results = map_index_queries(
   3021     self, indexers=indexers, method=method, tolerance=tolerance
   3022 )
   3024 if drop:
   3025     no_scalar_variables = {}

File ~/local/miniconda3/envs/met433/lib/python3.11/site-packages/xarray/core/indexing.py:190, in map_index_queries(obj, indexers, method, tolerance, **indexers_kwargs)
    188         results.append(IndexSelResult(labels))
    189     else:
--> 190         results.append(index.sel(labels, **options))
    192 merged = merge_sel_results(results)
    194 # drop dimension coordinates found in dimension indexers
    195 # (also drop multi-index if any)
    196 # (.sel() already ensures alignment)

File ~/local/miniconda3/envs/met433/lib/python3.11/site-packages/xarray/core/indexes.py:760, in PandasIndex.sel(self, labels, method, tolerance)
    758 else:
    759     if method is not None:
--> 760         indexer = get_indexer_nd(
    761             self.index, label_array, method, tolerance
    762         )
    763         if np.any(indexer < 0):
    764             raise KeyError(
    765                 f"not all values found in index {coord_name!r}"
    766             )

File ~/local/miniconda3/envs/met433/lib/python3.11/site-packages/xarray/core/indexes.py:565, in get_indexer_nd(index, labels, method, tolerance)
    563 if flat_labels.dtype == "float16":
    564     flat_labels = flat_labels.astype("float64")
--> 565 flat_indexer = index.get_indexer(flat_labels, method=method, tolerance=tolerance)
    566 indexer = flat_indexer.reshape(labels.shape)
    567 return indexer

File ~/local/miniconda3/envs/met433/lib/python3.11/site-packages/pandas/core/indexes/base.py:3874, in Index.get_indexer(self, target, method, limit, tolerance)
   3871 self._check_indexing_method(method, limit, tolerance)
   3873 if not self._index_as_unique:
-> 3874     raise InvalidIndexError(self._requires_unique_msg)
   3876 if len(target) == 0:
   3877     return np.array([], dtype=np.intp)

InvalidIndexError: Reindexing only valid with uniquely valued Index objects
```


Discussion:

Patch:
diff --git a/.github/workflows/code-analysis.yml b/.github/workflows/code-analysis.yml
index 301d0c0ffd..43d963a1f0 100644
--- a/.github/workflows/code-analysis.yml
+++ b/.github/workflows/code-analysis.yml
@@ -31,9 +31,9 @@ jobs:
       uses: actions/checkout@v4
 
     - name: Initialize CodeQL
-      uses: github/codeql-action/init@v3
+      uses: github/codeql-action/init@v2
       with:
         config-file: ./.github/codeql/codeql-config.yml
 
     - name: Perform CodeQL Analysis
-      uses: github/codeql-action/analyze@v3
+      uses: github/codeql-action/analyze@v2
diff --git a/ci/linting_requirements.txt b/ci/linting_requirements.txt
index 173ccafd40..1202dc8a08 100644
--- a/ci/linting_requirements.txt
+++ b/ci/linting_requirements.txt
@@ -1,4 +1,4 @@
-ruff==0.1.8
+ruff==0.1.7
 
 flake8==6.1.0
 pycodestyle==2.11.1
@@ -7,7 +7,7 @@ pyflakes==3.1.0
 flake8-continuation==1.0.5
 flake8-copyright==0.2.4
 flake8-isort==6.1.1
-isort==5.13.2
+isort==5.13.1
 flake8-requirements==2.0.1
 
 flake8-rst-docstrings==0.3.0
diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py
index eeba567cd6..0f36e33efd 100644
--- a/src/metpy/calc/thermo.py
+++ b/src/metpy/calc/thermo.py
@@ -2826,6 +2826,13 @@ def isentropic_interpolation_as_dataset(
     # Ensure matching coordinates by broadcasting
     all_args = xr.broadcast(temperature, *args)
 
+    # Check for duplicate isentropic levels, which can be problematic (Issue #3309)
+    unique, counts = np.unique(levels.m, return_counts=True)
+    if np.any(counts > 1):
+        _warnings.warn(f'Duplicate level(s) {unique[counts > 1]} provided. '
+                       'The output Dataset includes duplicate levels as a result. '
+                       'This may cause xarray to crash when working with this Dataset!')
+
     # Obtain result as list of Quantities
     ret = isentropic_interpolation(
         levels,


Test patch:
diff --git a/ci/test_requirements.txt b/ci/test_requirements.txt
index 3fa0a52627..11b3ab94ac 100644
--- a/ci/test_requirements.txt
+++ b/ci/test_requirements.txt
@@ -2,4 +2,4 @@ packaging==23.2
 pytest==7.4.3
 pytest-mpl==0.16.1
 netCDF4==1.6.5
-coverage==7.3.3
+coverage==7.3.2
diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py
index 118902b3f2..5e7aea34c9 100644
--- a/tests/calc/test_thermo.py
+++ b/tests/calc/test_thermo.py
@@ -1386,29 +1386,10 @@ def test_isentropic_pressure_4d():
     assert_almost_equal(isentprs[1][:, 1, ], truerh, 3)
 
 
-def test_isentropic_interpolation_dataarray():
-    """Test calculation of isentropic interpolation with xarray dataarrays."""
-    temp = xr.DataArray([[[296.]], [[292.]], [[290.]], [[288.]]] * units.K,
-                        dims=('isobaric', 'y', 'x'),
-                        coords={'isobaric': (('isobaric',), [1000., 950., 900., 850.],
-                                             {'units': 'hPa'}),
-                                'time': '2020-01-01T00:00Z'})
-
-    rh = xr.DataArray([[[100.]], [[80.]], [[40.]], [[20.]]] * units.percent,
-                      dims=('isobaric', 'y', 'x'), coords={
-        'isobaric': (('isobaric',), [1000., 950., 900., 850.], {'units': 'hPa'}),
-        'time': '2020-01-01T00:00Z'})
-
-    isentlev = [296., 297.] * units.kelvin
-    press, rh_interp = isentropic_interpolation(isentlev, temp.isobaric, temp, rh)
-
-    assert_array_almost_equal(press, np.array([[[1000.]], [[936.213]]]) * units.hPa, 3)
-    assert_array_almost_equal(rh_interp, np.array([[[100.]], [[69.19706]]]) * units.percent, 3)
-
-
-def test_isentropic_interpolation_as_dataset():
-    """Test calculation of isentropic interpolation with xarray."""
-    data = xr.Dataset(
+@pytest.fixture
+def xarray_isentropic_data():
+    """Generate test xarray dataset for interpolation functions."""
+    return xr.Dataset(
         {
             'temperature': (
                 ('isobaric', 'y', 'x'),
@@ -1424,8 +1405,25 @@ def test_isentropic_interpolation_as_dataset():
             'time': '2020-01-01T00:00Z'
         }
     )
+
+
+def test_isentropic_interpolation_dataarray(xarray_isentropic_data):
+    """Test calculation of isentropic interpolation with xarray dataarrays."""
+    rh = xarray_isentropic_data.rh
+    temp = xarray_isentropic_data.temperature
+
+    isentlev = [296., 297.] * units.kelvin
+    press, rh_interp = isentropic_interpolation(isentlev, temp.isobaric, temp, rh)
+
+    assert_array_almost_equal(press, np.array([[[1000.]], [[936.213]]]) * units.hPa, 3)
+    assert_array_almost_equal(rh_interp, np.array([[[100.]], [[69.19706]]]) * units.percent, 3)
+
+
+def test_isentropic_interpolation_as_dataset(xarray_isentropic_data):
+    """Test calculation of isentropic interpolation with xarray."""
     isentlev = [296., 297.] * units.kelvin
-    result = isentropic_interpolation_as_dataset(isentlev, data['temperature'], data['rh'])
+    result = isentropic_interpolation_as_dataset(isentlev, xarray_isentropic_data.temperature,
+                                                 xarray_isentropic_data.rh)
     expected = xr.Dataset(
         {
             'pressure': (
@@ -1458,6 +1456,14 @@ def test_isentropic_interpolation_as_dataset():
     assert result['isentropic_level'].attrs == expected['isentropic_level'].attrs
 
 
+def test_isentropic_interpolation_as_dataset_duplicate(xarray_isentropic_data):
+    """Test duplicate-level check in isentropic_interpolation_as_dataset."""
+    isentlev = [296., 296.] * units.kelvin
+    with pytest.warns(UserWarning):
+        _ = isentropic_interpolation_as_dataset(isentlev, xarray_isentropic_data.temperature,
+                                                xarray_isentropic_data.rh)
+
+
 @pytest.mark.parametrize('array_class', (units.Quantity, masked_array))
 def test_surface_based_cape_cin(array_class):
     """Test the surface-based CAPE and CIN calculation."""