Accessor function that accept polygon and return those cells (with data) that are within the area ?

[tik65536]

Hey hi ~

Would it be nice to have a function under the dggs accessor that accepts shapely polygon objects (or geojson) and returns data within the polygon (convert the polygon to cells and sel the data by the id) ?

Thanks,
Tik

[keewis]

That would be nice, indeed. It shouldn't be that hard to do, but we may need to extend the backend libraries (cdshealpix, h3ronpy) to perform the search for us: in an earlier attempt I've tried to "rasterize" a polygon, then select based on the resulting cell ids (see #9).

This does work, but is really inefficient: imagine you're trying to subset a high-res regional dataset with a huge polygon. The rasterization will give you all the cell ids within the polygon (for some definition of "within"), but we're really only interested in a small subset so most will be thrown away. What we'd want is something like this:

cdshealpix.nested.polygon_select(level, cell_ids, geom)

(and something similar in h3ronpy)

With that, we could provide the query API that was briefly discussed in #16

[surgura]

@VeckoTheGecko and me are looking into this, at least during the WCRP global km scale hackaton (2025).

We identified the following possible definitions of 'in the polygon':

1. center in the shape
2. full cell in shape
3. cell partially in shape

Could also provide:
a. ratio of how much of the cell is in the shape for 1,3
b. shapes cells split for 1,3

Order of interest
1 > 3 > 2
a: nice to have
b: I presume very low interest

Will first look into how to implement this efficiently, and were and what and if in the rust libraries.
See also
#9
https://github.com/keewis/healpix-geo -> soon a dependency of xdggs, so a good place to add functionality on top on healpix, in rust.

[surgura]

Ok. Before I did not fully grasp the problem yet, but now I understand it is also important to take into account that not the complete healpix sphere has data (usually people work with a specific region). It' is possible for a user to create a polygon that is partly out of the area with data. While this is not a problem, we do have the possibility here to optimize a bit by taking that into account and not do too much work in the empty regions.

What we are planning to do now:

Basic workflow remains the same:

• Transform polygon to a grid
• Calculate intersection between polygon grid and dataset grid

However, creating the polygon at the same depth of the data may be MANY cells to remember. If part of this is outside of the area with data it's quite wasteful. So as recommended in cds-astro/cds-healpix-rust#20 (comment)

• Transform polygon to MOC (efficient storage of an area in healpix)
• Transform dataset region to MOC
• Use intersection method from MOC library
• Convert back to ids.

Possibly we can think about the last step and retain the MOC structure? If the dataset is very large and high depth this can get huge as well.

[surgura]

See https://github.com/keewis/healpix-geo/pull/27