Is there any problem if we always keep attrs in the interpolation?

In my opinion it should be the default, yes.

BTW, do you have any idea about an example for the new interp feature?

A quick shot would be to use the existing sample dataset:

```python import matplotlib.pyplot as plt import xarray as xr import numpy as np

Raw data

ds = xr.tutorial.load_dataset('air_temperature') ds.air.isel(time=0).plot() plt.title('Raw data')

Interpolated data

new_lon = np.linspace(ds.lon[0], ds.lon[-1], ds.dims['lon'] * 4) new_lat = np.linspace(ds.lat[0], ds.lat[-1], ds.dims['lat'] * 4) dsi = ds.interp(lon=new_lon, lat=new_lat) dsi.air.isel(time=0).plot() plt.title('Interpolated data') ``` Which produces the plots:

As the interpolation routine can also be used for non-uniform gridded data, I don't think passing np.arange(NN) to scipy is a good idea (it will change the result if higher order method such as 'cubic' is used).

Right. Another feature of geospatial grids is that they are regularly spaced, which is not always the case here. Forget about my idea for now.

I tested your updated branch and we now come to the same results with my test data! I took the liberty to edit your comment because in your current implementation assume_sorted is False per default.

Another wish (sorry! ;): it would be great to implement the keep_attrs kwarg as well, because interpolation often conserves the units of data. Thanks!!!!

@fujiisoup I gave a go at your branch this evening, trying a few things on real data cases. First of all: this is great! Really looking forward to see this in xarray.

I noticed a small issue: currently this doesn't work with decreasing coordinates:

python nx = 5 ny = 4 lon = np.arange(nx) lat = np.arange(ny)[::-1] data = np.arange(nx*ny).reshape((ny, nx)) da = xr.DataArray(data, dims=('lat', 'lon'), coords={'lon':lon, 'lat':lat}) da.interp(lon=[1.1, 2.2, 3.3], lat=[1.2, 2.3]) -> ValueError from SciPy: The points in dimension 0 must be strictly ascending Traceback:

Decreasing coordinates are really frequent in climate data (latitudes are often decreasing, don't ask me why), so I wondered how I implemented this in salem, which brings me to the second point: I think we need an interp() equivalent to isel() too, which is the one I could make use of with salem (and, maybe, several other georeferencing tools).

The way salem works is that it uses the concept of geospatial grid where the points to interpolate are defined in a local coordinate system which is always positive and increasing. That is, if I want to use xarray's interpolation routines, I'd like to be able to tell xarray that I'd like to interpolate at "indices [1.1, 2.1, 3.1]" for example, where these can really be understood as indexes in the way da.isel() sees them.

I don't know if we really need a i_interp for this though, or if it could be solved with a kwarg to the existing API. The implementation for you would be straightforward, you'd just have to pass np.arange(NN) to scipy instead of the actual coordinates.

Thoughts?

I prefer not to support this now for keeping the interpolation rule simple.

Agreed. Users can turn their 1D non dimension coordinates into dimension coordinates if they want to.

Let us know when this is ready for review! I'm happy to have a look at it.

Honestly, I am not very happy with kwargs={} arguments (it would be more beautiful if they can be keyword arguments), but I think it is a necessary pain

I don't think it's that bad. It makes it clear that the arguments are passed over to the underlying function, and that xarray is not doing all the magic by itself.

I think the addition of the linear interpolation is still nice, but it would be nicer if we could also implement nan-skipping interpolations, which will enable us to use higher order interpolation with nan-including arrays.

Do you mean you want to use unstructured data interpolation routines? (https://docs.scipy.org/doc/scipy/reference/interpolate.html#multivariate-interpolation)

linear and nearest can handle NaN intuitively, but other methods need to invert some matrices

Why don't leave the NaN handling part to the actual function doing the interpolation? RegularGridInterpolator handles NaNs pretty much the same way @shoyer describes it.

@fujiisoup this is awesome! So glad you are looking into this.

To point 1: How many interpolate methods should we support?.

I think that your implementation is a very good first step. From the geosciences point of view I think that Spline interpolation would be a very good candidate too, but it only works in the 2-dimensional case which makes it a no-go for the general case...