odd multidimensional extrapolation behaviour with interp

What happened?

When using interp on a 3-dimensional DataArray to extract values on a 2-dimensional surface, there is no extrapolation along one of the dimensions of the 2D plane (but the other dimension works). In the MWE, I have created a 5x5x5 cube centred on the origin where x, y, and z run from -2 to 2, then interpolated along a 6x6 plane that is perpendicular to the y-axis, where the plane extends beyond the faces of the cube by 0.5 units. The values for which x are outside the cube extrapolate correctly, however the values for which z are outside the cube become nan.

What did you expect to happen?

Extrapolation should happen along both dimensions of the plane.

Minimal Complete Verifiable Example

>>>import numpy as np
>>>import xarray as xr
>>>x = np.linspace(-2,2,5)
>>>da = xr.DataArray(np.random.rand(5,5,5), coords={"x": x, "y": x, "z": x})
>>>xi = xr.DataArray(np.linspace(-2.5,2.5,6), dims="theta")
>>>yi = xr.DataArray(0.5*np.ones(6), dims="theta")
>>>zi = np.linspace(-2.5,2.5,6)
>>>da.interp(x=xi, y=yi, z=zi, kwargs={"fill_value": None})
<xarray.DataArray (theta: 6, z: 6)> Size: 288B
array([[       nan, 0.65509617, 0.65540961, 0.94953277, 0.96035659,
               nan],
       [       nan, 0.50551769, 0.50022194, 0.55698788, 0.47609523,
               nan],
       [       nan, 0.53349845, 0.55660868, 0.45080263, 0.36583418,
               nan],
       [       nan, 0.6058017 , 0.62612734, 0.49364122, 0.52795018,
               nan],
       [       nan, 0.67642216, 0.63931782, 0.56366723, 0.52998289,
               nan],
       [       nan, 0.87859657, 0.79462258, 0.79821649, 0.47355559,
               nan]])
Coordinates:
    x        (theta) float64 48B -2.5 -1.5 -0.5 0.5 1.5 2.5
    y        (theta) float64 48B 0.5 0.5 0.5 0.5 0.5 0.5
  * z        (z) float64 48B -2.5 -1.5 -0.5 0.5 1.5 2.5
Dimensions without coordinates: theta

MVCE confirmation

• Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• Complete example — the example is self-contained, including all data and the text of any traceback.
• Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• New issue — a search of GitHub Issues suggests this is not a duplicate.
• Recent environment — the issue occurs with the latest version of xarray and its dependencies.

Relevant log output

No response

Anything else we need to know?

No response

Environment

No source files or tests are named in the issue. Reproduce the supplied MWE with the reported xarray, NumPy, and SciPy versions, then locate the interpolation implementation and relevant tests; done means extrapolation works for both out-of-range plane dimensions without changing in-range results.