Surprising behaviour of Dataset/DataArray.interp() with NaN entries

I think this is due to documented 'undefined behaviour' of scipy.interpolate.interp1d, so not really a bug, but I think it would be more user-friendly if xarray gave an error in this case rather than producing an 'incorrect' result.

What happened:

If a DataArray contains a NaN value and is interpolated, output values that do not depend on the entry that was NaN may still be NaN.

What you expected to happen:

The docs for scipy.interpolate.interp1d say

Calling interp1d with NaNs present in input values results in undefined behaviour.

which explain the output below, and presumably mean it is not fixable on the xarray side (short of some ugly work-around). I think it would be good though to check for NaNs in DataArray/Dataset.interp(), and if they are present raise an exception (or possibly a warning?) about 'undefined behaviour'.

scipy.interpolate.interp2d has a similar note, while scipy.interpolate.interpn does not mention it (but has very limited information).

What I'd initially expected was an output would be valid at locations in the array that shouldn't depend on the NaN input: interpolating a 2d DataArray (with dims x and y) in the x-dimension, if only one y-index in the input has a NaN value, that y-index in the output might contain NaNs, but the others should be OK.

Minimal Complete Verifiable Example:

import numpy as np
import xarray as xr

da = xr.DataArray(np.ones([3, 4]), dims=("x", "y"))

da[0, 0] = float("nan")

newx = np.linspace(0., 3., 5)

interp_da = da.interp(x=newx)

print(interp_da)

On my system, this gives output:

<xarray.DataArray (x: 5, y: 4)>
array([[nan,  1.,  1.,  1.],
       [nan,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.],
       [nan, nan, nan, nan],
       [nan, nan, nan, nan]])
Coordinates:
  * x        (x) float64 0.0 0.75 1.5 2.25 3.0
Dimensions without coordinates: y

[Surprisingly, I get the same output even using method="nearest".]

You might expect at least the following, with NaN only at y=0:

<xarray.DataArray (x: 5, y: 4)>
array([[nan, 1., 1., 1.],
       [nan, 1., 1., 1.],
       [ 1., 1., 1., 1.],
       [nan, 1., 1., 1.],
       [nan, 1., 1., 1.]])
Coordinates:
  * x        (x) float64 0.0 0.75 1.5 2.25 3.0
Dimensions without coordinates: y

Environment:

Start by running the provided DataArray.interp() example and trace the DataArray/Dataset.interp() path into the scipy.interpolate.interp1d behavior described in the issue. Done means NaN-containing inputs no longer silently produce surprising results, with the chosen error or warning behavior covered by a regression test.