issue_comments

@shoyer I'm having trouble wrapping my head around this. The example above is essentially a 1D interpolation over y, a coordinate selection over x, and a broadcast over z. I don't think it should interpolate over x, right ?

I recently wrote a version of scipy.ndimage.map_coordinates for JAX in pure NumPt that I think could be straightforwardly ported into xarray.

On Thu, Feb 6, 2020 at 7:00 AM David Huard notifications@github.com wrote:

Just got bit by this as well. Computing monthly quantile correction factors, so I have an array with dimensions (month, quantile, lon, lat). I then want to apply these correction factors to a time series (time, lon, lat), so I compute the month and quantile of my time series, and want to interp into the quantile correction factors. This doesn't work because both the factors and the time series have (lat, lon) dimensions.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/pydata/xarray/issues/3252?email_source=notifications&email_token=AAJJFVWXXHT6ISDFN6ETRO3RBQQ2VA5CNFSM4IPBQQPKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEK7QWSA#issuecomment-582945608, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAJJFVUF4COJXNL377B3FZDRBQQ2VANCNFSM4IPBQQPA .

Just got bit by this as well. Computing monthly quantile correction factors, so I have an array with dimensions (month, quantile, lon, lat). I then want to apply these correction factors to a time series (time, lon, lat), so I compute the month and quantile of my time series, and want to interp into the quantile correction factors. This doesn't work because both the factors and the time series have (lat, lon) dimensions.

Ok. I realized this problem occurs only because x was a dimension of the new index idx and the original dataset. Perhaps sel should warn the user or raise an error when this happens.

Ok. I started playing around with this, but I am getting errors when indexing arrays with ND variables. data.sel(x=nd_x) works, but any subsequent operations complain that IndexVariable objects must be 1-dimensional: ```

import xarray as xr import numpy as np npdata = np.tile(np.arange(10), (5, 1)) ... data = xr.DataArray(npdata, dims=['y', 'x'], ... coords={'x': np.r_[:10], 'y': np.r_[:5]}) ... idx = xr.DataArray(npdata, dims=['z', 'x']) ... ... ans = data.sel(x=idx, method='bfill') ... assert set(ans.dims) == {'z', 'y', 'x'} ... print(ans) Traceback (most recent call last): File "<ipython-input-4-1d0ac0cac680>", line 8, in <module> print(ans) File "/Users/noah/workspace/software/xarray/xarray/core/common.py", line 129, in repr return formatting.array_repr(self) File "/Users/noah/workspace/software/xarray/xarray/core/formatting.py", line 463, in array_repr summary.append(repr(arr.coords)) File "/Users/noah/workspace/software/xarray/xarray/core/coordinates.py", line 78, in repr return formatting.coords_repr(self) File "/Users/noah/workspace/software/xarray/xarray/core/formatting.py", line 381, in coords_repr coords, title="Coordinates", summarizer=summarize_coord, col_width=col_width File "/Users/noah/workspace/software/xarray/xarray/core/formatting.py", line 361, in _mapping_repr summary += [summarizer(k, v, col_width) for k, v in mapping.items()] File "/Users/noah/workspace/software/xarray/xarray/core/formatting.py", line 361, in <listcomp> summary += [summarizer(k, v, col_width) for k, v in mapping.items()] File "/Users/noah/workspace/software/xarray/xarray/core/formatting.py", line 307, in summarize_coord coord = var.variable.to_index_variable() File "/Users/noah/workspace/software/xarray/xarray/core/variable.py", line 440, in to_index_variable self.dims, self._data, self._attrs, encoding=self._encoding, fastpath=True File "/Users/noah/workspace/software/xarray/xarray/core/variable.py", line 1943, in init raise ValueError("%s objects must be 1-dimensional" % type(self).name) ValueError: IndexVariable objects must be 1-dimensional ```

We could probably use this our own version for all linear and nearest neighbor interpolation. Then we won’t need scipy installed for that.

So when would interp use this manual interpolation method? Would it try it only if scipy fails or check the dimensions?

After reading the method='ffill' docs, I agree with you.

Linear interpolation for 1d -> nd is just a matter of averaging two indexing selections. If we leverage xarray's vectorized indexing operations to do the hard work, it should work automatically for dask arrays, sparse arrays and xarray's internal backend array types, all with any number of dimensions.

On Fri, Aug 23, 2019 at 2:17 PM Noah D Brenowitz notifications@github.com wrote:

In my experience, computing w efficiently is the tricky part. The function is slightly different, but metpy https://unidata.github.io/MetPy/latest/api/generated/metpy.interpolate.interpolate_1d.html uses a lot of tricks to make this work efficiently. A manual for-loop is much cleaner for this kind of stencil calculation IMO. What kind of duck arrays were you thinking of?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/pydata/xarray/issues/3252?email_source=notifications&email_token=AAJJFVQMRK7YOBW6LGJSXM3QGBAUHA5CNFSM4IPBQQPKYY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOD5BHDBY#issuecomment-524448135, or mute the thread https://github.com/notifications/unsubscribe-auth/AAJJFVUICIWBNW3EILTCQ53QGBAUHANCNFSM4IPBQQPA .

In my experience, computing w efficiently is the tricky part. The function is slightly different, but metpy uses a lot of tricks to make this work efficiently. A manual for-loop is much cleaner for this kind of stencil calculation IMO. What kind of duck arrays were you thinking of?

We could implement linear interpolation just as w * array.sel(..., method='ffill') + (1 - w) * array.sel(..., method='bfill'), where w is calculated from the indexes. This is probably about as efficient as scipy. Numba could be slightly faster (by avoiding a few memory allocations) but it will be hard to make it work as generally, especially on duck arrays.

I have some numba code which does this for linear interpolation. Does scipy support this pattern?