issues: 320275317

• [x] Closes #2079 (remove if there is no corresponding issue, which should only be the case for minor changes)
• [x] Tests added (for all bug fixes or enhancements)
• [x] Tests passed (for all non-documentation changes)
• [x] Fully documented, including whats-new.rst for all changes and api.rst for new API (remove if this change should not be visible to users, e.g., if it is an internal clean-up, or if this is part of a larger project that will be documented later)

I started working to add interpolate_at to xarray, as discussed in issue #2079 (but without caching).

I think I need to take care of more edge cases, but before finishing up this PR, I want to discuss what the best API is.

I would like to this method working similar to isel, which may support vectorized interpolation. Currently, this works as follwos

```python In [1]: import numpy as np ...: import xarray as xr ...: ...: da = xr.DataArray([0, 0.1, 0.2, 0.1], dims='x', coords={'x': [0, 1, 2, 3]}) ...:

In [2]: # simple linear interpolation ...: da.interpolate_at(x=[0.5, 1.5]) ...: Out[2]: <xarray.DataArray (x: 2)> array([0.05, 0.15]) Coordinates: * x (x) float64 0.5 1.5

In [3]: # with cubic spline interpolation ...: da.interpolate_at(x=[0.5, 1.5], method='cubic') ...: Out[3]: <xarray.DataArray (x: 2)> array([0.0375, 0.1625]) Coordinates: * x (x) float64 0.5 1.5

In [4]: # interpolation at one single position ...: da.interpolate_at(x=0.5) ...: Out[4]: <xarray.DataArray ()> array(0.05) Coordinates: x float64 0.5

In [5]: # interpolation with broadcasting ...: da.interpolate_at(x=xr.DataArray([[0.5, 1.0], [1.5, 2.0]], dims=['y', 'z'])) ...: Out[5]: <xarray.DataArray (y: 2, z: 2)> array([[0.05, 0.1 ], [0.15, 0.2 ]]) Coordinates: x (y, z) float64 0.5 1.0 1.5 2.0 Dimensions without coordinates: y, z

In [6]: da = xr.DataArray([[0, 0.1, 0.2], [1.0, 1.1, 1.2]], ...: dims=['x', 'y'], ...: coords={'x': [0, 1], 'y': [0, 10, 20]}) ...:

In [7]: # multidimensional interpolation ...: da.interpolate_at(x=[0.5, 1.5], y=[5, 15]) ...: Out[7]: <xarray.DataArray (x: 2, y: 2)> array([[0.55, 0.65], [ nan, nan]]) Coordinates: * x (x) float64 0.5 1.5 * y (y) int64 5 15

In [8]: # multidimensional interpolation with broadcasting ...: da.interpolate_at(x=xr.DataArray([0.5, 1.5], dims='z'), ...: y=xr.DataArray([5, 15], dims='z')) ...: Out[8]: <xarray.DataArray (z: 2)> array([0.55, nan]) Coordinates: x (z) float64 0.5 1.5 y (z) int64 5 15 Dimensions without coordinates: z ```

Design question

1. How many interpolate methods should we support? Currently, I only implemented scipy.interpolate.interp1d for 1dimensional interpolation and scipy.interpolate.RegularGridInterpolator for multidimensional interpolation. I think 90% usecases are linear, but there are more methods in scipy.

2. How do we handle nan? Currently this raises ValueError if nan is present. It may be possible to carry out the interpolation with skipping nan, but in this case the performance would be significantly drops because it cannot be vectorized.

3. Do we support interpolation along dimension without coordinate? In that case, do we attach new coordinate to the object?

4. How should we do if new coordinate has the dimensional coordinate for the dimension to be interpolated? e.g. in the following case, python da = xr.DataArray([0, 0.1, 0.2, 0.1], dims='x', coords={'x': [0, 1, 2, 3]}) rslt = da.interpolate_at(x=xr.DataArray([0.5, 1.5], dims=['x'], coords={'x': [1, 3]}) what would be rslt['x']?

I appreciate any comments.