From e703aafb66cd0169a201d5547480b63aca395eae Mon Sep 17 00:00:00 2001 From: claude Date: Tue, 21 Jul 2026 10:56:20 +0000 Subject: [PATCH 1/4] add HEALPix interchange guide (issue #63) --- README.md | 5 +- docs/healpix_interchange.md | 186 ++++++++++++++++++++++++++++++++++++ 2 files changed, 190 insertions(+), 1 deletion(-) create mode 100644 docs/healpix_interchange.md diff --git a/README.md b/README.md index 4720434..07145d4 100644 --- a/README.md +++ b/README.md @@ -23,7 +23,10 @@ The normative encoding and conventions — the packed-word bit layout, the decimal string grammar, the order 0–29 resolution table, the morton-hive store layout, and the coverage-MOC serializations, all frozen for the 1.x series — are documented in -[docs/specification.md](docs/specification.md). +[docs/specification.md](docs/specification.md). Moving a packed word to and from +the wider HEALPix ecosystem (`cdshealpix` / `healpy` `(order, nested-pixel)` +pairs) is covered in +[docs/healpix_interchange.md](docs/healpix_interchange.md). ## Performance diff --git a/docs/healpix_interchange.md b/docs/healpix_interchange.md new file mode 100644 index 0000000..32c8d4a --- /dev/null +++ b/docs/healpix_interchange.md @@ -0,0 +1,186 @@ +# HEALPix-ecosystem interchange + +A `mortie` packed 64-bit `morton_index` word *is* a HEALPix NESTED cell with its +order carried intrinsically (see [specification.md](specification.md) §1). This +guide shows how to move a word to and from the two common HEALPix carriers — +`(order, nested-pixel-id)` pairs as used by +[`cdshealpix`](https://cds-astro.github.io/cds-healpix-python/) and +[`healpy`](https://healpy.readthedocs.io/) — with runnable, cross-checked +snippets, then covers the three gotchas that bite at the boundary. + +Only **numpy** is a `mortie` runtime dependency. `cdshealpix` (with its +`astropy` dependency) ships in the `test` extra — it is the same oracle the +test-suite cross-checks against (`mortie/tests/test_morton_index.py`, +`test_coverage_hemisphere.py`). `healpy` is **not** in any `mortie` extra; +install it separately (`pip install healpy`) to run the healpy snippets below. + +## The conversion API + +| direction | function | returns | +|---|---|---| +| word → HEALPix | `mort2healpix(morton)` | `(nested_cell_id, order)` | +| word → normalized | `mort2norm(morton)` | `(normed, base_cell, order)` | +| normalized → word | `norm2mort(normed, base_cell, order)` | packed `uint64` word | +| order of a word | `infer_order_from_morton(morton)` | `int` | + +`mort2healpix` is the one-call bridge to the ecosystem: it returns the NESTED +`ipix` and the `order` (the HEALPix `depth`) that every other library keys on. +The nested id is the standard `base_cell * nside**2 + normed` composition with +`nside = 2**order`, so `mort2norm` / `norm2mort` are just the split/join around +that same relation. + +```python +import numpy as np +import mortie + +# geo2mort takes (lat, lon) and returns a length-1 array for scalar input, +# so index [0] for a single word. +m = mortie.geo2mort(-80.0, 120.0, order=6)[0] +print(repr(m)) # np.uint64(11570310392668225542) + +cell_id, order = mortie.mort2healpix(m) +print(cell_id, order) # 37010 6 + +normed, base_cell, order = mortie.mort2norm(m) +print(normed, base_cell, order) # 146 9 6 + +# nested id is base_cell * nside**2 + normed +nside = 2 ** order +print(cell_id == base_cell * nside**2 + normed) # True + +# round-trip back to the exact same word +m2 = mortie.norm2mort(normed, base_cell, order) +print(np.uint64(m2) == m) # True + +print(mortie.infer_order_from_morton(m)) # 6 +``` + +For an array of words `mort2healpix` returns `(cell_ids_array, order)` with a +single `order` when the words share one order (it raises on mixed orders — pass +one order at a time). + +## Round-trip against cdshealpix + +`cdshealpix.nested.lonlat_to_healpix(lon, lat, depth)` is the NESTED-`ipix` +oracle; `healpix_to_lonlat(ipix, depth)` inverts it. The cell ids match +`mort2healpix` exactly, and the centers agree to floating-point tolerance. + +```python +import numpy as np +import astropy.units as u +from cdshealpix.nested import lonlat_to_healpix, healpix_to_lonlat +import mortie + +lats = np.array([0.0, 41.8, -41.8, 80.0, -80.0, 12.3, -67.9]) +lons = np.array([0.0, 45.0, 135.0, 200.0, 305.0, 91.5, 270.2]) +order = 14 + +words = mortie.geo2mort(lats, lons, order=order) +cell_ids, o = mortie.mort2healpix(words) + +oracle = np.asarray(lonlat_to_healpix(lons * u.deg, lats * u.deg, depth=order), + dtype=np.int64) +print(np.array_equal(cell_ids, oracle)) # True +# cell_ids == [1275068416, 67108860, 2617245699, 800177021, +# 2962305922, 1556551876, 2997005193] + +# centers, mortie vs cdshealpix +clon, clat = healpix_to_lonlat(cell_ids, depth=order) +mlat, mlon = mortie.mort2geo(words) +print(float(np.max(np.abs(mlat - clat.to_value(u.deg))))) # ~1.98e-09 +print(float(np.max(np.abs(mlon - clon.to_value(u.deg))))) # 0.0 +``` + +Note the argument **order**: `geo2mort` takes `(lat, lon)`, while the cdshealpix +calls take `(lon, lat)` — a frequent source of transposed results. + +### Importing a foreign nested id + +To build a `mortie` word from a NESTED id produced elsewhere, split it into +`(normed, base_cell)` at that order and call `norm2mort`: + +```python +foreign = int(oracle[3]) # 800177021, a cdshealpix NESTED ipix at order 14 +nside_sq = (2 ** order) ** 2 +base_cell, normed = divmod(foreign, nside_sq) +word = mortie.norm2mort(normed, base_cell, order) +print(mortie.mort2healpix(word)) # (800177021, 14) +``` + +## Round-trip against healpy + +`healpy` keys on `nside = 2**order` (not the order itself) and needs +`nest=True`; pass `lonlat=True` to give it degrees in `(lon, lat)` order. + +```python +import numpy as np +import healpy as hp +import mortie + +lats = np.array([0.0, 41.8, -41.8, 80.0, -80.0, 12.3, -67.9]) +lons = np.array([0.0, 45.0, 135.0, 200.0, 305.0, 91.5, 270.2]) +order = 14 +nside = 2 ** order + +words = mortie.geo2mort(lats, lons, order=order) +cell_ids, _ = mortie.mort2healpix(words) + +hpix = hp.ang2pix(nside, lons, lats, nest=True, lonlat=True) +print(np.array_equal(cell_ids, hpix)) # True + +hlon, hlat = hp.pix2ang(nside, cell_ids, nest=True, lonlat=True) +mlat, mlon = mortie.mort2geo(words) +print(float(np.max(np.abs(mlat - hlat)))) # ~1.98e-09 +``` + +## Gotchas at the boundary + +### 1. The 4-bit prefix is `base_cell + 1`, and the word is unsigned + +The top nibble of a word stores the HEALPix base cell as `base + 1` (so the 12 +base cells occupy `1..=12`, and `0` is the empty/null sentinel). Base cells +**7–11 set bit 63**, making the word a large unsigned value. Store and exchange +the word as `uint64`; reinterpreting it as signed `int64` reads southern base +cells as negative and corrupts the id. + +```python +import numpy as np +import mortie + +w = mortie.geo2mort(-80.0, 120.0, order=3)[0] # lands in a southern base cell +print(repr(w)) # np.uint64(11565243843087433731) +print(int(w) >> 60) # 10 (prefix nibble == base_cell + 1) +_, base_cell, _ = mortie.mort2norm(w) +print(base_cell) # 9 +print(bool(int(w) >> 63 & 1)) # True (bit 63 set) +print(np.int64(w)) # -6881500230622117885 <- WRONG if read signed +``` + +The "negative = southern" signed form only exists in the decimal *string* +presentation ([specification.md](specification.md) §2); it is never a storage +form. Keep interchange columns `uint64`. + +### 2. Order is self-encoded — don't carry it out of band + +The order is intrinsic to the word (`infer_order_from_morton` reads it back), so +a word never needs a companion order column. But every ecosystem call *does* +need the order/`depth`/`nside` explicitly — always pass the `order` that +`mort2healpix` returned, and never mix orders in a single `mort2healpix` array +call (it raises on mixed orders). + +### 3. Points decode at order 29 + +A bare `geo2mort(lat, lon)` (no `order`) encodes an order-29 **point** — a +location with no area claim ([specification.md](specification.md) §1). It still +converts cleanly, but always lands at order 29: + +```python +p = mortie.geo2mort(-80.0, 120.0)[0] # order-29 point word +print(mortie.mort2healpix(p)) # (2604365600141906640, 29) +print(mortie.infer_order_from_morton(p)) # 29 +``` + +If you need an area cell at a specific resolution for interchange, pass an +explicit `order` to `geo2mort`. The empty/null sentinel (`geo2mort` of a +non-finite coordinate) is the all-zero word `0`, which is not a valid HEALPix +cell — filter it before handing ids to an ecosystem library. From a35253ce9400dc0d03c19aec90f21f7e7f4066b6 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 22 Jul 2026 10:42:58 +0000 Subject: [PATCH 2/4] fold review: version-stable prints, pinned-output note, precise base-cell wording (issue #63) --- docs/healpix_interchange.md | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/docs/healpix_interchange.md b/docs/healpix_interchange.md index 32c8d4a..0fd085c 100644 --- a/docs/healpix_interchange.md +++ b/docs/healpix_interchange.md @@ -36,7 +36,7 @@ import mortie # geo2mort takes (lat, lon) and returns a length-1 array for scalar input, # so index [0] for a single word. m = mortie.geo2mort(-80.0, 120.0, order=6)[0] -print(repr(m)) # np.uint64(11570310392668225542) +print(int(m)) # 11570310392668225542 cell_id, order = mortie.mort2healpix(m) print(cell_id, order) # 37010 6 @@ -91,6 +91,10 @@ print(float(np.max(np.abs(mlat - clat.to_value(u.deg))))) # ~1.98e-09 print(float(np.max(np.abs(mlon - clon.to_value(u.deg))))) # 0.0 ``` +The `cell_ids` list and the tolerance figures above are captured from the +current build — regenerate them by rerunning the snippet rather than editing +them by hand. + Note the argument **order**: `geo2mort` takes `(lat, lon)`, while the cdshealpix calls take `(lon, lat)` — a frequent source of transposed results. @@ -140,15 +144,15 @@ print(float(np.max(np.abs(mlat - hlat)))) # ~1.98e-09 The top nibble of a word stores the HEALPix base cell as `base + 1` (so the 12 base cells occupy `1..=12`, and `0` is the empty/null sentinel). Base cells **7–11 set bit 63**, making the word a large unsigned value. Store and exchange -the word as `uint64`; reinterpreting it as signed `int64` reads southern base -cells as negative and corrupts the id. +the word as `uint64`; reinterpreting it as signed `int64` makes base cells +7–11 read back negative and corrupts the id. ```python import numpy as np import mortie w = mortie.geo2mort(-80.0, 120.0, order=3)[0] # lands in a southern base cell -print(repr(w)) # np.uint64(11565243843087433731) +print(int(w)) # 11565243843087433731 print(int(w) >> 60) # 10 (prefix nibble == base_cell + 1) _, base_cell, _ = mortie.mort2norm(w) print(base_cell) # 9 From 0f7d602aa27a1d57a643a01759f84fffe0d25a54 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 22 Jul 2026 15:14:47 +0000 Subject: [PATCH 3/4] note healpix-rs as the backing crate (issue #63) --- docs/healpix_interchange.md | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/docs/healpix_interchange.md b/docs/healpix_interchange.md index 0fd085c..0ae38a1 100644 --- a/docs/healpix_interchange.md +++ b/docs/healpix_interchange.md @@ -8,7 +8,18 @@ guide shows how to move a word to and from the two common HEALPix carriers — [`healpy`](https://healpy.readthedocs.io/) — with runnable, cross-checked snippets, then covers the three gotchas that bite at the boundary. -Only **numpy** is a `mortie` runtime dependency. `cdshealpix` (with its +mortie does not reimplement the HEALPix geometry: its NESTED transforms +(`geo2mort`, `mort2geo`, `mort2healpix`, `norm2mort`, and the coverage kernels) +wrap the Rust [`healpix`](https://github.com/matt-cornell/healpix-rs) crate +(`matt-cornell/healpix-rs`, pinned in `Cargo.toml`) — the same crate used in +production for transforms to and from NESTED and the other wrapped HEALPix +routines. It is the only HEALPix implementation in the runtime path; there is no +C HEALPix library. `cdshealpix` and `healpy` appear below **only as independent +external oracles** — a second, unrelated implementation to cross-check that the +interchange is exact — not as runtime dependencies. + +Only **numpy** is a `mortie` runtime dependency (the `healpix` crate above is +compiled into the extension, not a Python dependency). `cdshealpix` (with its `astropy` dependency) ships in the `test` extra — it is the same oracle the test-suite cross-checks against (`mortie/tests/test_morton_index.py`, `test_coverage_hemisphere.py`). `healpy` is **not** in any `mortie` extra; From c4d598e134a6cbbaa603bbf2c21b7fba06cdd099 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 22 Jul 2026 19:30:27 +0000 Subject: [PATCH 4/4] order2res: accept array input, avoid int overflow --- mortie/tools.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/mortie/tools.py b/mortie/tools.py index f3e80a0..b3029d3 100644 --- a/mortie/tools.py +++ b/mortie/tools.py @@ -32,9 +32,16 @@ def order2res(order): equal-area), and the cell scale is the square root of that area. Derived from :data:`EARTH_RADIUS_KM` so code and the spec page (§3) share one Earth model (issue #119). + + ``order`` may be a scalar (returns a ``float``) or an array of orders such + as :func:`orders_of` yields (returns an ``ndarray``). """ - area = 4 * np.pi * EARTH_RADIUS_KM**2 / (12 * 4**order) # km2 - return float(np.sqrt(area)) + # Exponentiate in float so 4**order does not overflow an integer dtype at + # high orders (an ``orders_of`` uint8 array wraps 4**29 to 0 -> div-by-zero). + order = np.asarray(order, dtype=np.float64) + area = 4 * np.pi * EARTH_RADIUS_KM**2 / (12 * 4.0**order) # km2 + res = np.sqrt(area) + return float(res) if res.ndim == 0 else res def res2display(max_order=MAX_ORDER):