issue_comments: 121679580

So, the good news is that once we figure out the API for pointwise indexing, I think the nearest-neighbor part could be as simple as supplying method='nearest'.

The challenge is that we want to go from an DataArray that looks like this:

``` In [4]: arr = xray.DataArray([[1, 2], [3, 4]], dims=['x', 'y'])

In [5]: arr Out[5]: <xray.DataArray (x: 2, y: 2)> array([[1, 2], [3, 4]]) Coordinates: * x (x) int64 0 1 * y (y) int64 0 1 ```

To one that looks like that:

In [6]: xray.DataArray([1, 4], {'x': ('c', [0, 1]), 'y': ('c', [0, 1])}, dims='c') Out[6]: <xray.DataArray (c: 2)> array([1, 4]) Coordinates: y (c) int64 0 1 x (c) int64 0 1 * c (c) int64 0 1

Somehow, we need to figure out the name for the new dimension (c in this example).

My thought would be to have methods sel_points and isel_points that work similarly to sel and isel. This is straightforward if you already have xray 1D objects with a labeled dimension: arr.sel_points(x=x, y=y), where x and y are along the c dimension.

If you don't already have 1D xray objects, I suppose we could also allow arr.sel_points(x=('c', [0, 1]), y=('c', [0, 1])) or arr.sel_points('c', x=[0, 1], y=[0, 1]).