indexing existing cell ids

[keewis]

I'm trying to select existing cell ids using a geometry (similarly to what you'd do in a tabular database), where the existing cell ids represent a subset of the full sphere (a MOC, if I understand correctly?) and are ordered but not necessarily a range.

I could use the different *_search functions like polygon_search or cone_search, but this may have a significant overhead in memory if the given vector of cell ids only covers a small part of the search area. Thus, what I'm looking for is one or multiple functions that perform the cell-in-geometry check and return a boolean / coverage mode (I can then vectorize that in a separate library for use in python).

Looking at how the polygon_search method is currently written I don't think we can make use of that?

Otherwise, am I missing anything? Do you know a better way of implementing something like this? And is this the right library to put such a function (or group of functions)?

[fxpineau]

Have you already had a look at the MOC library and its Python wrapper MOCPy?

If I understand correctly what you want to do, and assuming your "existing cells" all have the same depth, I would:

• create a MOC from your list of existing cell ids
• create a MOC from your geometry (at the depth of your existing cells)
• perform the intersection of both MOCs

Then simply iterate over the cells in the resulting MOC, it is the list of existing cells which are in the geometry.

A BMOC is a list of cells at various depth (+ a binary flag).
A MOC is also a list of cells at various depth but represented as a set of range of cells at depth 29 (the maximum depth).
Both are compact so the memory overhead should not be a problem if you use relatively small areas or not too large depths.

To have an idea of the size in memory of a MOC, you can create FITS files from MOCCli, and look at the file size (minus the 5.625 KB of FITS Header overhead). E.g.:

moc from cone --help
moc from cone 15 0.0 +0.0 0.1 fits moc.fits

creates the MOC of a cone of center (0, 0) and radius = 0.1 at order 15.
The output FITS file is 12 KB large, so the data part should be about 6 KB.

[keewis]

I had heard of MOCs before but didn't really understand them, which means I also didn't know that these may use the unique scheme to compress the cells returned by the polygon rasterization (reducing the memory footprint significantly). If that's the case I may indeed be able to use the existing functions (and the moc library) to implement what I'm looking for.

[fxpineau]

A unique can be represented as a range of depth 29 cells.
A set of non-overlapping unique can thus be represented as a a set of non-overlapping depth 29 ranges.
If the set is ordered, it will be easier to manipulate or to compare with other such sets.
E.g. it is easy to merge consecutive ranges (to reduce the memory footprint).
We have algorithms to transform back a set of depth 29 non-overlapping ranges into a set of unique.

In the Rust MOC library (I think in all MOC libraries now, which was not true in the first versions of the Java library), a MOC is a set of such ordered, non overlapping depth 29 ranges.
The Rust MOC library uses the Rust HEALPix library to get BMOC for various geometries and transforms it (i is a set of ordered zuniq), dropping the flag, into an ordered set of depth 29 non-overlapping ranges.

So no worries for the memory footprint (unless you use very high depths with very wide areas).

P.S: a zuniq is another way to code cell number + the depth in a single integer, but made so that the natural ordering of the integer follows the z-order curve order (NESTED numbers are base on the z-order curve), i.e. the position of a zuniq in the list does not depend of its depth.
The natural ordering of unique numbers places high resolution cells first.

[keewis]

that's really interesting, thanks, I didn't know there were two different encoding schemes for nested + depth (and as always, my knowledge on MOCs is very limited).

My datasets are usually low to medium resolution (depths 6 - 13), so I guess I won't need to worry about memory footprints (at least not upfront). I'll experiment and get back to you.

[keewis]

@surgura tried to use RangeMOC instances for both in EOPF-DGGS/healpix-geo#27 and noticed that for polygons that have edges that follow a meridian (and thus have exactly equal longitude values) RangeMOC::from_polygon panics. However, as soon as the values differ slightly (e.g. 1e-8 degrees) the RangeMOC is returned properly.

@fxpineau, do you know what is going on there? (Also, I'm happy to create a dedicated issue on cds-moc-rust if you prefer).

[fxpineau]

Would you provide me with an example of polygon leading to a panic (I hope it happens only on debug mode, and not in release mode; if so it may come from a to restrictive asserts, else from a real bug).

I will add it to the tests and fix the issue.

[keewis]

an example would be what you get from a bounding box (ignoring the fact that bounding boxes don't have great-circle edges):

vec![
    (154.0, -44.0),
    (154.0, -10.0),
    (113.0, -10.0),
    (113.0, -44.0),
    (154.0, -44.0),
]

and I can indeed confirm that it does not panic in release mode, only in debug.

The assert is at

cds-healpix-rust/src/special_points_finder.rs

Lines 459 to 466 in ac8a3d2

	debug_assert!(
	p1.lon() < p2.lon(),
	"p1: ({}, {}); p2: ({}, {})",
	p1.lon().to_degrees(),
	p1.lat().to_degrees(),
	p2.lon().to_degrees(),
	p2.lat().to_degrees()
	);

and the error message is:

PanicException: p1: (154, -44); p2: (154, -44)

Edit: In this case it raises at the duplicated point (which may indicate a length-0 edge depending on how you interpret the list of coordinates), but this also happens for the other edges.

[fxpineau]

Before the debug_assert, I added an early escape when two points have the same latitude.

In the library, the last point automatically connects to the first one to close the polygon, so you should remove the duplicated point in your example.

[keewis]

the functionality in healpix-geo has been on halt while I'm working on the index I mentioned in cds-astro/cds-moc-rust#12, but I'll add a preprocessor that normalizes the geometry to not contain the closing vertex once that is done.