Repository: Unidata/MetPy

Issue:
interpolate.log_interpolate_1d fails with Xarray input of more than 3D (Numpy succeeds)
- Description: `metpy.interpolate.interpolate_1d` has `no implementation found for 'numpy.apply_along_axis'` error for cases where the input is `xarray.Dataarray` and has more than 3 dimensions (e.g. 4th dim of `time`). Function works fine when the same input  is `numpy.ndarray` as shown in the below code.

- Environment: Below sample is easy to create an environment for. The only soecific package other than common packages (xarray, nump, metpy) is `geocat-datafiles` that can be installed via:

`conda install -c ncar geocat-datafiles`

- Platform: MacOS 10.15.7
- Versions: Python 3.9.2, MetPy 1.0

```python
import xarray as xr
import numpy as np
import metpy.interpolate
import geocat.datafiles as gdf

# Open a netCDF data file using xarray default engine and load the data into xarrays
ds = xr.open_dataset(gdf.get("netcdf_files/atmos.nc"), decode_times=False)

# Extract the data needed
u = ds.U  # U component of wind
hyam = ds.hyam  # hybrid A coefficient
hybm = ds.hybm  # hybrid B coefficient
ps = ds.PS  # surface pressures in Pascals

# Define some constants
interp_axis = 1  # 'lev', interpolation axis
p0 = 100000.  # surface reference pressure in Pascals

# Specify output pressure levels
new_levels = np.array([1000, 950, 800, 700, 600, 500, 400, 300,
                       200])  # in millibars
new_levels = new_levels * 100  # convert to Pascals

# Calculate pressure levels at the hybrid levels
pressure = hyam * p0 + hybm * ps
# Make pressure shape same as data shape
pressure = pressure.transpose(*u.dims)

# Calling with Numpy is successful
output_np = metpy.interpolate.log_interpolate_1d(new_levels, pressure.values, u.values, axis=1)

# Calling with Xarray is failing!!!
output_xr = metpy.interpolate.log_interpolate_1d(new_levels, pressure, u, axis=1)

```




Discussion:

Patch:
diff --git a/src/metpy/interpolate/one_dimension.py b/src/metpy/interpolate/one_dimension.py
index bead0796de..f12bcea96f 100644
--- a/src/metpy/interpolate/one_dimension.py
+++ b/src/metpy/interpolate/one_dimension.py
@@ -239,4 +239,6 @@ def _strip_matching_units(*args):
     if all(hasattr(arr, 'units') for arr in args):
         return [arr.to(args[0].units).magnitude for arr in args]
     else:
-        return args
+        # Handle the case where we get mixed 'dimensionless' and bare array. This happens e.g.
+        # when you pass in a DataArray with no units for one arg.
+        return [arr.m_as('dimensionless') if hasattr(arr, 'units') else arr for arr in args]


Test patch:
diff --git a/tests/interpolate/test_one_dimension.py b/tests/interpolate/test_one_dimension.py
index eeee8e4b45..4c473bb4c5 100644
--- a/tests/interpolate/test_one_dimension.py
+++ b/tests/interpolate/test_one_dimension.py
@@ -5,6 +5,7 @@
 
 import numpy as np
 import pytest
+import xarray as xr
 
 from metpy.interpolate import interpolate_1d, interpolate_nans_1d, log_interpolate_1d
 from metpy.testing import assert_array_almost_equal
@@ -49,6 +50,16 @@ def test_log_interpolate_1d():
     assert_array_almost_equal(y_interp, y_interp_truth, 7)
 
 
+def test_log_interpolate_1d_mixed():
+    """Test log interpolation with a mix of compatible input types."""
+    x_log = xr.DataArray([1e3, 1e4, 1e5, 1e6])
+    y_log = np.log(x_log) * 2 + 3
+    x_interp = np.array([5e3, 5e4, 5e5])
+    y_interp_truth = np.array([20.0343863828, 24.6395565688, 29.2447267548])
+    y_interp = log_interpolate_1d(x_interp, x_log, y_log)
+    assert_array_almost_equal(y_interp, y_interp_truth, 7)
+
+
 def test_log_interpolate_1d_units():
     """Test interpolating with log x-scale with units."""
     x_log = np.array([1e3, 1e4, 1e5, 1e6]) * units.hPa