WKB/native-geometry

[espg]

References englacial/zagg#101

mortie's coverage API takes only lats/lons arrays, with no WKB/shapely entry point (verified across mortie/coverage.py).

We should:

a.) be able to ingest a WKB / WKT into our functions (coverage, etc).
b.) be able to emit something like a WKB for arrays of morton indicies; i.e., we should be able to MOC or flat encode coverage maps to binary compressed geometries.

Item (a) above I assume is pretty straightforward. What would we need to do to implement (b)? This would also presumably connect and overlap with englacial/zagg#72

[espg]

🤖 from Claude

Recommendation up front: split this into two independently-shippable tracks. (a) WKB/WKT ingest is a thin, shapely-gated optional adapter that reuses the existing lat/lon entry points — no runtime-dep change. (b) Emit is the real work, but most of its hard part already exists: mortie already exposes per-cell boundary vertices, so the missing piece is just lon/lat-degree assembly into (Multi)Polygon WKB, plus a decision on what geometry you want emitted. The central constraint that forks every option is the numpy-only runtime rule (CLAUDE.md §7; pyproject.toml:21-23 — numpy is the sole hard dep; shapely/geopandas live only in the docs extra at pyproject.toml:57-58; pandas/pyarrow are separate optional extras explicitly not runtime deps, pyproject.toml:30-45).

(a) WKB/WKT ingest

Today every coverage entry point takes parallel lat/lon arrays only: morton_coverage (mortie/coverage.py:90), morton_coverage_moc (mortie/coverage.py:195), linestring_coverage (mortie/linestring.py:38). _prep_rings (mortie/coverage.py:51) already accepts ring-lists for multipart/holes, so an adapter only has to turn a geometry into ring arrays. Options:

1. Shapely-gated optional adapter (recommended). A from_wkb / from_wkt helper that imports shapely lazily, extracts exterior+interior coords, and calls the existing functions. Lives behind the docs extra (or a new geo extra) and raises a clear "install mortie[geo]" if absent. Zero runtime-dep change, near-zero maintenance. Mirrors how arrow.py:42 already lazy-gates pyarrow.
2. Built-in numpy-only WKB/WKT reader. Parse point/linestring/polygon/multipolygon byte structure ourselves (WKB is a simple typed binary format; WKT needs a small tokenizer). No new dep, but it's real parser surface to own and test (endianness, SRID flag, Z/M dimensions, curve types we'd reject). Heavier than it looks once you cover the geometry zoo.
3. Keep I/O out of mortie entirely — document that callers pass already-parsed coordinate arrays, and provide the WKB↔arrays glue only in a notebook/example. Lowest cost, but punts the actual ask.

Tradeoff is dependency footprint vs parser-maintenance burden. (1) keeps mortie numpy-only at runtime while giving a real ingest path; (2) is the only fully-self-contained answer but adds a parser to maintain; (3) avoids both but doesn't really deliver (a).

(b) Emit — the harder half (but less missing than expected)

Emitting geometry from a morton cover needs (i) per-cell corner vertices and (ii) assembly into WKB. Finding: (i) already exists. boundaries(depth, pixel, step) (mortie/_healpix.py:50 → Rust rust_boundaries, src_rust/src/lib.rs:350, on cell_corners at src_rust/src/cell_geom.rs:75) returns a cell's corner unit-vectors (step>1 densifies each edge), and mort2healpix (mortie/tools.py:741) + infer_order_from_morton (mortie/tools.py:125) turn a (mixed-order) morton word into the (cell_id, order) that drives it. So emit is mostly: morton → (cell, order) → corners → lon/lat degrees (vec2ang, mortie/_healpix.py:77) → WKB ring bytes. The open design question is what to emit:

1. Per-cell boxes — one polygon per cell (a MultiPolygon/GeometryCollection). Trivial from the above, exact, but verbose and visually a grid of quads.
2. Dissolved boundary of the cover — a single (Multi)Polygon outline of the whole region. This is the "nice" geometry, but dissolving HEALPix-cell quads on the sphere is non-trivial (shared-edge cancellation, holes, multi-component) and is exactly the kind of polygon work shapely would normally do — which collides with numpy-only at runtime. Would likely be a shapely-gated helper, same as (a).
3. Lean on the existing dtype + Arrow/Parquet binary as the "compressed" interchange. mortie already serializes morton words to a binary, parquet-round-trippable form via the Arrow ExtensionType (mortie/arrow.py) and pandas ExtensionArray (mortie/morton_index.py); a MOC word self-encodes its order, so a mixed-order cover is already a compact uint64 array. If "binary compressed geometry" means "a small binary blob that reconstructs the region," the morton/MOC array is that, and Arrow/Parquet is the wire format — no per-cell polygon emit needed. This is the cheapest and most mortie-native answer, and is where this overlaps zagg#72 and the dtype threads (v1.0 docs: morton_index datatype user documentation #64, morton_index datatype skin: pandas + Arrow (phases 4 & 5 of #35) #51). It's a different artifact from WKB, though — geometry tools won't read it as geometry.

My read: (3) is the right default for storage/interchange; (1) is a cheap "give me real polygons" escape hatch we can ship numpy-only; (2) is the only one that needs heavy polygon machinery and should be shapely-gated like (a). This also intersects #61's mixed-order storage form — whatever canonical on-disk shape #61 lands on should be the thing (3) emits.

Risks / unknowns

• numpy-only runtime is the gating constraint for any dissolve/WKT-tokenizer path; both (a)(1) and (b)(2) want shapely, which can only be an extra.
• Antimeridian / pole wrapping: HEALPix cell quads touching the dateline or a pole need splitting/coordinate care or the WKB will self-cross. Note coverage already carries hard-won antimeridian/near-pole handling (src_rust/src/coverage.rs polar-bulge logic; mortie/tests/test_antimeridian_fix.py) — emit must not regress that.
• CRS assumptions: emit would be lon/lat degrees (EPSG:4326, vec2ang output). WKB carries no CRS unless we use the EWKB SRID flag — decide whether to set it.
• v1.0: mixed-order MOC semantics  #61 interplay: mixed-order storage form should drive (b)(3)'s output so we don't define two competing "compressed" representations.
• WKB/native-geometry #71 ↔ zagg: zagg#101 motivates ingest, zagg#72 overlaps the emit/interchange dtype — worth confirming which side owns the WKB boundary so we don't build it twice.

Questions (answerable in one pass)

1. Ingest dep policy: is shapely-as-an-extra acceptable for WKB/WKT I/O (option (a)(1)), or must ingest be a numpy-only built-in reader (a)(2)?
2. Emit shape: for (b), do you want a single dissolved-boundary (Multi)Polygon, per-cell box geometry, or both behind a flag?
3. "Compressed geometry" definition: is the existing morton/MOC Arrow/Parquet binary (b)(3) an acceptable substitute, or do you specifically need WKB bytes that geometry tools read as polygons?
4. WKB only, or WKT too, and should emitted WKB set an SRID (EWKB) for EPSG:4326 or stay plain?

---

Generated by Claude Code

[espg]

(a) WKB/WKT ingest

This is frustrating. We need to keep mortie as numpy only dependency, but I don't want us reimplementing that much code. We code scope a new subpackage surface, i.e., mortie[extras] where this could live, but not thrilled about that either. Even if we did this, I'm hesitant to use shapely here, given what a pain it is in geometric space-- much reasonable to use spherely if we where going to provide this within this library.

The lazy import has legs-- no mortie[extras], just checks if shapely/spherely exists on lazy import. I like that. What's the surface of this switching to spherely? Or checking for either spherely or shapely at run time?

My read: (3) is the right default for storage/interchange; (1) is a cheap "give me real polygons" escape hatch we can ship numpy-only; (2) is the only one that needs heavy polygon machinery and should be shapely-gated like (a). This also intersects #61 mixed-order storage form — whatever canonical on-disk shape #61 lands on should be the thing (3) emits.

Yeah, I agree with your read here.

CRS assumptions: emit would be lon/lat degrees (EPSG:4326, vec2ang output). WKB carries no CRS unless we use the EWKB SRID flag — decide whether to set it.

What's the pain of registering mortie as a projection in proj? These two already exist:

https://proj.org/en/stable/operations/projections/healpix.html
https://proj.org/en/stable/operations/projections/rhealpix.html

WKB only, or WKT too

Both WKB and WKT ideally

...and should emitted WKB set an SRID (EWKB) for EPSG:4326 or stay plain?

What's correct, best practice for this??

[espg]

🤖 from Claude

Answering the four follow-ups.

1. spherely vs shapely / "check for either at runtime"

What we actually need from a geometry lib here is narrow: a WKB/WKT codec (bytes ↔ rings) plus coordinate extraction. We do not need its spherical predicates — the spherical correctness (antimeridian/pole splitting) stays mortie's own job either way. That reframes the choice:

• shapely (≥2.0, GEOS-backed) is planar, but for pure I/O that's irrelevant — from_wkb/to_wkb/from_wkt/to_wkt and get_coordinates/.exterior.coords are rock-solid and ubiquitous. It's just a (de)serializer for us.
• spherely (s2geometry-backed) is the "right" spherical model, but it's pre-1.0; its WKB/WKT and coordinate-access surface is still stabilizing, and it pulls s2geometry.

Surface of supporting either: a thin mortie/geometry.py adapter that lazily try: import spherely then falls back to import shapely, normalizing to one internal interface (geometry → lon/lat-degree exterior+interior coord arrays; rings → WKB/WKT). The real cost is the two-backend shim: their coordinate-accessor APIs differ, so you maintain a small compat layer + tests run against whichever is installed, and CI's test extra must carry at least one. No mortie[geo] extra needed — plain lazy try/except ImportError with a clear "install shapely or spherely" message, exactly how mortie/arrow.py gates pyarrow. Runtime stays numpy-only.

My lean: make shapely≥2 the primary I/O backend (mature, the codec is all we lean on) and accept spherely if it's the one present, but only rely on the parts of its API that are stable. Since we use the lib only as a codec, "either at runtime" is reasonable — you're just paying for a 2-API shim.

2. Registering mortie as a PROJ projection

I'd advise not to — it conflates two different things. PROJ's healpix/rhealpix are continuous map projections (lon/lat ↔ projected x/y plane). mortie's morton word is a discrete cell index / DGGS, which PROJ has no model for — there's no way to express "morton cell N" as a PROJ CRS, so there's nothing to register for the index itself, and a new PROJ projection wouldn't add anything the existing healpix/rhealpix don't (the index↔cell mapping isn't a projection). Adding one is also expensive (C in PROJ, upstream contribution + release cycle). Recommendation: emit EPSG:4326 lon/lat degrees (what vec2ang already produces) and document that PROJ's healpix/rhealpix are there for anyone wanting the projected plane. A real DGGS identity is the OGC DGGS / custom-authority-code conversation — separate and much heavier.

3. WKB + WKT both

Agreed, and it's near-free: shapely/spherely give both codecs, so from_wkt/to_wkt sit right next to the WKB pair in the same adapter.

4. EWKB / SRID best practice

Plain ISO/OGC WKB carries no SRID (geometry only); EWKB (PostGIS extension) embeds one via a flagged header. Best practice in the modern geo stack:

• Default to plain WKB/WKT, no SRID. Every consumer reads it, no surprises. The CRS travels out-of-band — GeoParquet/GeoArrow store CRS in metadata, not per-geometry EWKB, and that's the dominant pattern.
• Offer EWKB/EWKT with SRID=4326 as opt-in (ewkb=True / srid= on emit) for PostGIS users who want it — but not the default. (shapely's own to_wkb(..., include_srid=False) defaults the same way.)

So: default plain, optional SRID-tagged, and document "output is lon/lat degrees, EPSG:4326" as the contract.

Net path

(a) ingest + emit both live behind one lazy mortie/geometry.py adapter (shapely-primary, spherely-accepted), numpy-only at runtime; (b) emit plain WKB/WKT by default, EWKB/EWKT(SRID=4326) opt-in; (c) no PROJ registration. If that reads right, I can turn it into a plan (then implement) — say the word and I'll write up the phased plan (ingest adapter → per-cell emit → optional dissolve behind the same gate). Refs #61 (mixed-order storage form), zagg#101/#72.

---

Generated by Claude Code

[espg]

My lean: make shapely≥2 the primary I/O backend (mature, the codec is all we lean on) and accept spherely if it's the one present, but only rely on the parts of its API that are stable. Since we use the lib only as a codec, "either at runtime" is reasonable — you're just paying for a 2-API shim.

sounds reasonable.

Registering mortie as a PROJ projection

we'll skip it

4. EWKB / SRID best practice

• Default to plain WKB/WKT, no SRID.
• Offer EWKB/EWKT with SRID=4326 as opt-in

approved for both of these.