Congratulations on this great new feature!

Merged!

Thanks, @shoyer, @fmaussion, @jhamman, for the detailed review.

OK, let's merge this when CI passes.

Thanks, @shoyer.

I think the last thing to do here is remove keep_attrs?

Done! (I thought I reverted this change already...)

interp_like sounds a good idea. I like it in another PR as this is already huge.

Do you want to add interp_like here or should we leave it for a later PR?

Hmm. I guess numeric_only may currently only be an argument for the reduce method. We probably should add it for methods like mean, too.

Thanks, @shoyer

We do the same thing for aggregations like mean(), behavior originally copied from pandas: http://pandas.pydata.org/pandas-docs/version/0.22/generated/pandas.DataFrame.mean.html

Thanks. I have not been surprised when mean drops the object type array. So it would be OK also for interp.

Maybe this PR is ready, then.

BTW, I noticed that our mean() does not accept the numeric_only keyword argument. Our docs do not say about numeric_only keyword either. Is it intended? (This behavior looks already accepted, and I think it is OK if we do not support this.)

Do we have other methods that behave similarly?

We do the same thing for aggregations like mean(), behavior originally copied from pandas: http://pandas.pydata.org/pandas-docs/version/0.22/generated/pandas.DataFrame.mean.html

We could add a numeric_only=True argument (like we do for reduction methods) to make this behavior more obvious and allow it to be controlled.

Done!

But, I am slightly worrying if silently dropping non-numeric array is surprising. Do we have other methods that behave similarly?

Or maybe we can raise a warning?

I would wait before trying to do interpolation with datetime types. There might be use cases for this but I can't think of them offhand. resample() would seem to cover most of these.

On the other hand, it could indeed be interesting to extend resample() to handling numbers. I proposed this a while ago for pandas (https://github.com/pandas-dev/pandas/issues/12828).

1. I suppose it probably does make sense to interpolate numeric coordinates, too. So this probably doesn't need to change.

Agreed.

1. I do think we should drop non-numeric variables. There just isn't any way to meaningfully interpolate them (especially strings). This means we will probably need to add a check in DataArray.interp to ensure that the variable dtype is numeric and otherwise raise an error.

Also agreed. In theory, datetime dtype could be safely interpolated. Maybe we can support this if there is a demand.

It might be necessary to consider the same design issue also for interpolate_na, bfill, ffill.

My opinions: 1. I suppose it probably does make sense to interpolate numeric coordinates, too. So this probably doesn't need to change. 2. I do think we should drop non-numeric variables. There just isn't any way to meaningfully interpolate them (especially strings). This means we will probably need to add a check in DataArray.interp to ensure that the variable dtype is numeric and otherwise raise an error.

This is very elegant and I think very close to being ready to merge.

There are two design issues that I would like to discuss first: 1. Should we interpolate both data_vars and coords or only data_vars? 2. How should we handle non-numeric variables?

The current behavior of this PR is: 1. Interpolate all variables, both data_vars and coords. 2. Non-numeric variables raise an error, because they can't be converted into float.

It might make sense to copy the behavior from Dataset.reduce: 1. Only interpolate data_vars. 2. Drop non-numeric variables instead of trying to interpolate them.

Thanks for the careful review. All done.

Sorry, I still need to find time to finish reviewing this — hopefully very soon! On Wed, May 23, 2018 at 8:57 AM Joe Hamman notifications@github.com wrote:

Is there more to do here or is this ready to go in? Seems to me its ready to go.

— You are receiving this because you were mentioned.

Reply to this email directly, view it on GitHub https://github.com/pydata/xarray/pull/2104#issuecomment-391400628, or mute the thread https://github.com/notifications/unsubscribe-auth/ABKS1kp1VJkez0bTtOqUssKWUMc2tUIMks5t1YbOgaJpZM4Typql .

Is there more to do here or is this ready to go in? Seems to me its ready to go.

Very nice! I noticed that the interpolation is performed among dimensions rather than coordinates in this PR. However the limitation to that is interpolation to/from curvilinear grids is not supported, which I think is a common enough use case, and would be extremely nice to have. Pretty sure scipy's interpolation tools work out of the box with curvilinear grids. Is an updated interface which works on coordinate variables rather than dimensions planned?

Yes, our intent is to eventually extend this same interface to handle these cases, but not in this initial PR.

I think this is ready for final review

@fujiisoup That illustration is great. Thanks for all this great work. Yes, adding more information to that would work well.

I think it would be useful to have a table in the docs showing the differences between / when to use interp, interpolate_na, reindex, sel(..., method='nearest'), resample.

Agreed. Thanks for the nice suggestion :+1: We have discussed a cheat sheet #1552 previously, but no work has been done yet.

I added a simple illustration presenting the difference between advanced-sel and -interp for this addition, and I think we can enrich this kind of figure to clarify the differences among these set of methods.

I think it would be useful to have a table in the docs showing the differences between / when to use interp, interpolate_na, reindex, sel(..., method='nearest'), resample.

@fmaussion , thanks for the example! It's a nicely simple but descriptive example. Updated docs.

After your suggestion about attrs, I found a bug in apply_ufunc (#2114). The tests are added but now failing because of this issue.

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:

Thank you for the double checking :+1:

it would be great to implement the keep_attrs kwarg as well, because interpolation often conserves the units of data.

Agreed. Is there any problem if we always keep attrs in the interpolation? .sel always keeps attrs and I think it would be better if users have less choice in each method.

BTW, do you have any idea about an example for the new interp feature? It would be nice if we have a good-looking example on the doc.

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!!!!

Thanks, @fmaussion I didn't realize that scipy.interpolate.interpn does not sort the original coordinates (interp1d does). Thanks for pointing this out!

The implementation for you would be straightforward, you'd just have to pass np.arange(NN) to scipy instead of the actual coordinates.

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). Instead, I would like to call sortby in our interp routine, so that the array passed to scipy is always sorted in the ascending order.

The new API I would propose is python da.interp(x=[0, 1, 2], method='linear', assume_sorted=False, kwargs={'bounds_error': False, 'fill_value': 'extrapolate'}) with a new keyword assume_sorted.

@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.

As I think I've mentioned before, I think we should save interpolation of unstructured grids for later -- if only to keep the complexity of this PR to the minimum.

I agreed that it would be better to concentrate on structured interpolation only in this PR. I concerned about the interpolation along 1-dimensional but non-dimensional coordinate, which @fmaussion suggested in a review comment.

It is actually the structured interpolation, but currently, I do not support this just because of the analogy to .sel(). If we want to support this, we may need to decide its behavior in particular the dimension or coordinate of the interpolated array. I prefer not to support this now for keeping the interpolation rule simple.

should we support the interpolation along the non-dimensional coordinate?

As I think I've mentioned before, I think we should save interpolation of unstructured grids for later -- if only to keep the complexity of this PR to the minimum.

Interpolation on unstructured grids needs some of auxiliary data structure to be done efficiently, so it might make sense to wait until after we support flexible index types such as a KDTree index saved on xarray objects. That said, we could also add it now (immediately after this PR), and add support for caching KDTrees later.

I think the API for unstructured grid interpolation is pretty straightforward and could even reuse the same interp/interpolate_at method. Unstructured grid interpolation would be triggered whenever multiple existing coordinates that share some of the same dimensions are interpolated.

what should be the default value of 'bounds_error' and fill_value. scipy's interpolator assumes bounds_error=True by default, but I think bounds_error=False is more xarray-like, but it may be unintuitive for users very familiar with scipy's behavior.

I agree, by default we should switch to bounds_error=False, and replace extrapolated values with NaN, which is already the default fill value for interp1d.

Thanks, @fmaussion

The current design question is

• should we support the interpolation along the non-dimensional coordinate?
• what should be the default value of 'bounds_error' and fill_value. scipy's interpolator assumes bounds_error=True by default, but I think bounds_error=False is more xarray-like, but it may be unintuitive for users very familiar with scipy's behavior.

Thanks, @fmaussion. I think it is almost ready.

The remaining things are + I am looking for a good-looking example for our docs on rtd. + Just before the merge, I'm going to rename missing.py to interp.py so that the filename is more descriptive.

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.

So I think it might make sense to save skipna=True for interpolation on unstructured grids.

I agreed. I updated some of this PR.

The proposed API is python da.interp(x=[0, 1, 2], method='linear', kwargs={'bounds_error': False, 'fill_value': 'extrapolate'})

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.

My inclination would be to limit this method to interpolation on regular grids, i.e., the methods supported by scipy.interpolate.interpn or interp1d.

I think a skipna=True option could be useful, but I don't know how we could do that for regular interpolation with multi-dimensional data. So I think it might make sense to save skipna=True for interpolation on unstructured grids.

Yes, I am thinking this for multidimensional interpolation of nan-including arrays. (For 1-dimensional interpolation, we can interpolate them row by row) Although this should be another work (PR), I am wondering what the best API is.

Is skipna=True option handy? If so, which should be the default?

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)

RegularGridInterpolator handles NaNs pretty much the same way @shoyer describes it.

I did not certainly realize how RegularGridInterpolator works, but it seems to work with nan including array too as you pointed out. But higher order interpolations (RegularGridInterpolator is linear), such as RectBivariateSpline, would not work (return an all nan array) when just 1 nan is included in the target array. It is the same with higher order interpolations in scipy.interpolate.interp1d such as cubic.

Why don't leave the NaN handling part to the actual function doing the interpolation?

It is the original behavior of scipy.interpolate functions (except for linear or nearest methods), which do some matrix inversions in its algorithm.

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.

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.

@shoyer , Yes, linear and nearest can handle NaN intuitively, but other methods need to invert some matrices. For example, cubic spline can not be used with nan-including arrays, ```python In [1]: x = np.arange(4) ...: y = np.array([0, 1, np.nan, 3]) ...:

In [2]: # 'linear' works as expected ...: interp1d(x, y, kind='linear')([0.5, 1.5]) ...: Out[2]: array([0.5, nan])

In [3]: # 'cubic' returns an all nan array ...: interp1d(x, y, kind='cubic')([0.5, 1.5]) ...: Out[3]: array([nan, nan]) ```

I still need to look at the code here, but I don't think we need separate code paths for NaN vs no-NaNs. For example, suppose we have a DataArray with a coordinate x=[0, 1, 2], then da.interp(x=0.5, method='linear') should be equal to 0.5 * (da.sel(x=0) + da.sel(x=1)), regardless of the values of the DataArray. If either da[0] or da[1] is NaN, then the result would be NaN, too.

Maybe the possible API would be python da.interp(x=[0, 1, 2], skipna=True) for nan-skipping interpolation and python da.interp(x=[0, 1, 2], skipna=False) for faster interpolation.

@jhamman,

I just gave your implementation a quick skim and I think there is currently some duplication of efforts so maybe we look there to see what can be combined.

Agreed. I will try to generalize your interpolator-adaptor more.

@shoyer

perhaps we should give it a shorter name. Maybe interp?

Agreed. interp sounds nice.

Which NaN are you referring to? 1. data on the arrays being indexed 2. Coordinates on the arrays being indexed 3. Points at which to interpolate

I meant 1.

Scipy.interpolate.interp1d does not support arrays with NaN. So if we want to support this, we need to interpolate the array line by line as interpolate_na does, i.e. no vectorization. For the case of multidimensional interpolation of arrays with NaN, it is equivalent with the interpolation on non-grided data. It is possible, but the performance will significantly drop.

First of all, this is awesome!

One question: since I think this new interpolation method will be used quite frequently (more often than interpolate_na and reindex), perhaps we should give it a shorter name. Maybe interp?

We should also consider adding interp_like or interpolate_like by analogy to reindex_like.

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

Yes, I agree. Interpolation should be vectorized similarly to isel/sel.

In particular, in the long term I think we should aim to make interpolate_at(..., method='nearest') and reindex(..., method='nearest') equivalent.

How many interpolate methods should we support?

I think this is fine to start (we can always add more later).

How do we handle nan?

Which NaN are you referring to?

1. data on the arrays being indexed
2. Coordinates on the arrays being indexed
3. Points at which to interpolate

Case (1) should definitely be supported. Especially if data is stored in dask arrays, we cannot necessarily check if there are NaNs.

Cases (2) and (3) are not important, because there are relatively few use-cases for NaN in coordinate arrays.

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

Currently we don't support this for reindex(), but I suppose we could do so unambiguous.

How should we do if new coordinate has the dimensional coordinate for the dimension to be interpolated?

I think the new coordinates should take priority, and the dimension coordinates on the new coordinate should be dropped. This is similar to what we do for sel():

``` In [10]: da = xr.DataArray([0, 0.1, 0.2, 0.1], dims='x', coords={'x': [0, 1, 2, 3]}) ...:

In [11]: da.sel(x=xr.DataArray([1, 2], dims=['x'], coords={'x': [1, 3]})) ...: Out[11]: <xarray.DataArray (x: 2)> array([0.1, 0.2]) Coordinates: * x (x) int64 1 2 ```

@fujiisoup - This is great. I wonder if we can draw more from the interpolator adapters I built in missing.py:

https://github.com/pydata/xarray/blob/b9f40cc1da9c45b3dd33a3434b69c3d8fce57138/xarray/core/missing.py#L284-L323

I just gave your implementation a quick skim and I think there is currently some duplication of efforts so maybe we look there to see what can be combined.

@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...