Repository: Unidata/MetPy

Issue:
Units issue with isentropic_interpolation_as_dataset and CONUS404 data
### What went wrong?

I am trying to use **metpy.calc.isentropic_interpolation_as_dataset** to interpolate CONUS404 data to isentropic levels. The CONUS404 data is produced using WRF v3.9.1.1. As such, I am using xarray and xWRF to read and process the data for use with MetPy. Unfortunately, I am receiving an error from this function stating:

`isentropic_interpolation given arguments with incorrect units: pressure requires "[pressure]" but given "dimensionless"`

Pressure is not an input for **isentropic_interpolation_as_dataset** and therefore I have no control over the pressure units referred to here. It looks like this function calls on the **isentropic_interpolation** function of metpy.calc, which does require pressure as an input, so perhaps however **isentropic_interpolation_as_dataset** is calculating pressure from the temperature values provided does not work with WRF/CONUS404 data. I have previously used **isentropic_interpolation_as_dataset** AND **isentropic_interpolation** without issue, but that was with ERA5 atmospheric reanalysis GRIB data, which is CF compliant and allows for the use of **metpy.parse_cf**. 

Any insight into how I can fix this issue and complete the interpolation would be greatly appreciated. Thank you!!

Note: You will notice that the "CLAT" coordinate is dropped when reading in variables. This is because isentropic_interpolation_as_dataset would not use variables with two latitude dimensions (i.e., XLAT and CLAT).

### Operating System

Linux

### Version

1.5.1

### Python Version

3.9.16

### Code to Reproduce

```python
import metpy
import metpy.calc as mpcalc
from metpy.units import units
import xarray as xr
import xwrf
import numpy as np

#------------------------------------------------------------------------------
#--- Open the NetCDF file -----------------------------------------------------
filename = "/glade/work/kshourd/29_June_2012/wrfout_d01_2012-06-29_12:00:00"
#--- xWRF postprocess & grid destaggering
ds = xr.open_dataset(filename).xwrf.postprocess()
ds = ds.xwrf.destagger()

lat = ds["XLAT"].drop_vars("CLAT")
lon = ds["XLONG"].drop_vars("CLAT")

#--- Extract variables
z = ds['geopotential_height'].drop_vars("CLAT")
z = z.metpy.quantify()
print(z)

ua = ds['wind_east'].drop_vars("CLAT")
ua = ua.metpy.quantify()
va = ds['wind_north'].drop_vars("CLAT")
va = va.metpy.quantify()

pressure = ds["air_pressure"].drop_vars("CLAT")
pressure = pressure.metpy.quantify()

theta = ds["air_potential_temperature"].drop_vars("CLAT")
theta = theta.metpy.quantify()
temp  = mpcalc.temperature_from_potential_temperature(pressure, theta)
temp  = temp.metpy.quantify()

# Compute the geostrophic wind
print("Geostrophic wind calculation...")
ugeo, vgeo = mpcalc.geostrophic_wind(z)
ugeo = ugeo.metpy.quantify()
vgeo = vgeo.metpy.quantify()

#--- Isentropic interpolation w/ MetPy
thlv     = np.arange(298,350,1)
isen_lev = thlv*units.kelvin

isentropic = mpcalc.isentropic_interpolation_as_dataset(isen_lev.squeeze(), temp.squeeze(), ugeo.squeeze(), vgeo.squeeze(), ua.squeeze(), va.squeeze())

print(list(isentropic.variables))
```


### Errors, Traceback, and Logs

```pytb
/glade/work/kshourd/testing123.py:53: UserWarning: More than one time coordinate present for variable .
  isentropic = mpcalc.isentropic_interpolation_as_dataset(isen_lev.squeeze(), temp.squeeze(), ugeo.squeeze(), vgeo.squeeze(), ua.squeeze(), va.squeeze())
Traceback (most recent call last):
  File "/glade/work/kshourd/testing123.py", line 53, in <module>
    isentropic = mpcalc.isentropic_interpolation_as_dataset(isen_lev.squeeze(), temp.squeeze(), ugeo.squeeze(), vgeo.squeeze(), ua.squeeze(), va.squeeze())
  File "/glade/work/kshourd/conda-envs/wrfplot/lib/python3.9/site-packages/metpy/calc/thermo.py", line 2765, in isentropic_interpolation_as_dataset
    ret = isentropic_interpolation(
  File "/glade/work/kshourd/conda-envs/wrfplot/lib/python3.9/site-packages/metpy/xarray.py", line 1544, in wrapper
    return func(*bound_args.args, **bound_args.kwargs)
  File "/glade/work/kshourd/conda-envs/wrfplot/lib/python3.9/site-packages/metpy/xarray.py", line 1328, in wrapper
    result = func(*bound_args.args, **bound_args.kwargs)
  File "/glade/work/kshourd/conda-envs/wrfplot/lib/python3.9/site-packages/metpy/units.py", line 324, in wrapper
    _check_units_inner_helper(func, sig, defaults, dims, *args, **kwargs)
  File "/glade/work/kshourd/conda-envs/wrfplot/lib/python3.9/site-packages/metpy/units.py", line 311, in _check_units_inner_helper
    raise ValueError(msg)
ValueError: This function changed in 1.0--double check that the function is being called properly.
`isentropic_interpolation` given arguments with incorrect units: `pressure` requires "[pressure]" but given "dimensionless"
```


Discussion:

Patch:
diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py
index 0f36e33efd..79571066f5 100644
--- a/src/metpy/calc/thermo.py
+++ b/src/metpy/calc/thermo.py
@@ -2624,7 +2624,7 @@ def isentropic_interpolation(levels, pressure, temperature, *args, vertical_dim=
         One-dimensional array of desired potential temperature surfaces
 
     pressure : array-like
-        One-dimensional array of pressure levels
+        Array of pressure
 
     temperature : array-like
         Array of temperature
@@ -2691,10 +2691,17 @@ def isentropic_interpolation(levels, pressure, temperature, *args, vertical_dim=
     pressure = pressure.to('hPa')
     temperature = temperature.to('kelvin')
 
+    # Construct slices for broadcasting with temperature (used for pressure & theta levels)
     slices = [np.newaxis] * temperature.ndim
     slices[vertical_dim] = slice(None)
     slices = tuple(slices)
-    pressure = units.Quantity(np.broadcast_to(pressure[slices].magnitude, temperature.shape),
+
+    # For 1-D pressure, we assume it's the vertical coordinate and know how it should broadcast
+    # to the same shape as temperature. Otherwise, just assume it's ready for broadcast, or
+    # it has the same shape and is a no-op.
+    if pressure.ndim == 1:
+        pressure = pressure[slices]
+    pressure = units.Quantity(np.broadcast_to(pressure.magnitude, temperature.shape),
                               pressure.units)
 
     # Sort input data
@@ -2779,7 +2786,8 @@ def isentropic_interpolation_as_dataset(
     *args,
     max_iters=50,
     eps=1e-6,
-    bottom_up_search=True
+    bottom_up_search=True,
+    pressure=None
 ):
     r"""Interpolate xarray data in isobaric coords to isentropic coords, returning a Dataset.
 
@@ -2799,6 +2807,9 @@ def isentropic_interpolation_as_dataset(
     bottom_up_search : bool, optional
         Controls whether to search for levels bottom-up (starting at lower indices),
         or top-down (starting at higher indices). Defaults to True, which is bottom-up search.
+    pressure : `xarray.DataArray`, optional
+        Array of pressure to use when the vertical coordinate for the passed in data is not
+        pressure (e.g. data using sigma coordinates)
 
     Returns
     -------
@@ -2833,10 +2844,20 @@ def isentropic_interpolation_as_dataset(
                        'The output Dataset includes duplicate levels as a result. '
                        'This may cause xarray to crash when working with this Dataset!')
 
+    if pressure is None:
+        pressure = all_args[0].metpy.vertical
+
+    if (units.get_dimensionality(pressure.metpy.units)
+            != units.get_dimensionality('[pressure]')):
+        raise ValueError('The pressure array/vertical coordinate for the passed in data does '
+                         'not appear to be pressure. In this case you need to pass in '
+                         'pressure values with proper units using the `pressure` keyword '
+                         'argument.')
+
     # Obtain result as list of Quantities
     ret = isentropic_interpolation(
         levels,
-        all_args[0].metpy.vertical,
+        pressure,
         all_args[0].metpy.unit_array,
         *(arg.metpy.unit_array for arg in all_args[1:]),
         vertical_dim=all_args[0].metpy.find_axis_number('vertical'),


Test patch:
diff --git a/tests/calc/test_thermo.py b/tests/calc/test_thermo.py
index 5e7aea34c9..72347a1fd8 100644
--- a/tests/calc/test_thermo.py
+++ b/tests/calc/test_thermo.py
@@ -1293,9 +1293,12 @@ def test_isentropic_pressure_p_increase_rh_out():
     assert_almost_equal(isentprs[1], truerh, 3)
 
 
-def test_isentropic_pressure_interp():
+@pytest.mark.parametrize('press_3d', [False, True])
+def test_isentropic_pressure_interp(press_3d):
     """Test calculation of isentropic pressure function."""
     lev = [100000., 95000., 90000., 85000.] * units.Pa
+    if press_3d:
+        lev = np.lib.stride_tricks.broadcast_to(lev[:, None, None], (4, 5, 5))
     tmp = np.ones((4, 5, 5))
     tmp[0, :] = 296.
     tmp[1, :] = 292.
@@ -1464,6 +1467,78 @@ def test_isentropic_interpolation_as_dataset_duplicate(xarray_isentropic_data):
                                                 xarray_isentropic_data.rh)
 
 
+@pytest.fixture
+def xarray_sigma_isentropic_data():
+    """Generate test xarray dataset on sigma vertical coords for interpolation functions."""
+    return xr.Dataset(
+        {
+            'temperature': (
+                ('z', 'y', 'x'),
+                [[[296.]], [[292.]], [[290.]], [[288.]]] * units.K
+            ),
+            'rh': (
+                ('z', 'y', 'x'),
+                [[[100.]], [[80.]], [[40.]], [[20.]]] * units.percent
+            ),
+            'pressure': (
+                ('z', 'y', 'x'),
+                [[[1000.]], [[950.]], [[900.]], [[850.]]] * units.hPa
+            )
+        },
+        coords={
+            'z': (('z',), [0.98, 0.928, 0.876, 0.825], {'units': 'dimensionless'}),
+            'time': '2020-01-01T00:00Z'
+        }
+    )
+
+
+def test_isen_interpolation_as_dataset_non_pressure_default(xarray_sigma_isentropic_data):
+    """Test isentropic interpolation with xarray data with non-pressure vertical coord."""
+    isentlev = [296., 297.] * units.kelvin
+    with pytest.raises(ValueError, match='vertical coordinate for the.*does not'):
+        isentropic_interpolation_as_dataset(isentlev,
+                                            xarray_sigma_isentropic_data.temperature,
+                                            xarray_sigma_isentropic_data.rh)
+
+
+def test_isen_interpolation_as_dataset_passing_pressre(xarray_sigma_isentropic_data):
+    """Test isentropic interpolation with xarray when passing a pressure array."""
+    isentlev = [296., 297.] * units.kelvin
+    result = isentropic_interpolation_as_dataset(
+        isentlev, xarray_sigma_isentropic_data.temperature,
+        xarray_sigma_isentropic_data.rh, pressure=xarray_sigma_isentropic_data.pressure)
+    expected = xr.Dataset(
+        {
+            'pressure': (
+                ('isentropic_level', 'y', 'x'),
+                [[[1000.]], [[936.213]]] * units.hPa,
+                {'standard_name': 'air_pressure'}
+            ),
+            'temperature': (
+                ('isentropic_level', 'y', 'x'),
+                [[[296.]], [[291.4579]]] * units.K,
+                {'standard_name': 'air_temperature'}
+            ),
+            'rh': (
+                ('isentropic_level', 'y', 'x'),
+                [[[100.]], [[69.19706]]] * units.percent
+            )
+        },
+        coords={
+            'isentropic_level': (
+                ('isentropic_level',),
+                [296., 297.],
+                {'units': 'kelvin', 'positive': 'up'}
+            ),
+            'time': '2020-01-01T00:00Z'
+        }
+    )
+    xr.testing.assert_allclose(result, expected)
+    assert result['pressure'].attrs == expected['pressure'].attrs
+    assert result['temperature'].attrs == expected['temperature'].attrs
+    assert result['isentropic_level'].attrs == expected['isentropic_level'].attrs
+
+
 @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."""