numpy-level point-kind geo2mort: geo2mort(lat, lon, order=29, points=True)

[espg]

What's needed

A vectorized, numpy-in/numpy-out point-kind encoder on the geo2mort surface:

geo2mort(lats, lons, order=29, points=True)  # -> np.ndarray[uint64], Kind::Point words

• points=False (default): exactly today's behavior — area-kind words, byte-identical, fully backward compatible.
• points=True: route through the kernel's existing point encoding (rust_ang2pix at order 29 → rust_mi_from_nested_point, i.e. Kind::Point via from_nested_point in src_rust/src/decimal_morton.rs) and return a plain contiguous uint64 ndarray — no pandas ExtensionArray wrapper.
• Point encoding is order-29-only, so points=True with an explicit order != 29 should raise ValueError, matching MortonIndexArray.from_latlon's existing contract (mortie/morton_index.py).
• Null/sentinel semantics: same as the area path — inputs are dense float64 lat/lon with no null mask; non-finite lat/lon should fail loudly (or produce the documented sentinel 0, whichever the area path does today — the two paths must agree). 0 remains the only reserved word (empty/invalid), never a valid encode result.
• Tests to include: (a) points=True words differ from area words at order 29 but clip2order coarsening is bit-identical at every order (this is the property zagg's grouping relies on); (b) a lone points=True word survives common_ancestor unchanged (kind preserved); (c) points=True, order != 29 raises; (d) dtype is uint64, C-contiguous, same shape as input.

Why (the hot-path numbers)

zagg's HealpixGrid.assign runs on every observation of every granule read (per-photon for ATL03 — millions of points per shard). It now needs point-kind words end-to-end (englacial/zagg#87: located t-digest channel via common_ancestor). The only public point encode today is the pandas ExtensionArray classmethod MortonIndexArray.from_latlon(..., points=True), and the public ways out of the wrapper are element-wise. Measured on 1M random points (mortie 0.8.4):

path	ms / 1M pts
geo2mort (area)	42.6
raw kernel rust_ang2pix + rust_mi_from_nested_point	46.0
from_latlon(points=True) construction alone	45.2
+ public unwrap (np.asarray / .to_numpy / .astype)	~490–560

The kernel is at parity with the area encode (~8%); the ~12x penalty is entirely the ExtensionArray→ndarray unwrap. zagg currently works around it by reading the array's internal uint64 storage through its own boundary adapter (zagg.grids.morton.morton_words), which is fast (~70 ms total) but leans on a private attribute. A numpy-level geo2mort(..., points=True) gives the hot path a public vectorized surface and drops the wrapper entirely.

Consumer

englacial/zagg#87 / englacial/zagg#165 — HealpixGrid.assign swaps one call (from_latlon(points=True) + unwrap → geo2mort(..., points=True)) the moment this lands in a release; that swap is the PR's tracked phase 6.

[espg]

🤖 from Claude

Plan: geo2mort(lats, lons, order=…, points=True) — a numpy-level point-kind encoder

Recommendation: add a fused single-pass Rust point path (mirroring the area path exactly) rather than composing the existing rust_ang2pix + rust_mi_from_nested_point bindings in Python. Both land at kernel parity (the issue's ~46 ms/1M), but the fused path avoids materializing the intermediate NESTED int64 array (and its int64→uint64 recast) on the hot path and keeps geo2mort's Rust surface a single vectorized entry point. The compose-in-Python route is the lower-effort fallback (option A below) if you'd rather not add a scalar.

Where the pieces are (grounded)

• Python surface: geo2mort(lats, lons, order=18) at mortie/tools.py:106-122 — thin wrapper over _rust_geo2mort (tools.py:11), always returns np.atleast_1d(result).
• Area binding: rust_geo2mort at src_rust/src/lib.rs:206-253 — order guard, scalar/array/broadcast handling, a py.allow_threads rayon pass over geo2mort::geo2mort_scalar, returns uint64 via into_pyarray_bound.
• Area scalar: geo2mort_scalar(lat, lon, order) at src_rust/src/geo2mort.rs:24-28 = get(order).hash(Degrees(lon, lat)) → from_nested(nest, order).
• Point kernel: decimal_morton::from_nested_point(nested) at src_rust/src/decimal_morton.rs:540 (order-29-only; rejects oversized nested), vectorized as rust_mi_from_nested_point at src_rust/src/lib.rs:866-889.
• Existing point path (the one to subsume): MortonIndexArray.from_latlon(..., points=True) at mortie/morton_index.py:151-181 does rust_ang2pix(MAX_ORDER, lon, lat) → np.ascontiguousarray(..., uint64) → rust_mi_from_nested_point(...), and raises ValueError when points and int(order) != MAX_ORDER (morton_index.py:164-168). That ValueError contract is the one this issue asks geo2mort to mirror.

One ergonomic wrinkle worth a decision (see Question 1)

geo2mort's default order is 18 (tools.py:106), whereas from_latlon's default is MAX_ORDER=29 (morton_index.py:151). So a bare geo2mort(lat, lon, points=True) would today resolve order=18 and, under the issue's "order != 29 raises" rule, would raise — losing the "no extra argument needed" ergonomic the issue calls for. The plan below assumes a order=None sentinel that resolves to 18 on the area path and 29 on the point path (so points=True needs no order arg, and points=True, order=<explicit non-29> still raises). Byte-for-byte area compatibility is preserved because the resolved default is still 18.

Phased implementation

• Phase 1 — fused Rust point path. Add geo2mort_point_scalar(lat, lon) -> u64 in src_rust/src/geo2mort.rs next to geo2mort_scalar: get(MAX_ORDER).hash(Degrees(lon, lat)) → decimal_morton::from_nested_point(nest). Extend rust_geo2mort (lib.rs:206) with a points: bool flag (default false) in its #[pyo3(signature = …)]: when points, require the resolved order to be 29 (else PyValueError, matching from_latlon), and route the scalar/array/broadcast branches through geo2mort_point_scalar. Same allow_threads rayon structure, same into_pyarray_bound → contiguous uint64 return. Rust unit tests: point word ≠ area word at order 29 for a fixed lat/lon; from_nested_point agreement with the ang2pix→from_nested_point compose; order guard.
• Phase 2 — Python surface + validation. Thread points=False and the order=None sentinel through geo2mort (tools.py:106): resolve order, raise the points=True & order!=29 ValueError in Python too (fast-fail before FFI, message aligned with morton_index.py:165), keep the np.ascontiguousarray(..., float64) coercion and the always-ndarray return. Rebuild via maturin develop --release.
• Phase 3 — tests (the four the issue specifies). In mortie/tests/: (a) points=True words differ from area words at order 29, but clip2order coarsening is bit-identical between the two at every order 0–29 (the property zagg's grouping relies on); (b) a lone points=True word is unchanged through common_ancestor (kind preserved); (c) points=True with explicit order != 29 raises ValueError; (d) result is uint64, C-contiguous, same shape as input (scalar-in → 1-element ndarray, matching today). Add a cross-check that geo2mort(..., points=True) is bit-identical to the current from_latlon(points=True) unwrap so the swap in zagg is provably equivalent.
• Phase 4 — null/non-finite parity. Pin that the point path's handling of non-finite lat/lon agrees with the area path (the issue's explicit requirement — the two must not diverge). Capture whatever the area path does today (fail loudly vs. sentinel 0) as the golden and assert the point path matches; document it in the geo2mort docstring. 0 stays the only reserved word.

Files likely touched

• src_rust/src/geo2mort.rs — new geo2mort_point_scalar + unit tests.
• src_rust/src/lib.rs — rust_geo2mort gains points flag + order guard.
• mortie/tools.py — geo2mort gains points + order=None sentinel + validation + docstring.
• mortie/tests/ — the four required tests + the from_latlon equivalence check.
• No new dependency; no change to the area code path's output.

Acceptance criteria

• geo2mort(lats, lons, points=True) returns a plain contiguous uint64 ndarray of Kind::Point words, same shape as input, no ExtensionArray wrapper.
• points=False (default) is byte-identical to today at every order — full backward compatibility.
• points=True with explicit order != 29 raises ValueError; the point/area paths agree on non-finite inputs; 0 remains the sole reserved word.
• The four issue-specified tests pass, plus bit-identical equivalence to from_latlon(points=True); cargo test/cargo clippy clean and pytest -v green.

Questions for review (answerable in one pass)

1. Order-default ergonomics. Adopt the order=None sentinel (area→18, point→29, so points=True needs no order arg and points=True, order=<non-29> raises)? The alternative — keep order=18 literal — means bare geo2mort(lat, lon, points=True) raises unless the caller also passes order=29, which is more surprising. Recommend the sentinel.
2. Fused kernel vs. compose (option A). Add the geo2mort_point_scalar fused path (recommended — one rayon pass, no intermediate NESTED array, mirrors the area binding), or just compose the existing rust_ang2pix(29,…) + rust_mi_from_nested_point bindings inside the Python wrapper (no new Rust, but materializes + recasts an intermediate int64 nested array per call)? Both hit kernel parity; the fused path is the cleaner hot-path surface.
3. Non-finite semantics. Confirm the intended contract is "match whatever the area path does today" (I'll pin it as golden), rather than defining new point-specific behavior here.

Marking waiting on your read of (1)–(3). No code/branch from me on this plan issue.

---

Generated by Claude Code

[espg]

geo2mort's default order is 18 (tools.py:106)

This is a holdover from from when order 18 was the max; now that 18 isn't the max, the default for geo2mort should be updated to order 29 / point.

Order-default ergonomics. Adopt the order=None sentinel (area→18, point→29, so points=True needs no order arg and points=True, order= raises)? The alternative — keep order=18 literal — means bare geo2mort(lat, lon, points=True) raises unless the caller also passes order=29, which is more surprising. Recommend the sentinel.

see above. order=None sentinel is fine, but is belt and suspenders once the default on geo2mort is updated from order 18 (which you should do)

Fused kernel vs. compose (option A)

yes, add the recommended geo2mort_point_scalar fused path

Non-finite semantics. Confirm the intended contract is "match whatever the area path does today" (I'll pin it as golden), rather than defining new point-specific behavior here.

I believe the above answers this