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 add support for querying by geometry #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 add support for querying by geometry #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 danlooo/DGGS.jl#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.

[keewis]

I've come to the conclusion that I'd try as much as possible to avoid materializing the cell boundaries into memory: the xvec index currently requires loading all geometries into memory at the same time, which may exceed the memory (and even if xvec gets a distributed rtree, materializing polygons is costly).

Instead, how about defining DGGSIndex.sel such that:

• ds.sel(cell_ids=lat_lon_coords) does point-wise indexing (where the object is possibly a Coordinates or Dataset object? Not sure, could also be a custom dataclass)
• ds.sel(cell_ids=geometry, method="intersects") does polygon / multi-polygon selection (the method defines whether this should be intersection, strict inclusion, etc.)
• ds.sel(cell_ids=bbox_tuple) does bbox selection (could also define a specific Bbox dataclass, or introduce a type="spherical / type="planar" kwarg instead that allows switching between great / small-circle selection, or switch depending on whether the geometry is a shapely or spherely object); the method is the same as above

For the latter two, the geometry / bbox would be rasterized using the backend library (e.g. h3ronpy has geometry_to_cells), and then we'd just have to do a cell id selection on a pandas index while ignoring non-existent cell ids. If the backend library does not support this I think we can just create pull requests to make them support that.

It does feel a little awkward to do cell_ids=<geometry or coords>, but unless we'd want to allow cells=<geometry/coords> instead that may be unavoidable? Otherwise we could also define ds.dggs.sel as something like

def sel(self, geometry=None, /, *, method=None, **coords):
    pass

It also won't resolve the concern for method="nearest" above, but I don't think we have to resolve that to make progress on the above.

[benbovy]

A few thoughts:

• Storing explicit DGGS cell geometries is at a certain expense but not orders of magnitude more than storing the DGGS cell ids : “only” 6x more memory in the (worst?) case of hexagonal cells. When using a projected CRS it may require more vertices (tesselate) to accurately represent curved cell edges, although this quickly becomes irrelevant for higher refinement levels (small cells). Otherwise cell corners is enough I guess with tools that work with geodetic edges (e.g. spherely).

• Spherely (via s2geometry) might eventually provide memory-optimized & serializable spatial indexes (see Querying geographies (spatial index) benbovy/spherely#72), although this will require some work.

• I think that in many cases Xvec will come handy (although sub-optimal) for spatial selection of DGGS cells using arbitrary input geometries. I also think that it might be better to support it transparently, i.e., provide some convenient API here to convert a DGGS data cube to a vector data cube -- or simply to add cells as POLYGON or POINT indexes geometry coordinates -- and let users know they can do spatial queries via Xvec. Perhaps it would be better to have concrete benchmark results or user feedback before trying to implement more optimal and flexible spatial indexing here?

• I'm not sure that rasterizing arbitrary geometries to DGGS cells is a trivial operation (see h3ronpy's polygon fill modes). It may not be always the most optimized approach either, e.g., rasterizing a large bbox at a high (fixed) refinement level would result in a very large number of cells.

• DGGS data cubes may have two main configurations: either a "raster-like" coverage of contiguous cells that are already sorted following a certain pattern (that may have some spatial meaning, e.g., healpix's RING vs. NESTED) or a "point-cloud-like" bunch of non-sorted, sparsely distributed cells (I'm here skipping "multi-resolution" / compact-form cell coverage for simplicity). The best approach for spatial indexing will likely differ depending on which of these configurations is considered.

[benbovy]

Regarding the API, I would avoid trying to support all kinds of data selection via Xarray’s sel() API.

I find ds.sel(cell_id=<cell ids | lat-lon coords | geometries>) potentially very confusing! (on a more general note I'm a bit concerned already by all kinds of behaviors that will result from more custom indexes available in the broader Xarray ecosystem). In fact I would restrict here the indexer type to actual cell ids, i.e., ds.sel(cell_id=<cell ids>), and why not support ids given on different refinement levels sometime in the future (e.g., select all cells given a cell ids at a parent level).

For more advanced queries, I’d rather see it exposed via something like ds.dggs.query(…) similarly to ds.xvec.query(…).

[keewis]

“only” 6x more memory in the (worst?) case of hexagonal cells

I think that's 12x: cell ids are one 64-bit integer per cell, while boundaries have 5-6 corners per cell with 2 64-bit floats per corner. For the case of healpix, this is 8x (4 corners with 2 floating point values each), which can make the difference between "comfortably fits into memory on a medium sized laptop" and "requires bigger hardware to even get created" (for healpix level 13, this is 6.4GB compared to 51.5GB). And that estimate is a lower bound because it doesn't count the overhead of the polygon objects / RTree.

might eventually provide memory-optimized & serializable spatial indexes

I guess I didn't really point this out, but since the DGGS cells are already arranged in a tree I think it should be possible to make use of that? Not sure how geometry rasterization works for most DGGS.

Either way the RTree will also have to provide the different containment / intersection modes, so there's nothing we would gain here.

rasterizing a large bbox at a high (fixed) refinement level would result in a very large number of cells.

indeed, but this can most likely be optimized, for example by extending the backend library to only return candidate cell ids if they are part of an input set. Or maybe my suggestion was just bad and there's a better way to search cell ids in a DGGS?

Perhaps it would be better to have concrete benchmark results or user feedback before trying to implement more optimal and flexible spatial indexing here?

That's the current state, which has turned out to be somewhat slow (see e.g. the healpix subsetting in the data interpolation script) compared to rasterization (the h3 subsetting in the same script)

However, while I think we can do better for searching, I'm very open to conversion methods (ds.dggs.to_xvec()?), and I'm also trying to push for a way to allow converting to uxarray in UXARRAY/uxarray#1055.

For more advanced queries, I’d rather see it exposed via something like ds.dggs.query(…) similarly to ds.xvec.query(…)

My inspiration was xvec's hook into ds.sel as that also only allows a single geometry as input, but I can see how cell_ids=geometry can be confusing. As mentioned above, I'd be fine with a separate method, though, whatever the name.

[benbovy]

I think that's 12x

Fair enough! (6x if 32-bit floating point precision is enough :-) ).

I guess my main point is that, although far from optimal (but probably not that bad either in many cases), converting DGGS as vector geometries and then use vector tools for spatial queries based on input arbitrary geometries is a practical solution that can work today and that would provide consistent results for any kind of DGGS. It also works with distance-based queries (find nearest cells).

The opposite approach (i.e., convert the query input geometry into DGGS cells + hash-table lookup) does certainly make sense too and might be better at indexing large datasets. The downside is that to my knowledge the current DGGS implementations wouldn't allow us supporting this approach in a consistent & DGGS-agnostic way. Also I don't think that extending those backend libraries would be an easy task (especially if we want to avoid doing so in a low-level language). That said, probably it is OK supporting this approach via xdggs even though details may vary from one DGGS to another and/or it is only for a subset of the DGGSs available? I.e., Xdggs would just provide a uniform API but would not guarantee a uniform behavior across DGGSs.

Both approaches are complementary IMO.

That's the current state, which has turned out to be somewhat slow

How does it scale? Do you have some results or plots we can look at?

My inspiration was xvec's hook into ds.sel as that also only allows a single geometry as input

Yeah, for xvec it makes more sense IMO since xvec index coordinate labels are geometry objects.

[keewis]

Do you have some results or plots we can look at?

Nothing systematic, no. If we need a benchmark I can do that once I have a bit more time, though.

Both approaches are complementary IMO.

Agreed. And additionally, going through xvec works today, while there's a lot more work necessary to support the other approach.

For another idea – which may have been mentioned before – is selecting by parent cell ids (and this one ds.sel(cell_ids=parent_cell_id) actually does make sense). For this to work, we need to figure out how to pass both the level and the cell id at the same time, or require an indexing scheme that combines both.

[jbusecke]

Hey folks, getting more into xdggs over here at the km-scale hackathon

Warning

I am pretty much a noob with anything geometry/unstructured, so def tell me if I am talking nonsense!

I was just talking to @surgura and @VeckoTheGecko who wan to work on this. I basically wanted to chime in here with some thoughts, questions and a WIP schematic of how this would work on a high level.

Ok diagram first:

So I wonder if things can be roughly abstracted into two stages:

• Constructing a geometry in the appropriate crs of the data (e.g. healpix vs h3) from a variety of inputs:
  • Some polygon input (which might need a source_crs input?)
  • Some notion of a lat/lon box (this was my original entry point to thinking about that. I was hoping that we could do something like ds.dggs.sel(lon=slice(0,10), lat=slice(-30,-10)), but I think @benbovy s comment makes a good point, so Ill think about .query() instead. I somehow feel that this would be more elegant to have here rather than in dggs.sel_latlon?
  • Some way of getting the parend cell id
  • I did not really understand the input type listed first in the post by @keewis comment
  • My assumption here is that all of these inputs can ultimately be converted to a single format (a geometry in the appropriate crs), which would nicely separate concerns from the following step
• Filtering cells based on input geometry
  • I think that it is fairly clear that we might need different methods here, but if we can internally assume that each method will accept only a single type of input, that might reduce complexity?

Curious to hear if this is completely off base or if yall think this is conceptually ok

On a more practical (hackathon) level. I am not sure how much I can contribute to the methods here, but I thought it might be useful to start building some high level test cases that construct synthetic inputs/desited outputs for a variety of input types and methods?

[surgura]

FYI @VeckoTheGecko and me are specifically looking at #134
There I created a little summary of what we plan ✌️

[keewis]

Some comments on the cases in the graph (this is easy to understand, thanks!):

• All geometries are assumed to be in spherical coordinates (for now), and would need to be reprojected manually or automatically (I'd tend towards manually)
• I tend to think of lat / lon slices as bounding boxes, so these could very well be expressed as a polygon. With that, you'd be left with a single generic case of selecting by polygon with two options on how to interpret the polygon edges: great circles, or rhumb lines / straight lines in a equirectangular projection (I'll call that edge_mode below).
• I'd very much like to avoid implementing the selection by parent cell ids by going through polygons: it is much easier and faster to figure out the parent cell id of each cell by dividing the cell ids in nested ordering by 4^delta_level and then comparing integers (for healpix, but for h3 there's h3ronpy.change_resolution)

Such that in the end, you'd have three methods for selecting:

# by cell id / parent cell id (uses integer operations and a pandas index or similar to select)
ds.sel(cell_ids=cell_ids | parent_cell_ids)
# point-wise indexing (converts to cell ids / containing cell and forwards to `Dataset.sel`)
ds.dggs.sel_latlon(lon=longitudes, lat=latitudes)  # may call that `.dggs.sel`? or `sel_points`?
# by geometry
ds.dggs.query(geometry, method=containment_mode, edges=edge_mode)

[benbovy]

Some notion of a lat/lon box (this was my original entry point to thinking about that. I was hoping that we could do something like ds.dggs.sel(lon=slice(0,10), lat=slice(-30,-10)), but I think @benbovy s #16 (comment) makes a good point, so Ill think about .query() instead. I somehow feel that this would be more elegant to have here rather than in dggs.sel_latlon?

I agree this would be more elegant. One way to support this I think could be to reuse Xarray's new and flexible CoordinateTransform (pydata/xarray#9543) to handle the conversion between DGGS cells and lat/lon coordinates (cell centroids) or geometry objects (cell polygons). This would allow exposing the two latter as Xarray lazy coordinates. We would also need to wrap everything into the same Xarray index, e.g., something like

>>> ds
<xarray.Dataset>
Dimensions:      (cell: 50000)
Coordinates:
  * cell_id       (cell) int64 ...
  * lat           (cell) float64 ...
  * lon           (cell) float64 ...
  * geom          (cell) object ...
Indexes:
  ┌ cell_id        DGGSIndex
  | lat
  | lon
  └ geom

where DGGSIndex becomes a "meta-index" that here encapsulates a PandasIndex for cell_ids and two CoordinateTransformIndex instances for (lat, lon) and geom respectively.

This makes possible various kinds of selections like ds.sel(cell_id=1092838), ds.sel(lon=slice(0,10), lat=slice(-30,-10)) or ds.sel(geom=shapely.bbox(...)), which are all internally handled by DGGSIndex.

More generally, I think that the cell_id coordinate is rarely meaningful for end-users and it is more useful for them to have the lat, lon and/or geom coordinates around. This solutions allows having those at hand without materializing the cell centroids and boundaries into memory.