Selecting point closest to (lon0, lat0) when lon,lat coordinates are 2D

[rsignell]

It's very common to want to extract a time series at a specified coordinate location, and I'm wondering whether xarray could support this directly without using xoak. Currently I'm using xoak, as in this reproducible example:

import xarray as xr
import intake
import xoak

intake_catalog_url = 'https://usgs-coawst.s3.amazonaws.com/useast-archive/coawst_intake.yml'
cat = intake.open_catalog(intake_catalog_url)
ds = cat['COAWST-USEAST'].to_dask()

lat, lon = 42.5, -70.0  # Gulf of Maine, 100km east of Boston, MA

da = ds['Hwave']
da.xoak.set_index(['lat_rho', 'lat_rho'], 'scipy_kdtree')
ds_point = xr.Dataset({"lon": ("point", [lon]), "lat": ("point", [lat])})
da.xoak.sel(lat_rho=ds_point.lat, lon_rho=ds_point.lon).sel(ocean_time='2012-10-01')

which produces:

Would it be possible to enable this in xarray or is this too specific a functionality to consider?

[welcome]

Thanks for opening your first issue here at xarray! Be sure to follow the issue template!
If you have an idea for a solution, we would really welcome a Pull Request with proposed changes.
See the Contributing Guide for more.
It may take us a while to respond here, but we really value your contribution. Contributors like you help make xarray better.
Thank you!

[keewis]

As far as I remember, we decided to not include a KDTree implementation directly in xarray (but I might be mistaken... it's hard to know exactly where to draw the line between "batteries included" and "third-party extension").

However, I believe @benbovy wanted to (at some point) change xoak to make use of the new index API. With that, your code snippet would become:

da = ds['Hwave'].set_xindex(["lat_rho", "lon_rho"], xoak.ScipyKDTree)
da.sel(lat_rho=ds_point["lat_rho"], lon_rho=ds_point["lon_rho"], time="2012-10-01")

or, using the existing API of xoak:

da = ds['Hwave'].xoak.set_index(["lat_rho", "lon_rho"], "scipy_kdtree")
da.sel(...)

[rsignell]

Thanks @keewis! That would be nicer, wouldn't it!?

[benbovy]

Yes, revamp xoak is part of the plan of my current work on indexes across the Xarray ecosystem.

Also addressing xarray-contrib/xoak#37, which will make your example nicer @rsignell.

[benbovy]

Given that Scipy is already an optional dependency of Xarray and assuming that a basic KDTreeIndex built on scipy.spatial.cKDTree wouldn't represent a big maintenance burden, it'd seem reasonable to me to have such index built in Xarray.

Xoak's scope is larger: it has experimental support for Dask and wraps more backend index structures from scikit-learn, pys2index, etc. In addition, Xoak could eventually have more geospatial features such as CRS-aware Xarray indexes built on top of Shapely or Spherely.