Disclose how much of the sphere a mesh covers

CHANGELOG.md

@@ -5,6 +5,31 @@ uses Semantic Versioning for public releases.
 
 ## Unreleased
 ### Fixed
+- Nothing in an area result said how much of the sphere it was summed over. A
+  5-degree mesh spanning 0-40E/0-40N with `sphere_radius: 6371000.0` returned
+  `total_area: 22936016715559.137 m^2` — 4.4967% of `4*pi*R^2` — with
+  `scientific_status: complete`, `physically_interpretable: true`, no warning
+  code, and a bare `postconditions: {status: not_evaluated, checks: []}`; the
+  identical call on a global mesh returned 1.0000 of the sphere with the same
+  status shape. `calculate_area` and `inspect_mesh` now both carry a
+  `mesh_coverage` block: `sphere_fraction`, `closed`, `euler_characteristic`,
+  `lon_extent`, `lat_extent`. A patch raises `MESH_NOT_GLOBAL` and drops to
+  `warning` but stays interpretable, because its total is a real physical
+  quantity and only the missing disclosure was wrong. The geometric and
+  topological halves are reported separately and are allowed to disagree: a
+  1-degree structured global grid stops half a cell short of each pole, so it
+  reads `sphere_fraction: 0.999963` with `closed: false` and 720 boundary
+  edges, which is honest on both counts and not a regional patch. Counting
+  edge incidences is a Python loop — 1.43 s at 196,608 faces, 5.99 s at
+  786,432 — so above 250,000 faces the topological half is skipped and
+  `closed` comes back `null` with `topology_skipped` giving the reason,
+  never `false`.
+- An abstained postcondition now says why it abstained. The area identity
+  holds only on a closed mesh, so a regional result came back
+  `{status: not_evaluated, checks: []}` and the payload never distinguished
+  that from a deployment running `verdict_policy: off`. The block now carries
+  `not_evaluated_because` when the server can name a reason, read off
+  `mesh_coverage` so naming it costs no second traversal of the mesh.
 - The response contract described a payload the server does not send. It
   declared a top-level `physically_interpretable` boolean that no code path
   emits — every producer nests that verdict inside `scientific_status` — and

docs/api.rst

@@ -16,6 +16,10 @@ These modules contain the pure computation logic, separate from MCP and I/O.
    :members:
    :undoc-members:
 
+.. automodule:: uxarray_mcp.domain.mesh_coverage
+   :members:
+   :undoc-members:
+
 .. automodule:: uxarray_mcp.domain.variable
    :members:
    :undoc-members:

docs/tools.md

@@ -312,6 +312,34 @@ Today `calculate_area` is the operation with a closed-form reference: the
 face areas of a closed mesh must sum to `4*pi*R^2`, or `4*pi` on a unit
 sphere. The check abstains — status `not_evaluated` — whenever it cannot be
 trusted: an open or regional mesh, a missing total, or an unreadable grid.
+When the server can name the reason, the block carries it as
+**`not_evaluated_because`**; a bare `not_evaluated` was indistinguishable
+from a deployment running `verdict_policy: off`.
+
+## Mesh coverage
+
+`calculate_area` and `inspect_mesh` both report a **`mesh_coverage`** block:
+`sphere_fraction`, `closed`, `euler_characteristic`, `lon_extent` and
+`lat_extent`. A 5-degree mesh spanning 0-40E/0-40N returns `total_area:
+22936016715559.137 m^2`, which is 4.4967% of `4*pi*R^2`; the same call on a
+global mesh returns 1.0000 of the sphere, and before this block nothing in
+either payload said which was which. A patch is warned about
+(`MESH_NOT_GLOBAL`) rather than refused: its total is a real physical
+quantity, and only the missing disclosure was wrong.
+
+The geometric and topological halves answer different questions and can
+honestly disagree. A 1-degree structured global grid reads `sphere_fraction:
+0.999963` with `closed: false` and `euler_characteristic: 0` — it stops half
+a cell short of each pole, so it covers essentially the whole sphere and is
+genuinely open. It is not warned about, because 3.7e-5 of the sphere is not
+a regional patch, but the area identity still abstains and says why.
+
+Counting edge incidences is a Python loop over every face — 1.43 s on a
+196,608-face HEALPix mesh, 5.99 s at the next zoom level — so above 250,000
+faces the topological half does not run. `closed` and
+`euler_characteristic` come back `null` with `topology_skipped` giving the
+reason, rather than a verdict nobody computed. `sphere_fraction` is
+vectorized and is always reported.
 
 `calculate_area` also declares which sphere it measured on. UXarray computes
 face areas on the unit sphere and never applies `sphere_radius`, so a global

src/uxarray_mcp/domain/__init__.py

@@ -6,6 +6,11 @@
 
 from .area import compute_area_stats
 from .mesh import is_healpix_spec, load_dataset, load_grid, parse_healpix_zoom
+from .mesh_coverage import (
+    compute_mesh_coverage,
+    mesh_coverage_warning_codes,
+    mesh_is_closed,
+)
 from .profile_coverage import (
     compute_profile_coverage,
     profile_coverage_warning_codes,
@@ -30,6 +35,9 @@
     "is_healpix_spec",
     "parse_healpix_zoom",
     "compute_area_stats",
+    "compute_mesh_coverage",
+    "mesh_coverage_warning_codes",
+    "mesh_is_closed",
     "compute_profile_coverage",
     "profile_coverage_warning_codes",
     "compute_target_coverage",

src/uxarray_mcp/domain/area.py

@@ -21,6 +21,8 @@
 
 from typing import Any
 
+from .mesh_coverage import compute_mesh_coverage
+
 #: UXarray's default when a grid declares nothing, and its unit-sphere basis.
 UNIT_SPHERE_RADIUS = 1.0
 
@@ -55,13 +57,18 @@ def compute_area_stats(grid: Any, sphere_radius: float | None = None) -> dict:
     if hasattr(face_areas, "attrs") and "units" in face_areas.attrs:
         area_units = face_areas.attrs["units"]
 
+    steradians = float(face_areas.sum())
     stats = {
-        "total_area": float(face_areas.sum()),
+        "total_area": steradians,
         "mean_area": float(face_areas.mean()),
         "min_area": float(face_areas.min()),
         "max_area": float(face_areas.max()),
         "area_units": area_units,
         "n_face": int(grid.n_face),
+        # Attached before scaling, and measured on the unit sphere whatever
+        # radius is applied below: a total is only readable as global or
+        # regional next to the fraction of the sphere it was summed over.
+        "mesh_coverage": compute_mesh_coverage(grid, steradians=steradians),
     }
     radius, source = resolve_sphere_radius(grid, sphere_radius)
     return apply_sphere_radius(stats, radius, source)

src/uxarray_mcp/domain/mesh_coverage.py

@@ -0,0 +1,242 @@
+"""How much of the sphere a mesh covers, and what shape it is.
+
+``calculate_area`` returns a regional patch's total the same way it returns
+a global one. Measured on a 5-degree mesh spanning 0-40E/0-40N with
+``sphere_radius=6371000.0``: ``total_area`` 22936016715559.137 m^2, which is
+4.4967% of ``4*pi*R^2``, delivered with ``physically_interpretable: True``,
+no warning codes, and ``postconditions: not_evaluated`` that never says why.
+The same call on a global mesh returns 1.0000 of the sphere. Nothing in the
+payload separated them.
+
+Two independent measurements are reported because they answer different
+questions and can honestly disagree:
+
+``sphere_fraction``
+    ``sum(face_areas) / (4*pi)`` on the unit sphere. Geometric: how much
+    surface is actually covered. Reported instead of the raw steradian sum
+    rather than alongside it -- the two differ by a mathematical constant,
+    and this block rides on every area and inspection result under a byte
+    budget (#83).
+``closed`` / ``euler_characteristic``
+    Topological: whether every edge is shared by exactly two faces, and
+    ``V - E + F``.
+
+A 1-degree structured global grid shows why both are needed. It reads
+``sphere_fraction`` 0.999962 and ``closed`` False, ``euler_characteristic``
+0 -- ``Grid.from_structured`` stops its nodes at +/-89.5 there rather than
+extending to the poles, so the mesh has two small polar holes. It covers
+essentially the whole sphere and is genuinely open, and a single verdict
+would have had to suppress one of those facts. The 2-degree grid of the
+same family does reach the poles: 1.000000 and closed, ``euler`` 2. The
+regional patch is a disk: ``euler`` 1.
+
+Cost is why the topology half is size-guarded. Counting edge incidences is
+a Python loop over every face; on a 196,608-face HEALPix mesh it takes
+1.43 s on top of 0.76 s for ``n_edge``, and the next zoom level -- 786,432
+faces -- costs 9.7 s for the pair. ``face_areas`` is vectorized and stays
+under 0.05 s across all of these, so the geometric half is always
+computed and the topological half abstains above the threshold, saying so
+rather than reporting ``closed: false`` for a mesh nobody looked at.
+"""
+
+from __future__ import annotations
+
+import math
+from typing import Any
+
+import numpy as np
+
+#: How far ``sphere_fraction`` may sit from 1.0 and still count as global.
+#:
+#: Quadrature error is three orders of magnitude smaller than this -- a
+#: 162-face global mesh integrates to 1.000002 -- so the slack is not for
+#: numerical noise. It is for meshes that are global in every sense a
+#: caller cares about but leave a pinhole somewhere: the 1-degree grid
+#: above misses 3.8e-5 of the sphere at its poles. A mesh missing more
+#: than 0.1% is missing something a caller would want named.
+GLOBAL_COVERAGE_TOLERANCE = 1e-3
+
+#: Face count above which the topology checks abstain rather than run.
+#:
+#: Set just above HEALPix zoom 7 (196,608 faces, 2.2 s for the pair) and
+#: below zoom 8 (786,432 faces, 9.7 s). An inspection call that takes ten
+#: seconds to report a boolean is not worth the boolean.
+TOPOLOGY_MAX_FACES = 250_000
+
+#: Decimal places used when matching node coordinates. Six is ~0.1 m on
+#: Earth's surface, far below any mesh spacing we deal with, and coarse
+#: enough to absorb the round-trip through NetCDF float64 text.
+_COORD_DECIMALS = 6
+
+#: Decimal places the reported numbers are rounded to. Full float64 repr
+#: costs ~20 characters each on a block that is re-sent every turn, and
+#: buys nothing a caller can use: 1e-6 of the sphere is 510 km^2, and 1e-6
+#: degree is ~0.1 m.
+_REPORT_DECIMALS = 6
+
+
+def _canonical_node_ids(grid: Any) -> list[int]:
+    """Map nodes onto identity by position, not by index.
+
+    A structured global grid stores the 0/360 seam twice and every pole
+    once per meridian, so counting edges on raw indices reports boundary
+    edges on a mesh that is geometrically closed. Merging nodes that sit
+    at the same point -- with all pole nodes collapsing to one, since
+    longitude is meaningless there -- makes the count reflect the surface
+    rather than the storage layout.
+    """
+    lon = np.asarray(grid.node_lon, dtype=float) % 360.0
+    lat = np.asarray(grid.node_lat, dtype=float)
+    seen: dict[str, int] = {}
+    ids: list[int] = []
+    for x, y in zip(lon, lat):
+        if abs(abs(y) - 90.0) < 1e-9:
+            key = f"pole{y:+.1f}"
+        else:
+            key = (
+                f"{round(x, _COORD_DECIMALS) % 360:.6f}_{round(y, _COORD_DECIMALS):.6f}"
+            )
+        ids.append(seen.setdefault(key, len(seen)))
+    return ids
+
+
+def mesh_is_closed(grid: Any) -> bool:
+    """True when every edge is shared by exactly two faces.
+
+    A closed mesh is the precondition for the ``4*pi*R^2`` identity. The
+    cheap version of this test -- comparing ``n_edge`` against Euler's
+    formula -- is wrong on meshes with holes, so count edge incidences
+    directly.
+    """
+    try:
+        connectivity = np.asarray(grid.face_node_connectivity)
+        node_ids = _canonical_node_ids(grid)
+    except Exception:  # pragma: no cover - mocked grids in unit tests
+        return False
+
+    n_node = len(node_ids)
+    incidence: dict[tuple[int, int], int] = {}
+    for face in connectivity:
+        nodes: list[int] = []
+        for raw in face:
+            index = int(raw)
+            if not 0 <= index < n_node:
+                continue  # fill value: a face with fewer nodes than the max
+            node = node_ids[index]
+            if not nodes or nodes[-1] != node:
+                nodes.append(node)
+        # A ring stored with a repeated first/last node is one edge, not two.
+        if len(nodes) > 1 and nodes[0] == nodes[-1]:
+            nodes.pop()
+        if len(nodes) < 3:
+            continue  # degenerate after merging coincident nodes
+        for index, node in enumerate(nodes):
+            other = nodes[(index + 1) % len(nodes)]
+            key = (min(node, other), max(node, other))
+            incidence[key] = incidence.get(key, 0) + 1
+    if not incidence:
+        return False
+    return all(count == 2 for count in incidence.values())
+
+
+def compute_mesh_coverage(
+    grid: Any,
+    *,
+    steradians: float | None = None,
+) -> dict[str, Any]:
+    """Measure how much of the sphere ``grid`` covers and what shape it is.
+
+    Parameters
+    ----------
+    grid : ux.Grid
+        Loaded UXarray grid.
+    steradians : float | None
+        ``sum(face_areas)`` on the unit sphere, when the caller has already
+        computed it. Passed in from ``compute_area_stats`` so the sum is not
+        paid for twice; recomputed here when absent. Not itself reported --
+        it becomes ``sphere_fraction``.
+
+    Returns
+    -------
+    dict
+        ``sphere_fraction`` (geometric), ``closed`` and
+        ``euler_characteristic`` (topological, ``None`` when skipped),
+        ``lon_extent`` and ``lat_extent``. ``topology_skipped`` appears only
+        when the mesh was too large to check, carrying the reason, so a
+        ``None`` verdict is never mistaken for a negative one.
+
+        ``lon_extent`` describes the mesh in the grid's own longitude
+        convention and is not a globality test: a global mesh stored on
+        [-180, 180] with 20-degree cells reads [-170, 170]. Use
+        ``sphere_fraction`` for that.
+    """
+    coverage: dict[str, Any] = {
+        "sphere_fraction": None,
+        "closed": None,
+        "euler_characteristic": None,
+        "lon_extent": None,
+        "lat_extent": None,
+    }
+
+    if steradians is None:
+        try:
+            steradians = float(np.asarray(grid.face_areas).sum())
+        except Exception:  # pragma: no cover - mocked grids in unit tests
+            steradians = None
+    if steradians is not None and math.isfinite(steradians):
+        coverage["sphere_fraction"] = round(
+            float(steradians) / (4.0 * math.pi), _REPORT_DECIMALS
+        )
+
+    try:
+        lon = np.asarray(grid.node_lon, dtype=float)
+        lat = np.asarray(grid.node_lat, dtype=float)
+        if lon.size and lat.size:
+            coverage["lon_extent"] = [
+                round(float(lon.min()), _REPORT_DECIMALS),
+                round(float(lon.max()), _REPORT_DECIMALS),
+            ]
+            coverage["lat_extent"] = [
+                round(float(lat.min()), _REPORT_DECIMALS),
+                round(float(lat.max()), _REPORT_DECIMALS),
+            ]
+    except Exception:  # pragma: no cover - mocked grids in unit tests
+        pass
+
+    try:
+        n_face = int(grid.n_face)
+    except Exception:  # pragma: no cover - mocked grids in unit tests
+        return coverage
+
+    if n_face > TOPOLOGY_MAX_FACES:
+        coverage["topology_skipped"] = (
+            f"{n_face} faces exceeds the {TOPOLOGY_MAX_FACES}-face limit for "
+            "counting edge incidences; closure was not checked."
+        )
+        return coverage
+
+    try:
+        coverage["euler_characteristic"] = int(grid.n_node) - int(grid.n_edge) + n_face
+    except Exception:  # pragma: no cover - mocked grids in unit tests
+        pass
+    coverage["closed"] = mesh_is_closed(grid)
+    return coverage
+
+
+def mesh_coverage_warning_codes(coverage: dict[str, Any]) -> list[str]:
+    """Stable codes for a mesh that is not the whole sphere.
+
+    Silent on a skipped topology check: not knowing whether a mesh is
+    closed is not evidence that it is open, and ``topology_skipped`` in
+    the block already says nobody looked.
+    """
+    fraction = coverage.get("sphere_fraction")
+    if fraction is None:
+        return []
+    if fraction < 1.0 - GLOBAL_COVERAGE_TOLERANCE:
+        return ["MESH_NOT_GLOBAL"]
+    if fraction > 1.0 + GLOBAL_COVERAGE_TOLERANCE:
+        # More surface than a sphere has means faces overlap or are stored
+        # twice. Quadrature cannot produce this at 1e-3.
+        return ["MESH_COVERAGE_EXCEEDS_SPHERE"]
+    return []