interpolate.log_interpolate_1d fails with Xarray input of more than 3D (Numpy succeeds)

[erogluorhan]

• 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

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)