Group A small-fix bundle from issue #108

mortie/geometry.py

@@ -226,8 +226,9 @@ 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`` and for
-        linear geometry.  Note ``morton_coverage_moc`` has no orientation
+        (see :func:`mortie.morton_coverage`).  Ignored when ``moc=True``;
+        ``normalize=False`` with linear geometry raises ``ValueError`` (a line
+        has no ring orientation).  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
@@ -260,6 +261,8 @@ def from_geometry(geom, order=18, moc=False, normalize=True,
         raise ValueError(
             "moc / tolerance / max_cells apply only to polygonal geometry"
         )
+    if not normalize:
+        raise ValueError("normalize applies 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]
@@ -644,6 +647,33 @@ def _ring_signed_area_lonlat(ring):
     return _spherical_signed_area(v)
 
 
+def _reject_hemisphere_cover(morton):
+    """Mirror of the Rust hemisphere guard (``src_rust/src/dissolve.rs``, issue
+    #108): exterior/hole classification keys off the sign of the mod-4π
+    spherical signed area, which is ambiguous once the cover nears 2π — fail
+    loud on the exact covered area (Σ π/(3·4^depth), cells are equal-area)
+    instead of silently swapping shells and holes.  Assumes disjoint,
+    non-duplicated cells (the dissolve precondition anyway — duplicate words
+    would break edge cancellation); duplicates double-count.  Returns the
+    exact covered area (steradians) for the wrap cross-check downstream.
+    """
+    from .tools import _rust_mort2nested
+
+    morton = np.atleast_1d(np.asarray(morton, dtype=np.uint64))
+    if morton.size == 0:
+        return 0.0
+    _, depths = _rust_mort2nested(np.ascontiguousarray(morton))
+    area = float(np.sum(np.pi / (3.0 * 4.0 ** depths.astype(np.float64))))
+    if area > 2.0 * np.pi * 0.98:
+        raise ValueError(
+            f"dissolved cover spans {area:.6f} sr — within 2% of a hemisphere "
+            "(2π sr) or beyond — so its exterior/hole winding is ambiguous; "
+            "split the cover into sub-hemisphere parts or pass dissolve=False "
+            "for per-cell polygons"
+        )
+    return area
+
+
 def _dissolved_rings_py(morton, step):
     """Reference (pure-Python) dissolve → ``(ext_pieces, holes)`` lon/lat rings.
 
@@ -655,6 +685,7 @@ def _dissolved_rings_py(morton, step):
     explicit ±90° pole vertices for a pole-enclosing region.  Each returned ring
     is a closed list of ``(lon, lat)`` degree pairs.
     """
+    cover_area = _reject_hemisphere_cover(morton)
     rings_xyz = _boundary_rings_xyz(morton, step)
     if not rings_xyz:
         return [], []
@@ -663,10 +694,23 @@ def _dissolved_rings_py(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.)
+    # The fan formula wraps mod 4π when a single ring encloses more than a
+    # hemisphere (possible even for a small cover, e.g. an equatorial band),
+    # which would flip the sign here; an honest |Σ| matches the exact covered
+    # area to within chord discretization (≲0.1 sr at step==1) while any wrap
+    # is off by ~4π, so a π tolerance separates them cleanly (mirrors the Rust
+    # cross-check in `src_rust/src/dissolve.rs::classify_and_split`).
     areas = [_spherical_signed_area(r) for r in rings_xyz]
-    if sum(areas) < 0.0:
+    total = sum(areas)
+    if abs(abs(total) - cover_area) > np.pi:
+        raise ValueError(
+            "dissolved cover has a boundary ring enclosing more than a "
+            f"hemisphere (|Σ ring areas| = {abs(total):.6f} sr vs covered area "
+            f"{cover_area:.6f} sr), so its exterior/hole winding cannot be "
+            "classified; split the cover into sub-hemisphere parts or pass "
+            "dissolve=False for per-cell polygons"
+        )
+    if total < 0.0:
         rings_xyz = [r[::-1] for r in rings_xyz]
         areas = [-a for a in areas]
 
@@ -785,7 +829,11 @@ def to_geometry(morton, dissolve=True, step=1):
     caps), exteriors crossing the antimeridian any even number of times, and
     antimeridian-crossing holes: crossing rings are cut at ±180° and reconnected
     by the GeoJSON convention — a single split ``MultiPolygon`` with explicit
-    ±90° pole vertices stitched down the antimeridian.
+    ±90° pole vertices stitched down the antimeridian.  A cover spanning near
+    or over a hemisphere (2π sr), or one with a boundary ring enclosing more
+    than a hemisphere (e.g. an equatorial band), raises ``ValueError`` — its
+    exterior/hole winding is ambiguous (issue #108); split such a cover or use
+    ``dissolve=False``.
     """
     mod = _require_shapely("geometry emit")
     if dissolve:

mortie/tests/test_geometry.py

@@ -206,6 +206,15 @@ def test_ingest_linear_rejects_polygon_only_args():
         mortie.from_wkt(wkt, order=6, moc=True)
 
 
+def test_ingest_linear_rejects_normalize_false():
+    # A line has no ring orientation: silently ignoring normalize=False would
+    # be a no-op, so the linear branch rejects it (issue #108).
+    for wkt in ("LINESTRING (0 0, 1 1, 2 0)",
+                "MULTILINESTRING ((0 0, 1 1), (2 2, 3 3))"):
+        with pytest.raises(ValueError, match="only to polygonal"):
+            mortie.from_wkt(wkt, order=6, normalize=False)
+
+
 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)
@@ -581,6 +590,40 @@ def test_dissolve_empty_cover():
     assert mp.geom_type == "MultiPolygon" and mp.is_empty
 
 
+def test_dissolve_hemisphere_cover_fails_loud():
+    # The dissolve keys exterior/hole off the sign of the mod-4π spherical
+    # signed area, ambiguous once the cover nears a hemisphere (2π sr): the
+    # guard on the exact covered area raises instead of silently swapping
+    # shells and holes (issue #108).  Rust runtime path and Python oracle agree.
+    # 24 order-1 cells (base cells 0-5) tile exactly half the sphere; the polar
+    # cap reaching past the equator is a polar-scale cover beyond it.
+    hemi = mortie.norm2mort(np.tile(np.arange(4), 6), np.repeat(np.arange(6), 4), 1)
+    over = _polar_cap(-2.0, 89.9, order=3)
+    for cov in (hemi, over):
+        with pytest.raises(ValueError, match="hemisphere"):
+            geometry.to_geometry(cov)
+        with pytest.raises(ValueError, match="hemisphere"):
+            geometry._dissolved_rings_py(cov, 1)
+    # The documented fallback stays available: per-cell emit, one quad per cell.
+    per_cell = geometry.to_geometry(hemi, dissolve=False)
+    assert shapely.get_num_geometries(per_cell) == hemi.size
+
+
+def test_dissolve_hemisphere_enclosing_ring_fails_loud():
+    # A thin equatorial band (~1.3 sr, far under the covered-area guard) has
+    # two boundary rings that each enclose more than a hemisphere, so the
+    # mod-4π fan sum wraps (Σ = A − 4π < 0) — without the Σ-vs-exact-area
+    # cross-check the global flip fires on a correctly-wound cover and the
+    # antimeridian stitcher dies with a cryptic error (issue #108 review).
+    lats, lons = np.meshgrid(np.arange(-3.5, 3.6, 0.5),
+                             np.arange(-180.0, 180.0, 1.0))
+    band = np.unique(mortie.geo2mort(lats.ravel(), lons.ravel(), order=4))
+    with pytest.raises(ValueError, match="enclosing more than a hemisphere"):
+        geometry.to_geometry(band)
+    with pytest.raises(ValueError, match="enclosing more than a hemisphere"):
+        geometry._dissolved_rings_py(band, 1)
+
+
 def _interleave(x, y, order):
     h = 0
     for i in range(order):

src_rust/src/dissolve.rs

@@ -41,10 +41,43 @@ pub struct ClassifiedRings {
     pub holes: Vec<Ring>,
 }
 
+// Hemisphere guard (issue #108): exterior/hole classification keys off the
+// sign of Σ ring signed areas, which is defined mod 4π — at 2π a cover reads
+// the same as its complement wound the other way, and past 2π the sign
+// silently inverts (the fan formula can also wrap per-ring at that scale).
+// The covered area itself is exact (equal-area cells), so gate on it: 2% of 2π
+// keeps a comfortable distance from the breakdown point while excluding only
+// covers within 2% of half the sphere.
+const HEMISPHERE_MARGIN: f64 = std::f64::consts::TAU * 0.02;
+
+/// Exact covered area (steradians) of a morton cover: Σ π/(3·4^depth).
+/// Assumes disjoint, non-duplicated cells (the dissolve precondition anyway —
+/// duplicate words would break edge cancellation); duplicates double-count.
+fn cover_area(morton: &[u64]) -> f64 {
+    morton
+        .iter()
+        .map(|&w| std::f64::consts::PI / (3.0 * 4f64.powi(mort2nested(w).1 as i32)))
+        .sum()
+}
+
 /// Dissolve a morton cover into classified planar (lon, lat) rings.
-pub fn dissolve(morton: &[u64], step: u32) -> ClassifiedRings {
+///
+/// Errs (with a message for a Python `ValueError`) when the cover spans
+/// near or over a hemisphere, or when a boundary ring encloses more than a
+/// hemisphere (the mod-4π fan sum wraps) — both make the winding-sign
+/// normalisation of [`classify_and_split`] untrustworthy.
+pub fn dissolve(morton: &[u64], step: u32) -> Result<ClassifiedRings, String> {
+    let area = cover_area(morton);
+    if area > std::f64::consts::TAU - HEMISPHERE_MARGIN {
+        return Err(format!(
+            "dissolved cover spans {area:.6} sr — within 2% of a hemisphere \
+             (2π sr) or beyond — so its exterior/hole winding is ambiguous; \
+             split the cover into sub-hemisphere parts or pass dissolve=False \
+             for per-cell polygons"
+        ));
+    }
     let rings = boundary_rings_xyz(morton, step);
-    classify_and_split(rings)
+    classify_and_split(rings, area)
 }
 
 // ── edge-cancellation: cover → boundary rings (unit vectors) ───────────────
@@ -387,18 +420,37 @@ fn ring_signed_area_lonlat(ring: &[(f64, f64)]) -> f64 {
     spherical_signed_area(&v)
 }
 
-fn classify_and_split(mut rings_xyz: Vec<Vec<Vec3>>) -> ClassifiedRings {
+fn classify_and_split(
+    mut rings_xyz: Vec<Vec<Vec3>>,
+    cover_area: f64,
+) -> Result<ClassifiedRings, String> {
     let mut out = ClassifiedRings {
         shells: Vec::new(),
         holes: Vec::new(),
     };
     if rings_xyz.is_empty() {
-        return out;
+        return Ok(out);
     }
     // Normalise global winding so the covered area (exteriors minus holes) is
     // positive — the boundary point order differs between step==1 and step>1.
+    // The fan formula wraps mod 4π when a single ring encloses more than a
+    // hemisphere (possible even for a small cover, e.g. an equatorial band),
+    // which would flip the sign here; an honest |Σ| matches the exact covered
+    // area to within chord discretization (≲0.1 sr at step==1) while any wrap
+    // is off by ~4π, so a π tolerance separates them cleanly.
     let mut areas: Vec<f64> = rings_xyz.iter().map(|r| spherical_signed_area(r)).collect();
-    if areas.iter().sum::<f64>() < 0.0 {
+    let total: f64 = areas.iter().sum();
+    if (total.abs() - cover_area).abs() > std::f64::consts::PI {
+        return Err(format!(
+            "dissolved cover has a boundary ring enclosing more than a \
+             hemisphere (|Σ ring areas| = {:.6} sr vs covered area {cover_area:.6} \
+             sr), so its exterior/hole winding cannot be classified; split the \
+             cover into sub-hemisphere parts or pass dissolve=False for \
+             per-cell polygons",
+            total.abs()
+        ));
+    }
+    if total < 0.0 {
         for r in rings_xyz.iter_mut() {
             r.reverse();
         }
@@ -442,7 +494,7 @@ fn classify_and_split(mut rings_xyz: Vec<Vec<Vec3>>) -> ClassifiedRings {
             }
         }
     }
-    out
+    Ok(out)
 }
 
 #[cfg(test)]
@@ -487,6 +539,50 @@ mod tests {
         assert!(net_winding(&ring).abs() > 180.0);
     }
 
+    #[test]
+    fn hemisphere_cover_fails_loud() {
+        // 24 order-1 cells (base cells 0-5) tile exactly half the sphere
+        // (area = 2π), where exterior/hole winding is ambiguous (issue #108).
+        let cover: Vec<u64> = (0..24u64)
+            .map(|nest| crate::morton::nested2mort(nest, 1))
+            .collect();
+        assert!((cover_area(&cover) - std::f64::consts::TAU).abs() < 1e-12);
+        let err = dissolve(&cover, 1)
+            .err()
+            .expect("hemisphere must be rejected");
+        assert!(err.contains("hemisphere"), "{err}");
+    }
+
+    #[test]
+    fn equatorial_band_fails_loud_not_cryptic() {
+        // A thin equatorial band (~1.3 sr, far under the area guard) has two
+        // boundary rings that each enclose more than a hemisphere, so the fan
+        // sum wraps mod 4π: the Σ-vs-exact-area cross-check must reject it
+        // with the curated message, not die in the antimeridian stitcher.
+        let mut band: Vec<u64> = Vec::new();
+        for ilon in 0..720 {
+            for ilat in -7..8 {
+                let lat = ilat as f64 * 0.5;
+                let lon = -180.0 + ilon as f64 * 0.5;
+                band.push(crate::geo2mort::geo2mort_scalar(lat, lon, 4));
+            }
+        }
+        band.sort_unstable();
+        band.dedup();
+        let err = dissolve(&band, 1).err().expect("band must be rejected");
+        assert!(err.contains("hemisphere"), "{err}");
+    }
+
+    #[test]
+    fn sub_hemisphere_cover_still_dissolves() {
+        // Base cells 0-3 (2/3 of a hemisphere) stay outside the guard.
+        let cover: Vec<u64> = (0..16u64)
+            .map(|nest| crate::morton::nested2mort(nest, 1))
+            .collect();
+        let got = dissolve(&cover, 1).unwrap();
+        assert!(!got.shells.is_empty());
+    }
+
     #[test]
     fn pole_cap_stitches_to_one_ring_with_pole_vertex() {
         // one segment running +180 -> -180 around the pole, stitched through -90.

src_rust/src/lib.rs

@@ -3,6 +3,9 @@
 //! This module provides Python bindings for fast morton encoding operations,
 //! replacing the numba-accelerated functions to eliminate Dask conflicts.
 
+// False positives on the pyo3/numpy `?` bridges, where the "useless" error conversion is load-bearing (issue #108).
+#![allow(clippy::useless_conversion)]
+
 pub mod arrow_ffi;
 pub mod buffer;
 pub mod cell_geom;
@@ -1133,14 +1136,17 @@ fn rust_mi_decimal_repr(py: Python<'_>, morton_array: PyReadonlyArray1<u64>) ->
 /// `(lon, lat)` degrees.  Crossing rings are cut at +/-180 and reconnected by
 /// the GeoJSON convention (explicit +/-90 pole vertices stitched down the
 /// antimeridian for a pole-enclosing region).  The Python side builds the
-/// backend Polygons and nests holes — see `mortie/geometry.py`.
+/// backend Polygons and nests holes — see `mortie/geometry.py`.  Raises
+/// `ValueError` for a cover spanning near or over a hemisphere, where the
+/// exterior/hole winding sign is ambiguous (issue #108).
 #[pyfunction]
 #[pyo3(signature = (morton, step=1))]
 fn rust_dissolve(py: Python<'_>, morton: PyReadonlyArray1<u64>, step: u32) -> PyResult<PyObject> {
     let data = morton.to_vec()?;
     let result = py.allow_threads(|| std::panic::catch_unwind(|| dissolve::dissolve(&data, step)));
     let classified = match result {
-        Ok(c) => c,
+        Ok(Ok(c)) => c,
+        Ok(Err(msg)) => return Err(PyValueError::new_err(msg)),
         Err(e) => return Err(PyValueError::new_err(panic_msg(e, "dissolve panicked"))),
     };
     let to_list = |rings: Vec<Vec<(f64, f64)>>| {