issue_comments: 241821366

A few recent developments relevant to this issue: - #974 discusses how we could add multi-dimensional indexing with broadcasting. This would subsume the need for separate methods like sel_points and allow also handle indexing grids with grids. - #947 adds first class support for MultiIndex coordinates into xarray. This is good model for how a KDTree could work.

So I'm now thinking an API more like this:

```

ds = ds.set_kdtree(spatial_index=['latitude', 'longitude'])

ds <xray.Dataset> Dimensions: (x: 4, y: 5) Coordinates: * x (x) int64 0 1 2 3 * y (y) int64 0 1 2 3 4 * spatial_index (x, y) KDTree - latitude (x, y) float64 0.49 0.5682 -0.3541 -0.9305 -0.9669 0.01558 ... - longitude (x, y) float64 0.3758 1.429 -1.698 -1.344 0.5237 0.6152 ... Data variables: temperature (x, y) float64 0.5735 -0.4871 0.4708 0.4907 -0.3318 0.2883 ...

result = ds.sel(latitude=other.latitude, longitude=other.longitude, ... method='nearest') ```

For building a tree with lat/lon remapped to spherical coordinates, we should write a method that converts lat and lon arrays into a tuple of x, y, z arrays (e.g., using apply_ufunc from #964). Then this looks like ds.set_kdtree(spatial_index=latlon_to_xyy(ds.latitude, ds.longitude)). Conceivably, we could add some sugar for this, e.g., ds.geo.set_kdtree(spatial_index=['latitude', 'longitude']).