Spatial indexing

[benbovy]

Currently a DGGSIndex only wraps a PandasIndex, which is limited for spatial indexing:

• dggs.sel_latlon converts the input latitude/longitude points as grid cells and then queries the pandas index of cells. If one of the lat/lon input points after being converted into a cell is not found in the index, the whole selection returns an error. This may be useful in some cases but not in (many) other case where we would want "true" nearest-neighbor lookup (i.e., find the nearest cell present in the index given any lat/lon point).
• The cell selection by polygon in #9 relies on the healpy.query_polygon function, which returns all grid cells (pixels) contained or intersecting the polygon, irrespective of whether or not those cells are actually present in the index. In #9 there's an extra step using numpy.isin() to avoid getting an error during the selection, which works but is probably not optimal. The more general issue is that other grid backends (H3) may not have this feature provided by healpy. Also, we might want more advanced query options (e.g., within vs. intersect, etc.).

A general solution for spatial indexing might be to follow the approach used in h3ron-polars, which provides a few spatial indexes (kd-tree, r-tree, packed-hilbert-rtree) that work with either the cell centroids or their envelope (rectangle). We could probably do the same here. We might be able to eventually reuse the indexes implemented in xoak once it is refactored to leverage Xarray custom indexes.

Alternatively, we could extract the cell geometries (#10) and then just rely on (or redirect users to) xvec for spatial queries. Xvec's GeometryIndex internally reuses GEOS' STRTree via shapely 2.0.

[danlooo]

I tried spatial indicies e.g. a k-d tree for the DGGRID DGGS in https://github.com/danlooo/DGGS.jl/issues/34 Searching for the nearest neighbor i.e. cell center lat/lon coordinate failed in reassembling the cell polygon geometry leading into many points that are assigned to a wrong cell. This is because the cell polygons are not regular anymore after the Snyder Equal Area Projection.