Repository: Unidata/MetPy

Issue:
interpolate_1d not properly working on ND arrays
While `interpolate_1d` advertises support for doing a 1D interpolation across a given axis for a multidimension array, this doesn't appear to be working correctly for `k_index`:

```python
p = units.Quantity([850, 700, 500], 'hPa')
t = np.random.randn(3, 151, 151) * units.degC
td = np.random.randn(3, 151, 151) * units.degC

mpcalc.k_index(p[:, None, None], t, td)
```

returns a single value. This is caused by `interpolate_1d` where:
```python
mpcalc.tools.interpolate_1d(units.Quantity([850, 700, 500],'hPa'), p[:, None, None], t, td)
```
also returns a single value.

Discussion:

Patch:
diff --git a/ci/requirements.txt b/ci/requirements.txt
index 0047e8bd09..94db527b93 100644
--- a/ci/requirements.txt
+++ b/ci/requirements.txt
@@ -1,5 +1,5 @@
 importlib_resources==5.10.0
-matplotlib==3.6.1
+matplotlib==3.6.0
 numpy==1.23.3
 pandas==1.5.0
 pooch==1.6.0
diff --git a/src/metpy/calc/thermo.py b/src/metpy/calc/thermo.py
index 7de2b3c59d..41bd29177b 100644
--- a/src/metpy/calc/thermo.py
+++ b/src/metpy/calc/thermo.py
@@ -2533,14 +2533,11 @@ def isentropic_interpolation(levels, pressure, temperature, *args, vertical_dim=
         fp = exner * (ka * t - a)
         return iter_log_p - (f / fp)
 
-    # Get dimensions in temperature
-    ndim = temperature.ndim
-
     # Convert units
     pressure = pressure.to('hPa')
     temperature = temperature.to('kelvin')
 
-    slices = [np.newaxis] * ndim
+    slices = [np.newaxis] * temperature.ndim
     slices[vertical_dim] = slice(None)
     slices = tuple(slices)
     pressure = units.Quantity(np.broadcast_to(pressure[slices].magnitude, temperature.shape),
@@ -2550,7 +2547,7 @@ def isentropic_interpolation(levels, pressure, temperature, *args, vertical_dim=
     sort_pressure = np.argsort(pressure.m, axis=vertical_dim)
     sort_pressure = np.swapaxes(np.swapaxes(sort_pressure, 0, vertical_dim)[::-1], 0,
                                 vertical_dim)
-    sorter = broadcast_indices(pressure, sort_pressure, ndim, vertical_dim)
+    sorter = broadcast_indices(sort_pressure, temperature.shape, vertical_dim)
     levs = pressure[sorter]
     tmpk = temperature[sorter]
 
diff --git a/src/metpy/calc/tools.py b/src/metpy/calc/tools.py
index f14a4199e5..262f22ef61 100644
--- a/src/metpy/calc/tools.py
+++ b/src/metpy/calc/tools.py
@@ -721,8 +721,8 @@ def find_bounding_indices(arr, values, axis, from_below=True):
         good[store_slice] = good_search
 
     # Create index values for broadcasting arrays
-    above = broadcast_indices(arr, indices, arr.ndim, axis)
-    below = broadcast_indices(arr, indices - 1, arr.ndim, axis)
+    above = broadcast_indices(indices, arr.shape, axis)
+    below = broadcast_indices(indices - 1, arr.shape, axis)
 
     return above, below, good
 
diff --git a/src/metpy/cbook.py b/src/metpy/cbook.py
index 33810d826c..c6801055b8 100644
--- a/src/metpy/cbook.py
+++ b/src/metpy/cbook.py
@@ -108,19 +108,23 @@ class Registry:
         return self._registry[name]
 
 
-def broadcast_indices(x, minv, ndim, axis):
+def broadcast_indices(indices, shape, axis):
     """Calculate index values to properly broadcast index array within data array.
 
-    See usage in interp.
+    The purpose of this function is work around the challenges trying to work with arrays of
+    indices that need to be "broadcast" against full slices for other dimensions.
+
+    See usage in `interpolate_1d` or `isentropic_interpolation`.
     """
     ret = []
+    ndim = len(shape)
     for dim in range(ndim):
         if dim == axis:
-            ret.append(minv)
+            ret.append(indices)
         else:
             broadcast_slice = [np.newaxis] * ndim
             broadcast_slice[dim] = slice(None)
-            dim_inds = np.arange(x.shape[dim])
+            dim_inds = np.arange(shape[dim])
             ret.append(dim_inds[tuple(broadcast_slice)])
     return tuple(ret)
 
diff --git a/src/metpy/interpolate/one_dimension.py b/src/metpy/interpolate/one_dimension.py
index e9a9f3ea63..09d4b9cf21 100644
--- a/src/metpy/interpolate/one_dimension.py
+++ b/src/metpy/interpolate/one_dimension.py
@@ -85,6 +85,7 @@ def interpolate_1d(x, xp, *args, axis=0, fill_value=np.nan, return_list_always=F
 
     Examples
     --------
+     >>> import metpy.interpolate
      >>> x = np.array([1., 2., 3., 4.])
      >>> y = np.array([1., 2., 3., 4.])
      >>> x_interp = np.array([2.5, 3.5])
@@ -102,24 +103,28 @@ def interpolate_1d(x, xp, *args, axis=0, fill_value=np.nan, return_list_always=F
     # Make x an array
     x = np.asanyarray(x).reshape(-1)
 
-    # Save number of dimensions in xp
-    ndim = xp.ndim
-
     # Sort input data
     sort_args = np.argsort(xp, axis=axis)
     sort_x = np.argsort(x)
 
+    # The shape after all arrays are broadcast to each other
+    # Can't use broadcast_shapes until numpy >=1.20 is our minimum
+    final_shape = np.broadcast(xp, *args).shape
+
     # indices for sorting
-    sorter = broadcast_indices(xp, sort_args, ndim, axis)
+    sorter = broadcast_indices(sort_args, final_shape, axis)
 
-    # sort xp
+    # sort xp -- need to make sure it's been manually broadcast due to our use of indices
+    # along all axes.
+    xp = np.broadcast_to(xp, final_shape)
     xp = xp[sorter]
-    # Ensure pressure in increasing order
+
+    # Ensure source arrays are also in sorted order
     variables = [arr[sorter] for arr in args]
 
     # Make x broadcast with xp
     x_array = x[sort_x]
-    expand = [np.newaxis] * ndim
+    expand = [np.newaxis] * len(final_shape)
     expand[axis] = slice(None)
     x_array = x_array[tuple(expand)]
 
@@ -140,8 +145,8 @@ def interpolate_1d(x, xp, *args, axis=0, fill_value=np.nan, return_list_always=F
         minv2[minv == 0] = 1
 
     # Get indices for broadcasting arrays
-    above = broadcast_indices(xp, minv2, ndim, axis)
-    below = broadcast_indices(xp, minv2 - 1, ndim, axis)
+    above = broadcast_indices(minv2, final_shape, axis)
+    below = broadcast_indices(minv2 - 1, final_shape, axis)
 
     if np.any(x_array < xp[below]):
         warnings.warn('Interpolation point out of data bounds encountered')


Test patch:
diff --git a/tests/interpolate/test_one_dimension.py b/tests/interpolate/test_one_dimension.py
index 2bb9e71e18..eeee8e4b45 100644
--- a/tests/interpolate/test_one_dimension.py
+++ b/tests/interpolate/test_one_dimension.py
@@ -182,3 +182,13 @@ def test_interpolate_masked_units():
     y_interp_truth = np.array([65., 75.]) * units.degC
     y_interp = interpolate_1d(x_interp, x, y)
     assert_array_almost_equal(y_interp, y_interp_truth, 7)
+
+
+def test_interpolate_broadcast():
+    """Test interpolate_1d with input levels needing broadcasting."""
+    p = units.Quantity([850, 700, 500], 'hPa')
+    t = units.Quantity(np.arange(60).reshape(3, 4, 5), 'degC')
+
+    t_level = interpolate_1d(units.Quantity(700, 'hPa'), p[:, None, None], t)
+    assert_array_almost_equal(t_level,
+                              units.Quantity(np.arange(20., 40.).reshape(1, 4, 5), 'degC'), 7)