Multi-dimensional extrapolation doesn't work

[lstngr]

What happened?

I am trying to interpolate and extrapolate a two dimensional DataArray. Passing the kwargs=dict(fill_value=None) argument to DataArray.interp still returns nan values.

The MCVE is taken from the documentation guide on interpolation. (Note that the documentation also shows nan values.)

What did you expect to happen?

I expected the routine to extrapolate the data when "out-of-bounds" coordinates were provided.

Minimal Complete Verifiable Example

import xarray as xr
import numpy as np

# multi-dimensional extrapolation
da = xr.DataArray(
     np.sin(0.3 * np.arange(12).reshape(4, 3)),
    [("time", np.arange(4)), ("space", [0.1, 0.2, 0.3])],
)
interped_da = da.interp(time=4, space=np.linspace(-0.1, 0.5, 10), kwargs={"fill_value": None})
print(interped_da)

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.

Relevant log output

<xarray.DataArray (space: 10)>
array([nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])
Coordinates:
    time     int64 4
  * space    (space) float64 -0.1 -0.03333 0.03333 0.1 ... 0.3 0.3667 0.4333 0.5

Anything else we need to know?

This is my first issue, I hope I didn't forget something! Thank you for the amazing package!

Environment

INSTALLED VERSIONS
------------------
commit: None
python: 3.9.12 | packaged by conda-forge | (main, Mar 24 2022, 23:22:55) 
[GCC 10.3.0]
python-bits: 64
OS: Linux
OS-release: 5.3.18-150300.59.63-default
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8
LOCALE: ('en_US', 'UTF-8')
libhdf5: 1.10.6
libnetcdf: None

xarray: 2022.3.0
pandas: 1.4.2
numpy: 1.21.5
scipy: 1.7.3
netCDF4: None
pydap: None
h5netcdf: None
h5py: 3.6.0
Nio: None
zarr: None
cftime: None
nc_time_axis: None
PseudoNetCDF: None
rasterio: None
cfgrib: None
iris: None
bottleneck: 1.3.4
dask: 2022.05.0
distributed: 2022.5.0
matplotlib: 3.5.2
cartopy: None
seaborn: None
numbagg: None
fsspec: 2022.3.0
cupy: None
pint: None
sparse: None
setuptools: 62.2.0
pip: 22.1
conda: None
pytest: None
IPython: 8.3.0
sphinx: 4.5.0

[keewis]

try da.interp(time=4, space=np.linspace(-0.1, 0.5, 10), kwargs={"fill_value": "extrapolate"}).

If I'm understanding the purpose of that example correctly that's a documentation bug, and we would welcome a PR to fix it.

[lstngr]

You're right, thanks for the help. It seems that orthogonal dimensions are now interpolated independently by calling SciPy's interp1d multiple times, which explains why passing 'extrapolate' works.
This behaviour has been introduced in #4155, in particular 7daad4f.

The documentation should be adapted to use @keewis's approach. I'd suggest including an example combining "advanced interpolation" and extrapolation, such as the one below:

da = xr.DataArray(
    np.sin(0.3 * np.arange(20).reshape(5, 4)),
    [("x", np.arange(5)), ("y", [0.1, 0.2, 0.3, 0.4])],
)
x = xr.DataArray([0.5, 1.5, 2.5, 3.5], dims="z")
y = xr.DataArray([0.15, 0.25, 0.35, 0.45], dims="z")

Without extrapolation:

da.interp(x=x, y=y)

<xarray.DataArray (z: 4)>
array([ 0.55626357,  0.63496063, -0.46643289,         nan])
Coordinates:
    x        (z) float64 0.5 1.5 2.5 3.5
    y        (z) float64 0.15 0.25 0.35 0.45
Dimensions without coordinates: z

With extrapolation:

da.interp(x=x, y=y, kwargs={'fill_value': None})

<xarray.DataArray (z: 4)>
array([ 0.55626357,  0.63496063, -0.46643289, -0.73507668])
Coordinates:
    x        (z) float64 0.5 1.5 2.5 3.5
    y        (z) float64 0.15 0.25 0.35 0.45
Dimensions without coordinates: z