From 33c4feffb75ccb6f8a216572a9ce0d4ad3270f76 Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sun, 28 Jun 2026 10:09:37 +0000
Subject: [PATCH 1/6] phase 1 of issue #71

---
 mortie/geometry.py            | 169 ++++++++++++++++++++++++++++++++++
 mortie/tests/test_geometry.py | 117 +++++++++++++++++++++++
 pyproject.toml                |   1 +
 3 files changed, 287 insertions(+)
 create mode 100644 mortie/geometry.py
 create mode 100644 mortie/tests/test_geometry.py

diff --git a/mortie/geometry.py b/mortie/geometry.py
new file mode 100644
index 0000000..27031ba
--- /dev/null
+++ b/mortie/geometry.py
@@ -0,0 +1,169 @@
+"""Lazy WKB/WKT geometry codec for mortie (issue #71).
+
+The runtime stays **numpy-only**: this module imports a geometry backend
+(``shapely>=2`` preferred, ``spherely`` accepted) lazily and uses it *only* as a
+codec — bytes/text ↔ ring coordinate arrays.  All spherical correctness
+(antimeridian / pole handling) stays mortie's own job; the backend is never
+asked for spatial predicates.  Importing :mod:`mortie` succeeds with neither
+backend installed; the geometry functions raise a clear :class:`ImportError`
+when first touched without one (the same lazy-gate pattern :mod:`mortie.arrow`
+uses for pyarrow).
+
+Coordinate convention: WKB/WKT store ``(x, y) = (lon, lat)`` degrees
+(EPSG:4326).  mortie's coverage entry points take ``(lats, lons)``, so this
+module flips the axes at the boundary and works in degrees throughout.
+"""
+
+import numpy as np
+
+# Cached backend: a ``(name, module)`` pair, resolved once on first use.
+_BACKEND = None
+
+# GEOS / shapely geometry type ids (shapely.get_type_id); spherely follows the
+# same numbering.  Only the ones we classify on are named.
+_TYPE_POINT = 0
+_TYPE_LINESTRING = 1
+_TYPE_LINEARRING = 2
+_TYPE_POLYGON = 3
+_TYPE_MULTIPOINT = 4
+_TYPE_MULTILINESTRING = 5
+_TYPE_MULTIPOLYGON = 6
+_TYPE_GEOMETRYCOLLECTION = 7
+
+_POLYGONAL = (_TYPE_POLYGON, _TYPE_MULTIPOLYGON)
+_LINEAR = (_TYPE_LINESTRING, _TYPE_MULTILINESTRING, _TYPE_LINEARRING)
+
+
+def _require_backend():
+    """Import a geometry backend lazily, raising a clear error if absent.
+
+    ``shapely>=2`` is the primary backend (its WKB/WKT codec is mature and is
+    all we lean on); ``spherely`` is accepted if it is the one present.  Returns
+    a ``(name, module)`` pair and caches it.
+    """
+    global _BACKEND
+    if _BACKEND is not None:
+        return _BACKEND
+    try:
+        import shapely
+
+        _BACKEND = ("shapely", shapely)
+        return _BACKEND
+    except ImportError:
+        pass
+    try:
+        import spherely
+
+        _BACKEND = ("spherely", spherely)
+        return _BACKEND
+    except ImportError:
+        pass
+    raise ImportError(
+        "mortie's WKB/WKT geometry I/O requires a geometry backend; install "
+        "`shapely>=2` (preferred) or `spherely` (e.g. `pip install shapely`). "
+        "mortie's runtime is numpy-only, so the backend is an optional extra."
+    )
+
+
+def _strip_ewkt_srid(text):
+    """Drop a leading ``SRID=<n>;`` prefix from an EWKT string, if present.
+
+    Plain WKT parsers reject the PostGIS EWKT prefix, so ingest tolerates it by
+    stripping it (the SRID is advisory; mortie's contract is always EPSG:4326).
+    """
+    s = text.lstrip()
+    if s[:5].upper() == "SRID=":
+        semi = s.find(";")
+        if semi != -1:
+            return s[semi + 1:]
+    return text
+
+
+def geometry_from_wkb(data):
+    """Decode WKB (or EWKB) bytes into a backend geometry object."""
+    _, mod = _require_backend()
+    return mod.from_wkb(data)
+
+
+def geometry_from_wkt(text):
+    """Decode WKT (or EWKT) text into a backend geometry object."""
+    _, mod = _require_backend()
+    return mod.from_wkt(_strip_ewkt_srid(text))
+
+
+def geometry_to_wkb(geom, srid=None):
+    """Encode a backend geometry to WKB bytes.
+
+    With ``srid`` set (e.g. ``4326``), emit **EWKB** carrying that SRID;
+    otherwise emit plain ISO/OGC WKB (the default, no embedded CRS).
+    """
+    _, mod = _require_backend()
+    if srid is not None:
+        geom = mod.set_srid(geom, int(srid))
+        return mod.to_wkb(geom, include_srid=True)
+    return mod.to_wkb(geom)
+
+
+def geometry_to_wkt(geom, srid=None):
+    """Encode a backend geometry to WKT text.
+
+    With ``srid`` set, emit **EWKT** (``SRID=<n>;<WKT>``); otherwise plain WKT.
+    """
+    _, mod = _require_backend()
+    text = mod.to_wkt(geom)
+    if srid is not None:
+        return f"SRID={int(srid)};{text}"
+    return text
+
+
+def _ring_latlon(mod, ring_geom):
+    """Extract a ring's vertices as ``(lat, lon)`` float64 arrays (degrees)."""
+    coords = np.asarray(mod.get_coordinates(ring_geom), dtype=np.float64)
+    # WKB/WKT store (x, y) = (lon, lat).
+    return coords[:, 1].copy(), coords[:, 0].copy()
+
+
+def _polygon_rings(mod, poly):
+    """All rings of one polygon as ``(lat, lon)`` pairs: exterior then holes."""
+    rings = [_ring_latlon(mod, mod.get_exterior_ring(poly))]
+    for i in range(int(mod.get_num_interior_rings(poly))):
+        rings.append(_ring_latlon(mod, mod.get_interior_ring(poly, i)))
+    return rings
+
+
+def decompose(geom):
+    """Decompose a backend geometry into mortie coverage inputs.
+
+    Returns ``(kind, parts)`` where:
+
+    * ``kind == "polygonal"`` and ``parts`` is a list of rings — exterior and
+      interior (hole) rings of every polygon, flattened.  mortie's even-odd
+      descent covers them in one pass, so disjoint outers union and nested
+      rings carve holes (matching :func:`mortie.morton_coverage`'s multipart
+      contract).
+    * ``kind == "linear"`` and ``parts`` is a list of lines, one per
+      (multi)linestring component.
+
+    Each entry is a ``(lat, lon)`` pair of float64 degree arrays.  Points,
+    geometry collections, and empty geometries are rejected — coverage has no
+    meaning for them.
+    """
+    _, mod = _require_backend()
+    type_id = int(mod.get_type_id(geom))
+
+    if type_id == _TYPE_POLYGON:
+        return "polygonal", _polygon_rings(mod, geom)
+    if type_id == _TYPE_MULTIPOLYGON:
+        rings = []
+        for poly in mod.get_parts(geom):
+            rings.extend(_polygon_rings(mod, poly))
+        return "polygonal", rings
+    if type_id in (_TYPE_LINESTRING, _TYPE_LINEARRING):
+        return "linear", [_ring_latlon(mod, geom)]
+    if type_id == _TYPE_MULTILINESTRING:
+        return "linear", [_ring_latlon(mod, ln) for ln in mod.get_parts(geom)]
+
+    raise ValueError(
+        f"unsupported geometry type for coverage (type id {type_id}); "
+        "expected Polygon, MultiPolygon, LineString, or MultiLineString"
+    )
diff --git a/mortie/tests/test_geometry.py b/mortie/tests/test_geometry.py
new file mode 100644
index 0000000..5800e5c
--- /dev/null
+++ b/mortie/tests/test_geometry.py
@@ -0,0 +1,117 @@
+"""Tests for the WKB/WKT geometry codec adapter (issue #71, phase 1).
+
+These pin :mod:`mortie.geometry`: the lazy backend gate, the WKB/WKT (and
+EWKB/EWKT) codec, and the decomposition of polygons / multipolygons / holes /
+linestrings into ``(lat, lon)`` ring arrays — the input shape the ingest path
+(phase 2) feeds to the existing coverage entry points.  The backend is used
+only as a codec; spherical correctness is mortie's own job and not exercised
+here.
+"""
+
+import numpy as np
+import pytest
+
+from mortie import geometry
+
+shapely = pytest.importorskip("shapely")
+
+
+def test_decompose_polygon_with_hole():
+    # A unit square with a square hole: exterior + one interior ring.
+    g = shapely.from_wkt(
+        "POLYGON ((0 0, 2 0, 2 2, 0 2, 0 0),"
+        "(0.5 0.5, 1.5 0.5, 1.5 1.5, 0.5 1.5, 0.5 0.5))"
+    )
+    kind, rings = geometry.decompose(g)
+    assert kind == "polygonal"
+    assert len(rings) == 2  # exterior + hole
+    ext_lat, ext_lon = rings[0]
+    # (x, y) = (lon, lat): the exterior spans lon/lat 0..2.
+    assert np.isclose(ext_lon.max(), 2.0) and np.isclose(ext_lat.max(), 2.0)
+    hole_lat, hole_lon = rings[1]
+    assert np.isclose(hole_lon.min(), 0.5) and np.isclose(hole_lat.min(), 0.5)
+
+
+def test_decompose_multipolygon_flattens_all_rings():
+    g = shapely.from_wkt(
+        "MULTIPOLYGON (((0 0, 1 0, 1 1, 0 1, 0 0)),"
+        "((5 5, 6 5, 6 6, 5 6, 5 5),(5.2 5.2, 5.8 5.2, 5.8 5.8, 5.2 5.8, 5.2 5.2)))"
+    )
+    kind, rings = geometry.decompose(g)
+    assert kind == "polygonal"
+    # poly1 (1 ring) + poly2 (exterior + 1 hole) = 3 rings, flattened.
+    assert len(rings) == 3
+
+
+def test_decompose_linestring_and_multilinestring():
+    ls = shapely.from_wkt("LINESTRING (0 0, 1 1, 2 0)")
+    kind, lines = geometry.decompose(ls)
+    assert kind == "linear"
+    assert len(lines) == 1
+    lat, lon = lines[0]
+    assert lat.shape == (3,) and lon.shape == (3,)
+
+    mls = shapely.from_wkt("MULTILINESTRING ((0 0, 1 1), (2 2, 3 3, 4 4))")
+    kind, lines = geometry.decompose(mls)
+    assert kind == "linear"
+    assert [ln[0].size for ln in lines] == [2, 3]
+
+
+def test_decompose_rejects_points_and_collections():
+    with pytest.raises(ValueError, match="unsupported geometry type"):
+        geometry.decompose(shapely.from_wkt("POINT (1 2)"))
+    with pytest.raises(ValueError, match="unsupported geometry type"):
+        geometry.decompose(
+            shapely.from_wkt("GEOMETRYCOLLECTION (POINT (1 2), LINESTRING (0 0, 1 1))")
+        )
+
+
+def test_wkb_wkt_codec_roundtrip():
+    wkt = "POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))"
+    g = geometry.geometry_from_wkt(wkt)
+    # WKB round-trip preserves the rings.
+    wkb = geometry.geometry_to_wkb(g)
+    assert isinstance(wkb, (bytes, bytearray))
+    g2 = geometry.geometry_from_wkb(wkb)
+    k1, r1 = geometry.decompose(g)
+    k2, r2 = geometry.decompose(g2)
+    assert k1 == k2
+    assert np.allclose(r1[0][0], r2[0][0]) and np.allclose(r1[0][1], r2[0][1])
+    # WKT round-trip.
+    g3 = geometry.geometry_from_wkt(geometry.geometry_to_wkt(g))
+    assert int(shapely.get_type_id(g3)) == int(shapely.get_type_id(g))
+
+
+def test_ewkb_ewkt_srid_optin():
+    g = geometry.geometry_from_wkt("POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))")
+
+    # EWKT carries the SRID prefix; plain WKT does not.
+    ewkt = geometry.geometry_to_wkt(g, srid=4326)
+    assert ewkt.startswith("SRID=4326;")
+    assert not geometry.geometry_to_wkt(g).startswith("SRID=")
+    # The EWKT prefix is tolerated on ingest (advisory; contract is EPSG:4326).
+    g_back = geometry.geometry_from_wkt(ewkt)
+    assert int(shapely.get_type_id(g_back)) == int(shapely.get_type_id(g))
+
+    # EWKB carries the SRID; from_wkb reads it back.
+    ewkb = geometry.geometry_to_wkb(g, srid=4326)
+    assert int(shapely.get_srid(geometry.geometry_from_wkb(ewkb))) == 4326
+    assert int(shapely.get_srid(geometry.geometry_from_wkb(geometry.geometry_to_wkb(g)))) == 0
+
+
+def test_backend_gate_message(monkeypatch):
+    # With no backend importable, a clear ImportError naming shapely/spherely.
+    import mortie.geometry as gm
+
+    monkeypatch.setattr(gm, "_BACKEND", None)
+    real_import = __import__
+
+    def _block(name, *args, **kwargs):
+        if name in ("shapely", "spherely"):
+            raise ImportError(f"blocked {name}")
+        return real_import(name, *args, **kwargs)
+
+    monkeypatch.setattr("builtins.__import__", _block)
+    with pytest.raises(ImportError, match="shapely"):
+        gm._require_backend()
+    # monkeypatch reverts _BACKEND and __import__; the next call re-resolves.
diff --git a/pyproject.toml b/pyproject.toml
index 51c2033..7773af1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -42,6 +42,7 @@ test = [
   "cdshealpix>=0.7",  # independent HEALPix oracle for the hemisphere/#11 tests
   "pandas>=2.0",      # exercise the morton_index ExtensionArray in CI
   "pyarrow>=14.0",    # exercise the morton_index Arrow ExtensionType in CI
+  "shapely>=2.0",     # WKB/WKT geometry codec for the issue #71 I/O (test-only)
 ]
 bench = [
   "pytest",

From 483a9df018e44c6865b07b09f02d864fca7caaff Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sun, 28 Jun 2026 10:16:06 +0000
Subject: [PATCH 2/6] phase 2 of issue #71

---
 mortie/__init__.py            |   9 +++
 mortie/geometry.py            | 139 ++++++++++++++++++++++++++++++----
 mortie/tests/test_geometry.py | 106 ++++++++++++++++++++++++++
 3 files changed, 241 insertions(+), 13 deletions(-)

diff --git a/mortie/__init__.py b/mortie/__init__.py
index f980e00..3e8410e 100644
--- a/mortie/__init__.py
+++ b/mortie/__init__.py
@@ -26,6 +26,11 @@
     morton_coverage_moc,
     split_base_cells,
 )
+from .geometry import (
+    from_geometry,
+    from_wkb,
+    from_wkt,
+)
 from .linestring import linestring_coverage
 
 # Import prefix trie functions
@@ -95,6 +100,10 @@
     'moc_min',
     'split_base_cells',
     'linestring_coverage',
+    'from_wkb',
+    'from_wkt',
+    'from_geometry',
+    'geometry',
     'MortonChild',
     'split_children',
     'split_children_geo',
diff --git a/mortie/geometry.py b/mortie/geometry.py
index 27031ba..3e7c7f7 100644
--- a/mortie/geometry.py
+++ b/mortie/geometry.py
@@ -21,17 +21,11 @@
 
 # GEOS / shapely geometry type ids (shapely.get_type_id); spherely follows the
 # same numbering.  Only the ones we classify on are named.
-_TYPE_POINT = 0
 _TYPE_LINESTRING = 1
 _TYPE_LINEARRING = 2
 _TYPE_POLYGON = 3
-_TYPE_MULTIPOINT = 4
 _TYPE_MULTILINESTRING = 5
 _TYPE_MULTIPOLYGON = 6
-_TYPE_GEOMETRYCOLLECTION = 7
-
-_POLYGONAL = (_TYPE_POLYGON, _TYPE_MULTIPOLYGON)
-_LINEAR = (_TYPE_LINESTRING, _TYPE_MULTILINESTRING, _TYPE_LINEARRING)
 
 
 def _require_backend():
@@ -65,6 +59,27 @@ def _require_backend():
     )
 
 
+def _require_shapely(what):
+    """Require the shapely backend for *what*, raising a clear error otherwise.
+
+    The raw WKB/WKT codec works on either backend, but ring decomposition and
+    SRID-tagged emit lean on shapely's geometry-introspection API
+    (``get_exterior_ring`` / ``get_parts`` / ``set_srid``), which spherely's
+    published surface does not yet expose.  Rather than fail with an opaque
+    ``AttributeError`` deep inside, refuse up front with guidance.  Whether to
+    invest in a spherely introspection shim is an open question for the issue
+    thread (see the PR's "Questions for review").
+    """
+    name, mod = _require_backend()
+    if name != "shapely":
+        raise NotImplementedError(
+            f"{what} currently requires the shapely>=2 backend; the active "
+            f"backend is {name!r}, which mortie uses only as a raw WKB/WKT "
+            "codec. Install shapely>=2 for this operation."
+        )
+    return mod
+
+
 def _strip_ewkt_srid(text):
     """Drop a leading ``SRID=<n>;`` prefix from an EWKT string, if present.
 
@@ -94,13 +109,15 @@ def geometry_from_wkt(text):
 def geometry_to_wkb(geom, srid=None):
     """Encode a backend geometry to WKB bytes.
 
-    With ``srid`` set (e.g. ``4326``), emit **EWKB** carrying that SRID;
-    otherwise emit plain ISO/OGC WKB (the default, no embedded CRS).
+    With ``srid`` set (e.g. ``4326``), emit **EWKB** carrying that SRID
+    (shapely backend only); otherwise emit plain ISO/OGC WKB (the default, no
+    embedded CRS) — works on either backend.
     """
-    _, mod = _require_backend()
     if srid is not None:
+        mod = _require_shapely("EWKB emit (srid=)")
         geom = mod.set_srid(geom, int(srid))
         return mod.to_wkb(geom, include_srid=True)
+    _, mod = _require_backend()
     return mod.to_wkb(geom)
 
 
@@ -144,11 +161,17 @@ def decompose(geom):
     * ``kind == "linear"`` and ``parts`` is a list of lines, one per
       (multi)linestring component.
 
-    Each entry is a ``(lat, lon)`` pair of float64 degree arrays.  Points,
-    geometry collections, and empty geometries are rejected — coverage has no
-    meaning for them.
+    Each entry is a ``(lat, lon)`` pair of float64 degree arrays.  Any Z
+    coordinate is dropped (mortie is 2-D lon/lat).  Points, geometry
+    collections, and empty geometries are rejected — coverage has no meaning
+    for them.
+
+    Requires the shapely backend (it leans on shapely's ring/parts
+    introspection); see :func:`_require_shapely`.
     """
-    _, mod = _require_backend()
+    mod = _require_shapely("geometry decomposition")
+    if bool(mod.is_empty(geom)):
+        raise ValueError("empty geometry has no coverage")
     type_id = int(mod.get_type_id(geom))
 
     if type_id == _TYPE_POLYGON:
@@ -167,3 +190,93 @@ def decompose(geom):
         f"unsupported geometry type for coverage (type id {type_id}); "
         "expected Polygon, MultiPolygon, LineString, or MultiLineString"
     )
+
+
+# ── ingest: geometry → morton coverage ─────────────────────────────────────
+
+
+def from_geometry(geom, order=18, moc=False, normalize=True,
+                  tolerance=None, max_cells=None):
+    """Cover a backend geometry with morton indices (issue #71).
+
+    The geometry is decomposed via :func:`decompose` and routed to mortie's
+    existing coverage entry points — so WKB/WKT ingest produces exactly the same
+    cover as calling those functions on the same ``(lats, lons)`` arrays.
+
+    * **Polygon / MultiPolygon** → :func:`mortie.morton_coverage` (flat) or, with
+      ``moc=True``, :func:`mortie.morton_coverage_moc` (compact mixed-order).
+      Holes and disjoint parts are handled by the one even-odd descent.
+    * **LineString / MultiLineString** → :func:`mortie.linestring_coverage`.
+
+    Parameters
+    ----------
+    geom : backend geometry
+        A shapely/spherely geometry object (e.g. from :func:`geometry_from_wkb`).
+    order : int, optional
+        HEALPix order (1–29).  Default 18.
+    moc : bool, optional
+        Polygonal only: return a compact MOC instead of a flat cover.
+    normalize : bool, optional
+        Flat polygon cover only: auto-correct ring orientation at ingest
+        (see :func:`mortie.morton_coverage`).  Ignored when ``moc=True``.
+    tolerance, max_cells : optional
+        Polygonal ``moc=True`` only: the adaptive stop criteria of
+        :func:`mortie.morton_coverage_moc` (mutually exclusive).
+
+    Returns
+    -------
+    numpy.ndarray or list of numpy.ndarray
+        Polygonal → 1-D ``uint64`` morton array.  LineString → 1-D array;
+        MultiLineString → list of arrays, one per line (the
+        :func:`mortie.linestring_coverage` contract).
+    """
+    from .coverage import morton_coverage, morton_coverage_moc
+    from .linestring import linestring_coverage
+
+    kind, parts = decompose(geom)
+
+    if kind == "polygonal":
+        lats = [p[0] for p in parts]
+        lons = [p[1] for p in parts]
+        if moc:
+            return morton_coverage_moc(
+                lats, lons, order=order, tolerance=tolerance, max_cells=max_cells
+            )
+        return morton_coverage(lats, lons, order=order, normalize=normalize)
+
+    # linear
+    if moc or tolerance is not None or max_cells is not None:
+        raise ValueError(
+            "moc / tolerance / max_cells apply only to polygonal geometry"
+        )
+    if len(parts) == 1:
+        return linestring_coverage(parts[0][0], parts[0][1], order=order)
+    lats = [p[0] for p in parts]
+    lons = [p[1] for p in parts]
+    return linestring_coverage(lats, lons, order=order)
+
+
+def from_wkb(data, order=18, moc=False, normalize=True,
+             tolerance=None, max_cells=None):
+    """Cover a geometry given as WKB (or EWKB) bytes.
+
+    Thin wrapper: decode with :func:`geometry_from_wkb`, then
+    :func:`from_geometry`.  See :func:`from_geometry` for the parameters.
+    """
+    return from_geometry(
+        geometry_from_wkb(data), order=order, moc=moc, normalize=normalize,
+        tolerance=tolerance, max_cells=max_cells,
+    )
+
+
+def from_wkt(text, order=18, moc=False, normalize=True,
+             tolerance=None, max_cells=None):
+    """Cover a geometry given as WKT (or EWKT) text.
+
+    Thin wrapper: decode with :func:`geometry_from_wkt`, then
+    :func:`from_geometry`.  See :func:`from_geometry` for the parameters.
+    """
+    return from_geometry(
+        geometry_from_wkt(text), order=order, moc=moc, normalize=normalize,
+        tolerance=tolerance, max_cells=max_cells,
+    )
diff --git a/mortie/tests/test_geometry.py b/mortie/tests/test_geometry.py
index 5800e5c..78b15b1 100644
--- a/mortie/tests/test_geometry.py
+++ b/mortie/tests/test_geometry.py
@@ -11,6 +11,7 @@
 import numpy as np
 import pytest
 
+import mortie
 from mortie import geometry
 
 shapely = pytest.importorskip("shapely")
@@ -66,6 +67,21 @@ def test_decompose_rejects_points_and_collections():
         )
 
 
+def test_decompose_rejects_empty_geometry():
+    for wkt in ("POLYGON EMPTY", "LINESTRING EMPTY", "MULTIPOLYGON EMPTY"):
+        with pytest.raises(ValueError, match="empty geometry"):
+            geometry.decompose(shapely.from_wkt(wkt))
+
+
+def test_decompose_drops_z_coordinate():
+    # A 3-D polygon ingests as its 2-D lon/lat footprint (Z is dropped).
+    g = shapely.from_wkt("POLYGON Z ((0 0 5, 1 0 5, 1 1 5, 0 1 5, 0 0 5))")
+    kind, rings = geometry.decompose(g)
+    assert kind == "polygonal"
+    lat, lon = rings[0]
+    assert lat.ndim == 1 and lon.ndim == 1  # no third column leaked through
+
+
 def test_wkb_wkt_codec_roundtrip():
     wkt = "POLYGON ((0 0, 1 0, 1 1, 0 1, 0 0))"
     g = geometry.geometry_from_wkt(wkt)
@@ -99,6 +115,96 @@ def test_ewkb_ewkt_srid_optin():
     assert int(shapely.get_srid(geometry.geometry_from_wkb(geometry.geometry_to_wkb(g)))) == 0
 
 
+# ── Phase 2: ingest reproduces the array-path coverage ─────────────────────
+
+# A small polygon well away from the poles / antimeridian.
+_LATS = [40.0, 50.0, 50.0, 40.0]
+_LONS = [-120.0, -120.0, -110.0, -110.0]
+
+
+def _poly_wkt(lats, lons):
+    pts = ", ".join(f"{lo} {la}" for la, lo in zip(lats, lons))
+    first = f"{lons[0]} {lats[0]}"
+    return f"POLYGON (({pts}, {first}))"
+
+
+def test_ingest_polygon_matches_array_path():
+    want = mortie.morton_coverage(_LATS, _LONS, order=6)
+    wkt = _poly_wkt(_LATS, _LONS)
+    got_wkt = mortie.from_wkt(wkt, order=6)
+    got_wkb = mortie.from_wkb(geometry.geometry_to_wkb(shapely.from_wkt(wkt)), order=6)
+    assert np.array_equal(got_wkt, want)
+    assert np.array_equal(got_wkb, want)
+
+
+def test_ingest_polygon_with_hole_matches_array_path():
+    outer_lat, outer_lon = _LATS, _LONS
+    hole_lat = [43.0, 47.0, 47.0, 43.0]
+    hole_lon = [-117.0, -117.0, -113.0, -113.0]
+    want = mortie.morton_coverage(
+        [outer_lat, hole_lat], [outer_lon, hole_lon], order=6
+    )
+    wkt = (
+        "POLYGON (("
+        + ", ".join(f"{lo} {la}" for la, lo in zip(outer_lat, outer_lon))
+        + f", {outer_lon[0]} {outer_lat[0]}),("
+        + ", ".join(f"{lo} {la}" for la, lo in zip(hole_lat, hole_lon))
+        + f", {hole_lon[0]} {hole_lat[0]}))"
+    )
+    assert np.array_equal(mortie.from_wkt(wkt, order=6), want)
+
+
+def test_ingest_multipolygon_matches_array_path():
+    lats2 = [10.0, 20.0, 20.0, 10.0]
+    lons2 = [-80.0, -80.0, -70.0, -70.0]
+    want = mortie.morton_coverage([_LATS, lats2], [_LONS, lons2], order=6)
+    wkt = (
+        "MULTIPOLYGON ((("
+        + ", ".join(f"{lo} {la}" for la, lo in zip(_LATS, _LONS))
+        + f", {_LONS[0]} {_LATS[0]})),(("
+        + ", ".join(f"{lo} {la}" for la, lo in zip(lats2, lons2))
+        + f", {lons2[0]} {lats2[0]})))"
+    )
+    assert np.array_equal(mortie.from_wkt(wkt, order=6), want)
+
+
+def test_ingest_polygon_moc_matches_array_path():
+    want = mortie.morton_coverage_moc(_LATS, _LONS, order=8)
+    wkt = _poly_wkt(_LATS, _LONS)
+    assert np.array_equal(mortie.from_wkt(wkt, order=8, moc=True), want)
+
+
+def test_ingest_linestring_matches_array_path():
+    lats = [40.0, 50.0, 45.0]
+    lons = [-120.0, -110.0, -100.0]
+    want = mortie.linestring_coverage(lats, lons, order=6)
+    wkt = "LINESTRING (" + ", ".join(f"{lo} {la}" for la, lo in zip(lats, lons)) + ")"
+    got = mortie.from_wkt(wkt, order=6)
+    assert np.array_equal(got, want)
+
+
+def test_ingest_multilinestring_matches_array_path():
+    lats = [[40.0, 50.0], [10.0, 20.0, 15.0]]
+    lons = [[-120.0, -110.0], [-80.0, -70.0, -60.0]]
+    want = mortie.linestring_coverage(lats, lons, order=6)
+    wkt = (
+        "MULTILINESTRING (("
+        + ", ".join(f"{lo} {la}" for la, lo in zip(lats[0], lons[0]))
+        + "),("
+        + ", ".join(f"{lo} {la}" for la, lo in zip(lats[1], lons[1]))
+        + "))"
+    )
+    got = mortie.from_wkt(wkt, order=6)
+    assert isinstance(got, list) and len(got) == 2
+    assert all(np.array_equal(g, w) for g, w in zip(got, want))
+
+
+def test_ingest_linear_rejects_polygon_only_args():
+    wkt = "LINESTRING (0 0, 1 1, 2 0)"
+    with pytest.raises(ValueError, match="only to polygonal"):
+        mortie.from_wkt(wkt, order=6, moc=True)
+
+
 def test_backend_gate_message(monkeypatch):
     # With no backend importable, a clear ImportError naming shapely/spherely.
     import mortie.geometry as gm

From 373efe6096b582fc412e8a232fe4cae44727545c Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sun, 28 Jun 2026 10:21:06 +0000
Subject: [PATCH 3/6] phase 3 of issue #71

---
 mortie/__init__.py            |  6 +++
 mortie/geometry.py            | 86 ++++++++++++++++++++++++++++++++++-
 mortie/tests/test_geometry.py | 61 +++++++++++++++++++++++++
 3 files changed, 152 insertions(+), 1 deletion(-)

diff --git a/mortie/__init__.py b/mortie/__init__.py
index 3e8410e..aa6bfc1 100644
--- a/mortie/__init__.py
+++ b/mortie/__init__.py
@@ -30,6 +30,9 @@
     from_geometry,
     from_wkb,
     from_wkt,
+    to_geometry,
+    to_wkb,
+    to_wkt,
 )
 from .linestring import linestring_coverage
 
@@ -103,6 +106,9 @@
     'from_wkb',
     'from_wkt',
     'from_geometry',
+    'to_wkb',
+    'to_wkt',
+    'to_geometry',
     'geometry',
     'MortonChild',
     'split_children',
diff --git a/mortie/geometry.py b/mortie/geometry.py
index 3e7c7f7..8c834b2 100644
--- a/mortie/geometry.py
+++ b/mortie/geometry.py
@@ -218,7 +218,10 @@ def from_geometry(geom, order=18, moc=False, normalize=True,
         Polygonal only: return a compact MOC instead of a flat cover.
     normalize : bool, optional
         Flat polygon cover only: auto-correct ring orientation at ingest
-        (see :func:`mortie.morton_coverage`).  Ignored when ``moc=True``.
+        (see :func:`mortie.morton_coverage`).  Ignored when ``moc=True`` and for
+        linear geometry.  Note ``morton_coverage_moc`` has no orientation
+        auto-correct, so with ``moc=True`` the ring winding is taken **as
+        authored** — for hemisphere-plus polygons wind exteriors CCW / holes CW.
     tolerance, max_cells : optional
         Polygonal ``moc=True`` only: the adaptive stop criteria of
         :func:`mortie.morton_coverage_moc` (mutually exclusive).
@@ -280,3 +283,84 @@ def from_wkt(text, order=18, moc=False, normalize=True,
         geometry_from_wkt(text), order=order, moc=moc, normalize=normalize,
         tolerance=tolerance, max_cells=max_cells,
     )
+
+
+# ── emit: morton coverage → geometry ───────────────────────────────────────
+
+
+def _per_cell_polygons(mod, morton, step):
+    """Build one backend Polygon per cell of *morton* (lon/lat degrees).
+
+    Reuses :func:`mortie.mort2polygon` for the corner→lon/lat boundary (with its
+    antimeridian normalization), grouping by order so a mixed-order MOC cover is
+    handled — ``mort2polygon`` itself requires a single order per call.
+    """
+    from .tools import _rust_mort2nested, mort2polygon
+
+    morton = np.atleast_1d(np.asarray(morton, dtype=np.uint64))
+    if morton.size == 0:
+        return []
+
+    _, depths = _rust_mort2nested(np.ascontiguousarray(morton))
+    polys = []
+    for d in np.unique(depths):
+        grp = morton[depths == d]
+        if grp.size == 1:
+            rings_ll = [mort2polygon(int(grp[0]), step=step)]
+        else:
+            rings_ll = mort2polygon(grp, step=step)
+        for ring in rings_ll:
+            # mort2polygon yields closed [lat, lon] pairs; WKB wants (lon, lat).
+            polys.append(mod.Polygon([(lon, lat) for lat, lon in ring]))
+    return polys
+
+
+def to_geometry(morton, dissolve=True, step=1):
+    """Convert a morton cover to a backend geometry (issue #71).
+
+    Parameters
+    ----------
+    morton : array_like of uint64
+        A morton cover (flat or mixed-order MOC; each word self-encodes order).
+    dissolve : bool, optional
+        ``True`` (default) emits the single dissolved outline of the whole cover
+        (exterior rings, holes, and disjoint components).  ``False`` emits a
+        per-cell ``MultiPolygon`` — one quad per cell.
+    step : int, optional
+        Boundary points per cell edge (default 1 = 4 corners / straight chords).
+        ``step>1`` densifies each edge to follow the curved HEALPix boundary.
+
+    Returns
+    -------
+    backend geometry
+        A shapely (or spherely) ``MultiPolygon`` in EPSG:4326 lon/lat degrees.
+
+    Notes
+    -----
+    Emit requires the shapely backend (it constructs geometry objects).
+    """
+    mod = _require_shapely("geometry emit")
+    if dissolve:
+        raise NotImplementedError(
+            "dissolve=True (the dissolved-outline emit) lands in phase 4 of "
+            "issue #71; pass dissolve=False for the per-cell MultiPolygon"
+        )
+    return mod.MultiPolygon(_per_cell_polygons(mod, morton, step))
+
+
+def to_wkb(morton, dissolve=True, step=1, srid=None):
+    """Emit a morton cover as WKB (or EWKB) bytes.
+
+    See :func:`to_geometry` for ``dissolve`` / ``step``.  With ``srid`` set
+    (e.g. ``4326``), emit EWKB carrying that SRID; otherwise plain WKB.
+    """
+    return geometry_to_wkb(to_geometry(morton, dissolve=dissolve, step=step), srid=srid)
+
+
+def to_wkt(morton, dissolve=True, step=1, srid=None):
+    """Emit a morton cover as WKT (or EWKT) text.
+
+    See :func:`to_geometry` for ``dissolve`` / ``step``.  With ``srid`` set,
+    emit EWKT (``SRID=<n>;<WKT>``); otherwise plain WKT.
+    """
+    return geometry_to_wkt(to_geometry(morton, dissolve=dissolve, step=step), srid=srid)
diff --git a/mortie/tests/test_geometry.py b/mortie/tests/test_geometry.py
index 78b15b1..30fe74d 100644
--- a/mortie/tests/test_geometry.py
+++ b/mortie/tests/test_geometry.py
@@ -205,6 +205,67 @@ def test_ingest_linear_rejects_polygon_only_args():
         mortie.from_wkt(wkt, order=6, moc=True)
 
 
+def test_ingest_moc_via_wkb_and_clockwise_spelling():
+    # moc ingest works through WKB (not just WKT)...
+    want = mortie.morton_coverage_moc(_LATS, _LONS, order=8)
+    wkb = geometry.geometry_to_wkb(shapely.from_wkt(_poly_wkt(_LATS, _LONS)))
+    assert np.array_equal(mortie.from_wkb(wkb, order=8, moc=True), want)
+    # ...and a clockwise ring gives the same sub-hemisphere cover as CCW
+    # (normalize=True default makes ordinary polygons orientation-insensitive).
+    ccw = _poly_wkt(list(reversed(_LATS)), list(reversed(_LONS)))
+    cw = _poly_wkt(_LATS, _LONS)
+    assert np.array_equal(
+        mortie.from_wkt(cw, order=6), mortie.from_wkt(ccw, order=6)
+    )
+
+
+# ── Phase 3: per-cell emit (dissolve=False) ────────────────────────────────
+
+
+def test_emit_per_cell_one_polygon_per_cell():
+    cov = mortie.morton_coverage(_LATS, _LONS, order=6)
+    g = geometry.to_geometry(cov, dissolve=False)
+    assert g.geom_type == "MultiPolygon"
+    assert shapely.get_num_geometries(g) == cov.size
+
+
+def test_emit_per_cell_mixed_order_moc():
+    moc = mortie.morton_coverage_moc(_LATS, _LONS, order=8)
+    g = geometry.to_geometry(moc, dissolve=False)
+    # Each MOC cell (any order) emits exactly one quad.
+    assert shapely.get_num_geometries(g) == moc.size
+
+
+def test_emit_wkb_wkt_roundtrip_matches_cell_corners():
+    cov = mortie.morton_coverage(_LATS, _LONS, order=6)
+    wkb = geometry.to_wkb(cov, dissolve=False)
+    back = shapely.from_wkb(wkb)
+    assert shapely.get_num_geometries(back) == cov.size
+    # The first emitted cell's exterior matches mort2polygon's lon/lat corners.
+    poly0 = shapely.get_geometry(back, 0)
+    ring_lonlat = shapely.get_coordinates(shapely.get_exterior_ring(poly0))
+    want = np.array([[lon, lat] for lat, lon in mortie.mort2polygon(int(cov[0]))])
+    assert np.allclose(np.sort(ring_lonlat, axis=0), np.sort(want, axis=0))
+    # WKT path parses too.
+    assert shapely.from_wkt(geometry.to_wkt(cov, dissolve=False)).geom_type \
+        == "MultiPolygon"
+
+
+def test_emit_srid_optin_and_empty_cover():
+    cov = mortie.morton_coverage(_LATS, _LONS, order=6)
+    assert geometry.to_wkt(cov, dissolve=False, srid=4326).startswith("SRID=4326;")
+    assert int(shapely.get_srid(shapely.from_wkb(
+        geometry.to_wkb(cov, dissolve=False, srid=4326)))) == 4326
+    empty = geometry.to_geometry(np.array([], dtype=np.uint64), dissolve=False)
+    assert empty.geom_type == "MultiPolygon" and empty.is_empty
+
+
+def test_emit_dissolve_default_pending_phase4():
+    cov = mortie.morton_coverage(_LATS, _LONS, order=6)
+    with pytest.raises(NotImplementedError, match="phase 4"):
+        geometry.to_wkb(cov)
+
+
 def test_backend_gate_message(monkeypatch):
     # With no backend importable, a clear ImportError naming shapely/spherely.
     import mortie.geometry as gm

From 765f30e6320ac715adbb04673c11d81c65b8007e Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Sun, 28 Jun 2026 10:25:54 +0000
Subject: [PATCH 4/6] fold phase 3 self-review on #71

---
 mortie/tests/test_geometry.py | 49 +++++++++++++++++++++++++++++++++--
 1 file changed, 47 insertions(+), 2 deletions(-)

diff --git a/mortie/tests/test_geometry.py b/mortie/tests/test_geometry.py
index 30fe74d..021e160 100644
--- a/mortie/tests/test_geometry.py
+++ b/mortie/tests/test_geometry.py
@@ -112,7 +112,8 @@ def test_ewkb_ewkt_srid_optin():
     # EWKB carries the SRID; from_wkb reads it back.
     ewkb = geometry.geometry_to_wkb(g, srid=4326)
     assert int(shapely.get_srid(geometry.geometry_from_wkb(ewkb))) == 4326
-    assert int(shapely.get_srid(geometry.geometry_from_wkb(geometry.geometry_to_wkb(g)))) == 0
+    plain = geometry.geometry_from_wkb(geometry.geometry_to_wkb(g))
+    assert int(shapely.get_srid(plain)) == 0
 
 
 # ── Phase 2: ingest reproduces the array-path coverage ─────────────────────
@@ -242,15 +243,59 @@ def test_emit_wkb_wkt_roundtrip_matches_cell_corners():
     back = shapely.from_wkb(wkb)
     assert shapely.get_num_geometries(back) == cov.size
     # The first emitted cell's exterior matches mort2polygon's lon/lat corners.
+    # (cov is a single-order flat cover, so emit order tracks cov order here.)
     poly0 = shapely.get_geometry(back, 0)
     ring_lonlat = shapely.get_coordinates(shapely.get_exterior_ring(poly0))
     want = np.array([[lon, lat] for lat, lon in mortie.mort2polygon(int(cov[0]))])
-    assert np.allclose(np.sort(ring_lonlat, axis=0), np.sort(want, axis=0))
+    # Compare as ordered (lon, lat) PAIRS (lexsort on rows keeps the pairing, so
+    # a per-vertex lon/lat swap would be caught — column-wise sorting would not).
+    got_rows = ring_lonlat[np.lexsort((ring_lonlat[:, 1], ring_lonlat[:, 0]))]
+    want_rows = want[np.lexsort((want[:, 1], want[:, 0]))]
+    assert got_rows.shape == want_rows.shape
+    assert np.allclose(got_rows, want_rows)
     # WKT path parses too.
     assert shapely.from_wkt(geometry.to_wkt(cov, dissolve=False)).geom_type \
         == "MultiPolygon"
 
 
+def test_emit_single_cell_cover():
+    # The grp.size == 1 scalar branch of _per_cell_polygons.
+    cov = mortie.morton_coverage(_LATS, _LONS, order=6)[:1]
+    g = geometry.to_geometry(cov, dissolve=False)
+    assert shapely.get_num_geometries(g) == 1
+    assert g.is_valid
+
+
+def test_emit_step_densifies_edges():
+    cov = mortie.morton_coverage(_LATS, _LONS, order=6)[:1]
+    g1 = geometry.to_geometry(cov, dissolve=False, step=1)
+    g8 = geometry.to_geometry(cov, dissolve=False, step=8)
+    n1 = shapely.get_coordinates(shapely.get_exterior_ring(
+        shapely.get_geometry(g1, 0))).shape[0]
+    n8 = shapely.get_coordinates(shapely.get_exterior_ring(
+        shapely.get_geometry(g8, 0))).shape[0]
+    # step=1 → 4 corners (+closing); step=8 → 32 boundary points (+closing).
+    assert n1 == 5 and n8 == 33
+
+
+def test_emit_antimeridian_and_polar_cells_are_valid():
+    # A cover straddling the antimeridian and one over the north pole; per-cell
+    # emit (with mort2polygon's antimeridian normalization) must stay valid
+    # (no self-intersection) — the plan's emit acceptance criterion.
+    am = mortie.morton_coverage(
+        [10.0, 20.0, 20.0, 10.0], [179.0, 179.0, -179.0, -179.0], order=5
+    )
+    polar = mortie.morton_coverage(
+        [85.0, 85.0, 89.0, 89.0], [-90.0, 90.0, 90.0, -90.0], order=5
+    )
+    for cov in (am, polar):
+        g = geometry.to_geometry(cov, dissolve=False)
+        assert g.geom_type == "MultiPolygon" and shapely.get_num_geometries(g) > 0
+        # Every emitted cell quad is a valid (non-self-intersecting) polygon.
+        for i in range(shapely.get_num_geometries(g)):
+            assert shapely.get_geometry(g, i).is_valid
+
+
 def test_emit_srid_optin_and_empty_cover():
     cov = mortie.morton_coverage(_LATS, _LONS, order=6)
     assert geometry.to_wkt(cov, dissolve=False, srid=4326).startswith("SRID=4326;")

From 9809406553cafdec8e77cffe7663276e6bb0a8ad Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Mon, 29 Jun 2026 10:27:27 +0000
Subject: [PATCH 5/6] phase 4 of issue #71

---
 mortie/geometry.py            | 253 +++++++++++++++++++++++++++++++++-
 mortie/tests/test_geometry.py | 158 ++++++++++++++++++++-
 2 files changed, 402 insertions(+), 9 deletions(-)

diff --git a/mortie/geometry.py b/mortie/geometry.py
index 8c834b2..f37b8ce 100644
--- a/mortie/geometry.py
+++ b/mortie/geometry.py
@@ -19,6 +19,12 @@
 # Cached backend: a ``(name, module)`` pair, resolved once on first use.
 _BACKEND = None
 
+# Snap scale for vertex identity in the dissolve edge-cancellation (rounding
+# unit-vector components to 1e-10 makes a shared HEALPix corner — which both
+# adjacent cells compute identically — a single integer-keyed vertex, so their
+# shared edge cancels exactly without a floating tolerance search).
+_DISSOLVE_SNAP = 1e10
+
 # GEOS / shapely geometry type ids (shapely.get_type_id); spherely follows the
 # same numbering.  Only the ones we classify on are named.
 _TYPE_LINESTRING = 1
@@ -315,6 +321,239 @@ def _per_cell_polygons(mod, morton, step):
     return polys
 
 
+# ── emit: dissolved-boundary outline (phase 4) ─────────────────────────────
+#
+# The dissolved outline is built natively (no backend spatial predicate): every
+# cell contributes its boundary as a loop of directed edges; interior edges that
+# two adjacent cells share are traversed in opposite directions and cancel, and
+# the surviving edges chain into the outline rings.  Correctness is mortie's own
+# job throughout — the backend is only asked to *construct* the final Polygon.
+
+
+def _xyz_to_latlon(vecs):
+    """Unit vectors ``(M, 3)`` → ``(lat, lon)`` degree arrays, lon in (-180, 180]."""
+    z = np.clip(vecs[:, 2], -1.0, 1.0)
+    lat = np.degrees(np.arcsin(z))
+    lon = np.degrees(np.arctan2(vecs[:, 1], vecs[:, 0]))
+    return lat, lon
+
+
+def _spherical_signed_area(ring_xyz):
+    """Signed area (steradians) of a spherical polygon of unit vectors.
+
+    Positive = the region lies to the left of the directed boundary (an exterior
+    ring); negative = the boundary winds the other way (a hole).  Uses the
+    van Oosterom–Strackee signed-solid-angle sum over a fan from vertex 0.
+    """
+    v = ring_xyz
+    if v.shape[0] < 3:
+        return 0.0
+    a = v[0]
+    b = v[1:-1]
+    c = v[2:]
+    num = np.einsum("j,ij->i", a, np.cross(b, c))
+    den = 1.0 + b @ a + np.einsum("ij,ij->i", b, c) + c @ a
+    return float(np.sum(2.0 * np.arctan2(num, den)))
+
+
+def _boundary_rings_xyz(morton, step):
+    """Edge-cancellation dissolve of a cover → list of boundary rings.
+
+    Each ring is an ``(M, 3)`` array of unit vectors (open — first vertex not
+    repeated).  A mixed-order MOC is densified to its finest order first so every
+    cell carries unit-length edges that cancel against their neighbours.  ``step``
+    samples ``step`` points per cell edge (``step>1`` traces the curved HEALPix
+    boundary); shared sub-edge points coincide between neighbours and cancel too.
+    """
+    from . import _healpix as hp
+    from .coverage import moc_to_order
+    from .tools import _rust_mort2nested
+
+    morton = np.atleast_1d(np.asarray(morton, dtype=np.uint64))
+    if morton.size == 0:
+        return []
+
+    nested, depths = _rust_mort2nested(np.ascontiguousarray(morton))
+    udepths = np.unique(depths)
+    if udepths.size > 1:
+        morton = np.asarray(moc_to_order(morton, int(udepths.max())), dtype=np.uint64)
+        nested, depths = _rust_mort2nested(np.ascontiguousarray(morton))
+    order = int(depths[0])
+    nest = np.ascontiguousarray(nested.astype(np.int64))
+
+    bnd = hp.boundaries(order, nest, step=step)
+    if bnd.ndim == 2:
+        bnd = bnd[np.newaxis, ...]
+    pts = np.transpose(bnd, (0, 2, 1))  # (N, K, 3), K = 4*step in boundary order
+    n_cells, k = pts.shape[0], pts.shape[1]
+    flat = pts.reshape(-1, 3)
+
+    # Integer-snap every boundary point to a vertex id; a shared corner/sub-edge
+    # point collapses to one id, so adjacent cells reference the same vertex.
+    snapped = np.round(flat * _DISSOLVE_SNAP).astype(np.int64)
+    _, first_idx, inv = np.unique(
+        snapped, axis=0, return_index=True, return_inverse=True
+    )
+    id_xyz = flat[first_idx]  # representative unit vector per vertex id
+    inv = inv.reshape(n_cells, k)
+
+    # Directed edges (vertex id → vertex id) around every cell boundary.
+    starts = inv.ravel()
+    ends = np.roll(inv, -1, axis=1).ravel()
+    keep = starts != ends  # drop any degenerate zero-length edge
+    edges = list(zip(starts[keep].tolist(), ends[keep].tolist()))
+
+    # An interior edge appears as (a, b) in one cell and (b, a) in its neighbour;
+    # the surviving boundary is the net direction at each undirected edge.
+    from collections import Counter, defaultdict
+
+    counts = Counter(edges)
+    out = defaultdict(list)
+    for (a, b), c in counts.items():
+        net = c - counts.get((b, a), 0)
+        for _ in range(net):
+            out[a].append(b)
+
+    # Chain the surviving directed edges into closed rings.
+    rings = []
+    for start in list(out.keys()):
+        while out[start]:
+            cur = start
+            chain = [cur]
+            while True:
+                nxt = out[cur].pop()
+                chain.append(nxt)
+                cur = nxt
+                if cur == start:
+                    break
+            rings.append(id_xyz[np.asarray(chain[:-1])])
+    return rings
+
+
+def _count_antimeridian_crossings(lon):
+    """How many ring edges jump >180° in longitude (closed ring of lon degrees)."""
+    closed = np.concatenate([lon, lon[:1]])
+    return int(np.sum(np.abs(np.diff(closed)) > 180.0))
+
+
+def _split_at_antimeridian(coords):
+    """Split a lon/lat ring at the ±180° meridian into clean closed pieces.
+
+    ``coords`` is an open list of ``(lon, lat)`` with an even number of
+    antimeridian crossings (no enclosed pole — that case is rejected upstream).
+    Each crossing edge is cut at the meridian (latitude linearly interpolated),
+    and each resulting segment is closed along its own ±180° side.  Returns a
+    list of closed ``(lon, lat)`` rings, all within a single hemisphere of
+    longitude so the planar polygon is unambiguous.
+    """
+    n = len(coords)
+    segments = []
+    cur = []
+    for i in range(n):
+        lo0, la0 = coords[i]
+        lo1, la1 = coords[(i + 1) % n]
+        cur.append((lo0, la0))
+        if abs(lo1 - lo0) > 180.0:
+            lo1u = lo1 - 360.0 if lo1 > lo0 else lo1 + 360.0
+            boundary = 180.0 if lo1u > lo0 else -180.0
+            frac = (boundary - lo0) / (lo1u - lo0)
+            la_x = la0 + frac * (la1 - la0)
+            cur.append((boundary, la_x))
+            segments.append(cur)
+            cur = [(-boundary, la_x)]
+    if segments:
+        segments[0] = cur + segments[0]  # the wrap-around segment closes the first
+    else:
+        segments = [cur]
+    return [seg + [seg[0]] for seg in segments]
+
+
+def _point_in_ring(x, y, ring):
+    """Even-odd ray-cast point-in-polygon (``ring`` = closed list of (x, y))."""
+    inside = False
+    n = len(ring)
+    j = n - 1
+    for i in range(n):
+        xi, yi = ring[i]
+        xj, yj = ring[j]
+        if (yi > y) != (yj > y):
+            x_cross = xi + (y - yi) / (yj - yi) * (xj - xi)
+            if x < x_cross:
+                inside = not inside
+        j = i
+    return inside
+
+
+def _planar_signed_area(ring):
+    """Shoelace signed area of a closed list of ``(x, y)`` (for size ordering)."""
+    a = np.asarray(ring, dtype=np.float64)
+    x, y = a[:, 0], a[:, 1]
+    return 0.5 * float(np.sum(x * np.roll(y, -1) - np.roll(x, -1) * y))
+
+
+def _dissolved_polygons(mod, morton, step):
+    """Build the dissolved outline of *morton* as a list of backend Polygons.
+
+    Exterior and hole rings come from the edge-cancellation engine; holes are
+    nested into the exterior that contains them.  Covers whose outline encloses a
+    pole, and antimeridian-crossing holes, are rejected with a clear
+    :class:`NotImplementedError` (the spherical→planar pole/hole split is the
+    remaining sub-piece of issue #71 — see the PR thread).
+    """
+    rings_xyz = _boundary_rings_xyz(morton, step)
+    if not rings_xyz:
+        return []
+
+    # Normalise global winding: the cover's net signed area (exteriors minus
+    # holes) is the covered area, always positive.  HEALPix orders boundary
+    # points one way for step==1 and the other for step>1, so key the
+    # exterior/hole sign off this invariant rather than a fixed convention.
+    areas = [_spherical_signed_area(r) for r in rings_xyz]
+    if sum(areas) < 0.0:
+        rings_xyz = [r[::-1] for r in rings_xyz]
+        areas = [-a for a in areas]
+
+    ext_pieces = []
+    holes = []
+    for ring, area in zip(rings_xyz, areas):
+        lat, lon = _xyz_to_latlon(ring)
+        crossings = _count_antimeridian_crossings(lon)
+        if crossings % 2 == 1:
+            raise NotImplementedError(
+                "dissolved emit of a pole-enclosing cover (e.g. a polar cap) is "
+                "not yet supported; pass dissolve=False for the per-cell "
+                "MultiPolygon (issue #71 phase 4 follow-up)"
+            )
+        ll = list(zip(lon.tolist(), lat.tolist()))
+        if area >= 0.0:
+            ext_pieces.extend(_split_at_antimeridian(ll))
+        elif crossings > 0:
+            raise NotImplementedError(
+                "dissolved emit with an antimeridian-crossing hole is not yet "
+                "supported; pass dissolve=False (issue #71 phase 4 follow-up)"
+            )
+        else:
+            holes.append(ll + [ll[0]])
+
+    # Nest each hole into the smallest exterior piece that contains it.
+    hole_groups = [[] for _ in ext_pieces]
+    ext_areas = [abs(_planar_signed_area(p)) for p in ext_pieces]
+    for hole in holes:
+        cx = float(np.mean([p[0] for p in hole]))
+        cy = float(np.mean([p[1] for p in hole]))
+        best = None
+        for idx, piece in enumerate(ext_pieces):
+            if _point_in_ring(cx, cy, piece) and (
+                best is None or ext_areas[idx] < ext_areas[best]
+            ):
+                best = idx
+        hole_groups[best if best is not None else 0].append(hole)
+
+    return [
+        mod.Polygon(ext_pieces[i], hole_groups[i]) for i in range(len(ext_pieces))
+    ]
+
+
 def to_geometry(morton, dissolve=True, step=1):
     """Convert a morton cover to a backend geometry (issue #71).
 
@@ -324,7 +563,8 @@ def to_geometry(morton, dissolve=True, step=1):
         A morton cover (flat or mixed-order MOC; each word self-encodes order).
     dissolve : bool, optional
         ``True`` (default) emits the single dissolved outline of the whole cover
-        (exterior rings, holes, and disjoint components).  ``False`` emits a
+        (exterior rings, holes, and disjoint components), built natively by
+        edge-cancellation — no backend spatial predicate.  ``False`` emits a
         per-cell ``MultiPolygon`` — one quad per cell.
     step : int, optional
         Boundary points per cell edge (default 1 = 4 corners / straight chords).
@@ -337,14 +577,15 @@ def to_geometry(morton, dissolve=True, step=1):
 
     Notes
     -----
-    Emit requires the shapely backend (it constructs geometry objects).
+    Emit requires the shapely backend (it constructs geometry objects).  The
+    dissolved emit (``dissolve=True``) does not yet support covers whose outline
+    encloses a pole or holes that cross the antimeridian — both raise
+    :class:`NotImplementedError`; use ``dissolve=False`` for those (issue #71
+    phase 4 follow-up).
     """
     mod = _require_shapely("geometry emit")
     if dissolve:
-        raise NotImplementedError(
-            "dissolve=True (the dissolved-outline emit) lands in phase 4 of "
-            "issue #71; pass dissolve=False for the per-cell MultiPolygon"
-        )
+        return mod.MultiPolygon(_dissolved_polygons(mod, morton, step))
     return mod.MultiPolygon(_per_cell_polygons(mod, morton, step))
 
 
diff --git a/mortie/tests/test_geometry.py b/mortie/tests/test_geometry.py
index 021e160..4591800 100644
--- a/mortie/tests/test_geometry.py
+++ b/mortie/tests/test_geometry.py
@@ -305,10 +305,17 @@ def test_emit_srid_optin_and_empty_cover():
     assert empty.geom_type == "MultiPolygon" and empty.is_empty
 
 
-def test_emit_dissolve_default_pending_phase4():
+def test_emit_dissolve_is_the_default():
+    # dissolve=True is the default: to_wkb(cov) with no flag emits the dissolved
+    # outline (one ring for a contiguous box), not the per-cell quads.
     cov = mortie.morton_coverage(_LATS, _LONS, order=6)
-    with pytest.raises(NotImplementedError, match="phase 4"):
-        geometry.to_wkb(cov)
+    dissolved = shapely.from_wkb(geometry.to_wkb(cov))
+    assert dissolved.geom_type == "MultiPolygon"
+    assert shapely.get_num_geometries(dissolved) == 1
+    # The dissolved outline has far fewer vertices than the per-cell emit.
+    per_cell = geometry.to_geometry(cov, dissolve=False)
+    assert shapely.get_num_coordinates(dissolved) < \
+        shapely.get_num_coordinates(per_cell)
 
 
 def test_backend_gate_message(monkeypatch):
@@ -327,3 +334,148 @@ def _block(name, *args, **kwargs):
     with pytest.raises(ImportError, match="shapely"):
         gm._require_backend()
     # monkeypatch reverts _BACKEND and __import__; the next call re-resolves.
+
+
+# ── dissolved-outline emit (phase 4) ───────────────────────────────────────
+
+
+def _union_oracle(cov):
+    """shapely.unary_union of the per-cell quads — the independent dissolve
+    reference, valid away from the antimeridian / poles (planar union)."""
+    polys = [
+        shapely.Polygon([(lon, lat) for lat, lon in mortie.mort2polygon(int(w))])
+        for w in cov
+    ]
+    return shapely.unary_union(polys)
+
+
+def _ring_spherical_area(coords):
+    """Signed steradian area of a closed lon/lat ring (backend-independent — the
+    one oracle that still works across the antimeridian / poles)."""
+    a = np.asarray(coords[:-1], dtype=np.float64)
+    rlat = np.radians(a[:, 1])
+    rlon = np.radians(a[:, 0])
+    v = np.column_stack(
+        [np.cos(rlat) * np.cos(rlon), np.cos(rlat) * np.sin(rlon), np.sin(rlat)]
+    )
+    p0, b, c = v[0], v[1:-1], v[2:]
+    num = np.einsum("j,ij->i", p0, np.cross(b, c))
+    den = 1.0 + b @ p0 + np.einsum("ij,ij->i", b, c) + c @ p0
+    return float(np.sum(2.0 * np.arctan2(num, den)))
+
+
+def _cover_area(cov):
+    _, depths = mortie.tools._rust_mort2nested(
+        np.ascontiguousarray(np.asarray(cov, dtype=np.uint64))
+    )
+    return float(np.sum(4.0 * np.pi / (12.0 * (4.0 ** depths))))
+
+
+def test_dissolve_matches_union_oracle():
+    cov = mortie.morton_coverage(_LATS, _LONS, order=6)
+    mp = geometry.to_geometry(cov)
+    assert mp.is_valid and shapely.get_num_geometries(mp) == 1
+    # Native edge-cancellation dissolve == shapely's planar union, to precision.
+    assert mp.symmetric_difference(_union_oracle(cov)).area < 1e-9
+
+
+def test_dissolve_polygon_with_hole():
+    big = mortie.morton_coverage(
+        [30.0, 30.0, 60.0, 60.0], [-130.0, -100.0, -100.0, -130.0], order=5
+    )
+    inner = mortie.morton_coverage(
+        [42.0, 42.0, 48.0, 48.0], [-120.0, -110.0, -110.0, -120.0], order=5
+    )
+    cov = np.array(sorted(set(big.tolist()) - set(inner.tolist())), dtype=np.uint64)
+    mp = geometry.to_geometry(cov)
+    assert mp.is_valid and shapely.get_num_geometries(mp) == 1
+    # The carved-out interior is emitted as exactly one hole.
+    assert shapely.get_num_interior_rings(shapely.get_geometry(mp, 0)) == 1
+    assert mp.symmetric_difference(_union_oracle(cov)).area < 1e-9
+
+
+def test_dissolve_disjoint_components():
+    a = mortie.morton_coverage(
+        [40.0, 40.0, 45.0, 45.0], [-120.0, -115.0, -115.0, -120.0], order=6
+    )
+    b = mortie.morton_coverage(
+        [40.0, 40.0, 45.0, 45.0], [-100.0, -95.0, -95.0, -100.0], order=6
+    )
+    cov = np.unique(np.concatenate([a, b]))
+    mp = geometry.to_geometry(cov)
+    assert mp.is_valid and shapely.get_num_geometries(mp) == 2
+    assert mp.symmetric_difference(_union_oracle(cov)).area < 1e-9
+
+
+def test_dissolve_step_densifies_and_matches():
+    cov = mortie.morton_coverage(_LATS, _LONS, order=6)
+    g1 = geometry.to_geometry(cov, step=1)
+    g4 = geometry.to_geometry(cov, step=4)
+    # step traces the curved HEALPix edge: more boundary vertices, still valid.
+    # (The curved outline differs from the straight-chord union by the chord
+    # error, so area conservation — not planar symdiff — is the right oracle.)
+    assert g4.is_valid
+    assert shapely.get_num_coordinates(g4) > shapely.get_num_coordinates(g1)
+    ring = list(shapely.get_coordinates(
+        shapely.get_exterior_ring(shapely.get_geometry(g4, 0))))
+    assert abs(_ring_spherical_area(ring) - _cover_area(cov)) < 1e-3
+
+
+def test_dissolve_mixed_order_moc():
+    moc = mortie.morton_coverage_moc(_LATS, _LONS, order=8)
+    mp = geometry.to_geometry(moc)
+    assert mp.is_valid and shapely.get_num_geometries(mp) == 1
+    # The MOC is densified to its finest order, so the dissolved outline encloses
+    # the same area as the MOC's cells (independent of the planar union, which a
+    # coarse-vs-fine chord mismatch would perturb).
+    assert abs(_ring_spherical_area(
+        list(shapely.get_coordinates(shapely.get_exterior_ring(
+            shapely.get_geometry(mp, 0))))
+    ) - _cover_area(mortie.coverage.moc_to_order(moc, 8))) < 1e-3
+
+
+def test_dissolve_antimeridian_split():
+    # A box straddling +/-180: the outline crosses the antimeridian (an even
+    # number of times, no pole) and must split into valid pieces.
+    cov = mortie.morton_coverage(
+        [10.0, 10.0, 20.0, 20.0], [170.0, -170.0, -170.0, 170.0], order=5
+    )
+    mp = geometry.to_geometry(cov)
+    assert mp.is_valid and shapely.get_num_geometries(mp) >= 2
+    # No emitted piece spans more than a hemisphere of longitude (the hallmark of
+    # a correctly split antimeridian polygon).
+    out_area = 0.0
+    for i in range(shapely.get_num_geometries(mp)):
+        ring = list(shapely.get_coordinates(
+            shapely.get_exterior_ring(shapely.get_geometry(mp, i))))
+        lons = np.asarray(ring)[:, 0]
+        assert lons.max() - lons.min() <= 180.0 + 1e-9
+        out_area += _ring_spherical_area(ring)
+    # The split conserves the covered area exactly (no sliver lost at the seam).
+    assert abs(out_area - _cover_area(cov)) < 1e-3
+
+
+def test_dissolve_pole_enclosing_raises():
+    # A ring of cells around the north pole encloses it (odd antimeridian
+    # crossings) — not yet supported; must raise a clear, specific error.
+    lats, lons = [], []
+    for lo in range(-180, 180, 20):
+        lats += [82.0, 82.0, 89.9, 89.9]
+        lons += [lo, lo + 20, lo + 20, lo]
+    cap = np.unique(mortie.morton_coverage(lats, lons, order=4))
+    with pytest.raises(NotImplementedError, match="pole-enclosing"):
+        geometry.to_geometry(cap)
+
+
+def test_dissolve_wkb_wkt_srid_roundtrip():
+    cov = mortie.morton_coverage(_LATS, _LONS, order=6)
+    assert shapely.from_wkb(geometry.to_wkb(cov)).geom_type == "MultiPolygon"
+    assert shapely.from_wkt(geometry.to_wkt(cov)).geom_type == "MultiPolygon"
+    assert int(shapely.get_srid(
+        shapely.from_wkb(geometry.to_wkb(cov, srid=4326)))) == 4326
+    assert geometry.to_wkt(cov, srid=4326).startswith("SRID=4326;")
+
+
+def test_dissolve_empty_cover():
+    mp = geometry.to_geometry(np.array([], dtype=np.uint64))
+    assert mp.geom_type == "MultiPolygon" and mp.is_empty

From d1729a54d29248b89d44d3710fecbdd83cfc973b Mon Sep 17 00:00:00 2001
From: Claude <noreply@anthropic.com>
Date: Mon, 29 Jun 2026 10:40:37 +0000
Subject: [PATCH 6/6] fold phase 4 self-review on #71

---
 mortie/geometry.py            | 160 ++++++++++++++++++++++++----------
 mortie/tests/test_geometry.py |  39 +++++++++
 2 files changed, 154 insertions(+), 45 deletions(-)

diff --git a/mortie/geometry.py b/mortie/geometry.py
index f37b8ce..bff0c2a 100644
--- a/mortie/geometry.py
+++ b/mortie/geometry.py
@@ -14,6 +14,8 @@
 module flips the axes at the boundary and works in degrees throughout.
 """
 
+import math
+
 import numpy as np
 
 # Cached backend: a ``(name, module)`` pair, resolved once on first use.
@@ -405,46 +407,95 @@ def _boundary_rings_xyz(morton, step):
 
     # An interior edge appears as (a, b) in one cell and (b, a) in its neighbour;
     # the surviving boundary is the net direction at each undirected edge.
-    from collections import Counter, defaultdict
+    from collections import Counter
 
     counts = Counter(edges)
-    out = defaultdict(list)
+    survivors = []
     for (a, b), c in counts.items():
         net = c - counts.get((b, a), 0)
-        for _ in range(net):
-            out[a].append(b)
+        survivors.extend([(a, b)] * net)
+    return _chain_rings(survivors, id_xyz)
+
+
+def _tangent_azimuth(p, q):
+    """Azimuth (radians) from unit vector *p* toward unit vector *q*, in p's
+    tangent plane (north-referenced).  Used to order edges around a vertex."""
+    d = q - np.dot(q, p) * p
+    nd = np.linalg.norm(d)
+    if nd < 1e-15:
+        return 0.0
+    d = d / nd
+    east = np.cross([0.0, 0.0, 1.0], p)
+    ne = np.linalg.norm(east)
+    east = np.array([1.0, 0.0, 0.0]) if ne < 1e-9 else east / ne
+    north = np.cross(p, east)
+    return math.atan2(float(np.dot(d, east)), float(np.dot(d, north)))
+
+
+def _chain_rings(survivors, id_xyz):
+    """Chain surviving directed boundary edges into closed rings.
+
+    At a non-manifold vertex (out-degree > 1 — e.g. two cells touching only at a
+    corner) the next edge is chosen by angular order: the surviving edge whose
+    departure azimuth is the smallest turn anticlockwise from the reversed
+    arrival direction.  This right-hand-rule traversal yields *simple* rings
+    (the cells' boundaries stay separate rather than crossing into a bowtie),
+    independent of the cover's global winding.
+    """
+    from collections import defaultdict
+
+    az = {e: _tangent_azimuth(id_xyz[e[0]], id_xyz[e[1]]) for e in survivors}
+    records = [[a, b, True] for a, b in survivors]
+    by_start = defaultdict(list)
+    for rec in records:
+        by_start[rec[0]].append(rec)
 
-    # Chain the surviving directed edges into closed rings.
     rings = []
-    for start in list(out.keys()):
-        while out[start]:
-            cur = start
-            chain = [cur]
-            while True:
-                nxt = out[cur].pop()
-                chain.append(nxt)
-                cur = nxt
-                if cur == start:
-                    break
-            rings.append(id_xyz[np.asarray(chain[:-1])])
+    for seed in records:
+        if not seed[2]:
+            continue
+        seed_start = seed[0]
+        cur = seed
+        chain = []
+        while cur is not None and cur[2]:
+            cur[2] = False
+            chain.append(cur[0])
+            v = cur[1]
+            if v == seed_start:
+                break  # returned to the start vertex — ring closed
+            cand = [r for r in by_start[v] if r[2]]
+            if not cand:
+                break
+            if len(cand) == 1:
+                cur = cand[0]
+            else:
+                # Smallest turn anticlockwise from the reversed arrival keeps the
+                # walk on the same face (no crossing) at a non-manifold vertex.
+                back = _tangent_azimuth(id_xyz[v], id_xyz[cur[0]])
+                cur = min(cand, key=lambda r: (az[(r[0], r[1])] - back) % (2 * math.pi))
+        rings.append(id_xyz[np.asarray(chain)])
     return rings
 
 
-def _count_antimeridian_crossings(lon):
-    """How many ring edges jump >180° in longitude (closed ring of lon degrees)."""
-    closed = np.concatenate([lon, lon[:1]])
-    return int(np.sum(np.abs(np.diff(closed)) > 180.0))
+def _antimeridian_winding(lon):
+    """Net signed longitude winding (degrees) and antimeridian-crossing count of
+    a closed ring of longitudes.  Net ≈ ±360 ⟺ the ring encircles a pole."""
+    deltas = np.diff(np.concatenate([lon, lon[:1]]))
+    crossings = int(np.sum(np.abs(deltas) > 180.0))
+    net = float(np.sum((deltas + 180.0) % 360.0 - 180.0))
+    return net, crossings
 
 
 def _split_at_antimeridian(coords):
-    """Split a lon/lat ring at the ±180° meridian into clean closed pieces.
-
-    ``coords`` is an open list of ``(lon, lat)`` with an even number of
-    antimeridian crossings (no enclosed pole — that case is rejected upstream).
-    Each crossing edge is cut at the meridian (latitude linearly interpolated),
-    and each resulting segment is closed along its own ±180° side.  Returns a
-    list of closed ``(lon, lat)`` rings, all within a single hemisphere of
-    longitude so the planar polygon is unambiguous.
+    """Split a lon/lat ring crossing the ±180° meridian exactly twice into two
+    clean closed pieces.
+
+    ``coords`` is an open list of ``(lon, lat)`` with exactly two antimeridian
+    crossings (the only case routed here — pole-enclosing rings and rings with
+    more than two crossings are rejected upstream).  Each crossing edge is cut at
+    the meridian (latitude linearly interpolated) and each segment is closed
+    along its own ±180° side, so both pieces sit within a single hemisphere of
+    longitude and the planar polygon is unambiguous.
     """
     n = len(coords)
     segments = []
@@ -495,10 +546,11 @@ def _dissolved_polygons(mod, morton, step):
     """Build the dissolved outline of *morton* as a list of backend Polygons.
 
     Exterior and hole rings come from the edge-cancellation engine; holes are
-    nested into the exterior that contains them.  Covers whose outline encloses a
-    pole, and antimeridian-crossing holes, are rejected with a clear
-    :class:`NotImplementedError` (the spherical→planar pole/hole split is the
-    remaining sub-piece of issue #71 — see the PR thread).
+    nested into the exterior that contains them.  The following are rejected with
+    a clear :class:`NotImplementedError` — the spherical→planar pole/hole split
+    is the remaining sub-piece of issue #71 (see the PR thread): covers whose
+    outline encloses a pole, an exterior crossing the antimeridian more than
+    twice, and antimeridian-crossing holes.
     """
     rings_xyz = _boundary_rings_xyz(morton, step)
     if not rings_xyz:
@@ -508,6 +560,8 @@ def _dissolved_polygons(mod, morton, step):
     # holes) is the covered area, always positive.  HEALPix orders boundary
     # points one way for step==1 and the other for step>1, so key the
     # exterior/hole sign off this invariant rather than a fixed convention.
+    # (Spherical signed area is defined mod 4π, so this assumes the cover stays
+    # well under a hemisphere — true for every realistic emit input.)
     areas = [_spherical_signed_area(r) for r in rings_xyz]
     if sum(areas) < 0.0:
         rings_xyz = [r[::-1] for r in rings_xyz]
@@ -517,37 +571,53 @@ def _dissolved_polygons(mod, morton, step):
     holes = []
     for ring, area in zip(rings_xyz, areas):
         lat, lon = _xyz_to_latlon(ring)
-        crossings = _count_antimeridian_crossings(lon)
-        if crossings % 2 == 1:
+        net_winding, crossings = _antimeridian_winding(lon)
+        if abs(net_winding) > 180.0:  # net ≈ ±360° ⟺ the ring encircles a pole
             raise NotImplementedError(
                 "dissolved emit of a pole-enclosing cover (e.g. a polar cap) is "
                 "not yet supported; pass dissolve=False for the per-cell "
                 "MultiPolygon (issue #71 phase 4 follow-up)"
             )
         ll = list(zip(lon.tolist(), lat.tolist()))
-        if area >= 0.0:
+        if area < 0.0:
+            if crossings > 0:
+                raise NotImplementedError(
+                    "dissolved emit with an antimeridian-crossing hole is not "
+                    "yet supported; pass dissolve=False (issue #71 phase 4 "
+                    "follow-up)"
+                )
+            holes.append(ll + [ll[0]])
+        elif crossings == 0:
+            ext_pieces.append(ll + [ll[0]])
+        elif crossings == 2:
             ext_pieces.extend(_split_at_antimeridian(ll))
-        elif crossings > 0:
+        else:
             raise NotImplementedError(
-                "dissolved emit with an antimeridian-crossing hole is not yet "
-                "supported; pass dissolve=False (issue #71 phase 4 follow-up)"
+                "dissolved emit of an exterior crossing the antimeridian more "
+                "than twice is not yet supported; pass dissolve=False (issue "
+                "#71 phase 4 follow-up)"
             )
-        else:
-            holes.append(ll + [ll[0]])
 
-    # Nest each hole into the smallest exterior piece that contains it.
+    # Nest each hole into the smallest exterior piece that contains it.  A hole
+    # vertex lies strictly inside its surrounding exterior, so test a vertex
+    # (a guaranteed-interior point) rather than the centroid, which a concave or
+    # split ring can push outside the region.
     hole_groups = [[] for _ in ext_pieces]
     ext_areas = [abs(_planar_signed_area(p)) for p in ext_pieces]
     for hole in holes:
-        cx = float(np.mean([p[0] for p in hole]))
-        cy = float(np.mean([p[1] for p in hole]))
+        hx, hy = hole[0]
         best = None
         for idx, piece in enumerate(ext_pieces):
-            if _point_in_ring(cx, cy, piece) and (
+            if _point_in_ring(hx, hy, piece) and (
                 best is None or ext_areas[idx] < ext_areas[best]
             ):
                 best = idx
-        hole_groups[best if best is not None else 0].append(hole)
+        if best is None:
+            raise NotImplementedError(
+                "dissolved emit could not nest a hole into any exterior (an "
+                "unsupported self-touching outline); pass dissolve=False"
+            )
+        hole_groups[best].append(hole)
 
     return [
         mod.Polygon(ext_pieces[i], hole_groups[i]) for i in range(len(ext_pieces))
diff --git a/mortie/tests/test_geometry.py b/mortie/tests/test_geometry.py
index 4591800..4a41d91 100644
--- a/mortie/tests/test_geometry.py
+++ b/mortie/tests/test_geometry.py
@@ -479,3 +479,42 @@ def test_dissolve_wkb_wkt_srid_roundtrip():
 def test_dissolve_empty_cover():
     mp = geometry.to_geometry(np.array([], dtype=np.uint64))
     assert mp.geom_type == "MultiPolygon" and mp.is_empty
+
+
+def _interleave(x, y, order):
+    h = 0
+    for i in range(order):
+        h |= ((x >> i) & 1) << (2 * i)
+        h |= ((y >> i) & 1) << (2 * i + 1)
+    return h
+
+
+def test_dissolve_corner_touch_yields_simple_rings():
+    # Two cells in one base face that touch ONLY at a corner (non-manifold
+    # boundary vertex).  Angular ring-chaining must keep them as two separate
+    # valid polygons rather than crossing into a self-touching bowtie.
+    from mortie import _rustie
+
+    order, face, side = 4, 4, 2 ** 4
+    nested = np.array(
+        [face * side * side + _interleave(5, 5, order),
+         face * side * side + _interleave(6, 6, order)],
+        dtype=np.uint64,
+    )
+    cov = np.unique(np.asarray(
+        _rustie.rust_mi_from_nested(nested, order), dtype=np.uint64))
+    mp = geometry.to_geometry(cov)
+    assert mp.is_valid and shapely.get_num_geometries(mp) == 2
+    for i in range(2):
+        assert shapely.get_geometry(mp, i).is_valid
+
+
+def test_dissolve_step_cancels_seams_no_spurious_holes():
+    # With step>1 the shared sub-edge points between neighbours must still cancel,
+    # or failed cancellation would leave sliver interior rings.  A contiguous box
+    # must dissolve to a single hole-free outline at every step.
+    cov = mortie.morton_coverage(_LATS, _LONS, order=6)
+    for step in (2, 4, 8):
+        mp = geometry.to_geometry(cov, step=step)
+        assert mp.is_valid and shapely.get_num_geometries(mp) == 1
+        assert shapely.get_num_interior_rings(shapely.get_geometry(mp, 0)) == 0