One boundary back-end, on the corner graph

.gitignore

@@ -4,3 +4,7 @@ __pycache__
 
 # Version file generated from the git tag by hatch-vcs at build time.
 regionate/_version.py
+
+# Downloaded example data (large; fetched by examples/load_example_*.py)
+data/*.nc
+data/*.zarr

CLAUDE.md

@@ -4,7 +4,7 @@ This file provides guidance to Claude Code (claude.ai/code) when working with co
 
 ## What this is
 
-`regionate` builds **xgcm-grid-consistent** regional masks and boundaries for ocean/climate model output. It supports **arbitrary `xgcm.Grid` topologies** — single-tile periodic and bipolar-fold (Arctic) MOM6 grids, and genuinely multi-tile grids defined by `face_connections` (e.g. ECCOv4r4 lat-lon-cap, cubed-sphere) — by driving all grid logic from the topology-aware `sectionate` API rather than hard-coded MOM6 specifics. Given a geographic polygon, it snaps the polygon to a discrete model grid, producing a boolean cell mask plus the staggered (u,v) velocity faces that trace the region's boundary — so that volume/mass/heat budgets integrated over the masked region are exactly consistent with fluxes through the boundary faces. It leans heavily on its sibling package [`sectionate`](https://github.com/MOM6-community/sectionate) for the section/face-tracing math (requires the topology-driven API on [`hdrake/sectionate@topology-driven-neighbors`](https://github.com/hdrake/sectionate/tree/topology-driven-neighbors), formerly PRs #47/#48), and on the bipolar north-fold boundary (`padding={..., "Y": {"fold": ...}}`, formerly xgcm#711) plus the multi-tile `face_connections` padding fix (xgcm#712) and the vector-pad fix (xgcm#749), all released in `xgcm >= 0.10.1` (on PyPI).
+`regionate` builds **xgcm-grid-consistent** regional masks and boundaries for ocean/climate model output. It supports **arbitrary `xgcm.Grid` topologies** — single-tile periodic and bipolar-fold (Arctic) MOM6 grids, and genuinely multi-tile grids defined by `face_connections` (e.g. ECCOv4r4 lat-lon-cap, cubed-sphere) — by driving all grid logic from the topology-aware `sectionate` API rather than hard-coded MOM6 specifics. Given a geographic polygon, it snaps the polygon to a discrete model grid, producing a boolean cell mask plus the staggered (u,v) velocity faces that trace the region's boundary — so that volume/mass/heat budgets integrated over the masked region are exactly consistent with fluxes through the boundary faces. It leans heavily on its sibling package [`sectionate`](https://github.com/MOM6-community/sectionate) for the section/face-tracing math (requires the topology-driven API on [`hdrake/sectionate@corner-node-topology`](https://github.com/hdrake/sectionate/tree/corner-node-topology), which is PR #47 plus the corner-topology rewrite), and on the bipolar north-fold boundary (`padding={..., "Y": {"fold": ...}}`, formerly xgcm#711) plus the multi-tile `face_connections` padding fix (xgcm#712) and the vector-pad fix (xgcm#749), all released in `xgcm >= 0.10.1` (on PyPI).
 
 ## AI Usage Policy
 
@@ -34,7 +34,7 @@ pytest -k <name>                                # run a single test by name
 pip install -e .                                # editable install (after creating env)
 ```
 
-Dev environment is conda-based (see README). CI (`.github/workflows/ci.yml`) installs `ci/environment.yml`, installs `sectionate` from `hdrake/sectionate@topology-driven-neighbors` (still unreleased), does `pip install -e .` (which pulls `xgcm >= 0.10.1` straight from PyPI — 0.10.1 ships the north-fold boundary and the multi-tile padding fixes), then runs `pytest` across Python 3.11–3.14. There is no linter configured. Optional end-to-end tests against real MOM6/ECCO grids are skipped unless `REGIONATE_REALDATA_TESTS=1`.
+Dev environment is conda-based (see README). CI (`.github/workflows/ci.yml`) installs `ci/environment.yml`, installs `sectionate` from `hdrake/sectionate@corner-node-topology` (still unreleased), does `pip install -e .` (which pulls `xgcm >= 0.10.1` straight from PyPI — 0.10.1 ships the north-fold boundary and the multi-tile padding fixes), then runs `pytest` across Python 3.11–3.14. There is no linter configured. Optional end-to-end tests against real MOM6/ECCO grids are skipped unless `REGIONATE_REALDATA_TESTS=1`.
 
 ## Versioning
 
@@ -60,7 +60,7 @@ Supporting modules:
 - `regions.py` — **`Regions`**, a dict-like collection mapping names → `Region`/`GriddedRegion`; `overlaps.py` handles intersections between members.
 - `grid_conform.py` — the core grid-snapping logic (polygon → grid indices → mask). `mask_from_grid_boundaries` rasterizes a boundary by splitting it at the antimeridian into a `[-180,180]` MultiPolygon and OR-ing per-piece `regionmask` masks (pole-encircling boundaries are extended to the South Pole first).
 - `geometry.py` — antimeridian / MultiPolygon shapely helpers (`split_at_antimeridian`, `normalize_lon`, ...) used by `grid_conform.py`.
-- `boundaries.py` — `grid_boundaries_from_mask` (inverse: recover boundary faces from a mask). A shared front-end `_trace_and_drop` traces the mask per face with `contourpy` and drops boundary segments interior to a seam (both separated cells in-mask in the topology-aware halo `_pad_center`); two back-ends then stitch the surviving arcs: `_single_tile_boundaries_from_mask` stitches by physical seam coincidence (handles walls, **periodic** axes, and the bipolar **fold** uniformly; `f_c=None`), and `_multitile_boundaries_from_mask` stitches on the grid's outer (shared-corner) corner-node graph from `sectionate.gridutils.outer_topology`: every traced corner resolves to a physical corner node (uniform across rotated/reversed seams, cube-vertex junctions, the pole, and grid cuts/folds), edges traced from both sides of an undeclared cut/fold annihilate as interior, surviving segments chain into closed loops by node identity, and each node is emitted as the native corner that stores it (`f_c_list`; `None` for single-tile). Boundaries through points stored on no face (LLC90's 4th Arctic-cap vertex, single-sided contact with its lon=-115 Antarctic cut) raise rather than fabricate indices. A region wrapping any seam yields one seam-consistent boundary (an annulus/strip may yield several loops); `_pad_center` pads only genuine seams (periodic/fold/face_connections) with real neighbours and everything else with NaN, so walls are never mistaken for seams. For the budget's *area side* on multi-tile staggered grids, use `outer_topology(grid).padded_transports(u, v)` — exact and xgcm-independent, and the only path that handles edges stored on no face (walls, grid cuts, a cap's un-stored vertex). (xgcm#749 fixes the bare-`DataArray` `grid.diff(..., other_component=)` path across rotated/reversed seams — the dict form was always exact — but even a correct pad cannot supply a halo for edges stored on no face, so `padded_transports` remains preferred here.)
+- `boundaries.py` — `grid_boundaries_from_mask` (inverse: recover boundary faces from a mask). `_trace_and_drop` traces the mask per face with `contourpy` and drops boundary segments interior to a seam (both separated cells in-mask in the topology-aware halo `_pad_center`). **One** back-end, `_boundaries_from_arcs`, then stitches the surviving arcs on the grid's corner topology (`sectionate.topology.corner_topology`): every traced corner resolves to a physical corner *node* — uniform across a periodic wrap, a bipolar fold, rotated and reversed tile seams and cube-vertex junctions alike — edges traced from both sides of a seam annihilate as interior, surviving segments chain into closed loops by node identity, and each node is emitted as the native corner that stores it, with its face index. `f_c` is present on every grid; a single-tile grid is one face, so it is zeros there. The corner frame is the 'outer' lattice (`o = 1`) whatever the native staggering, which is what the corner topology is indexed on — it used to depend on the staggering, and a `'right'` grid was traced one corner off. Boundaries through corners stored on no face raise rather than fabricate indices. `_pad_center` pads only genuine seams (periodic/fold/`face_connections`) with real neighbours and everything else with NaN, so walls are never mistaken for seams — but note it reads the grid's *metadata* only, so it does not see identifications declared through `sectionate.topology.declare_identifications` (LLC90's southern boundary fold); `connected_components` therefore treats such a seam as unconnected. On LLC90 those cells are Antarctic land, so it is currently unobservable, but it is a region-identity bug waiting for a declared seam that carries flow. For the budget's *area side* on multi-tile staggered grids, use `corner_topology(grid).padded_transports(u, v)` — exact, and the only path that handles edges stored on no face.
 - `integrate.py` — `check_global_coverage` validates that a `Regions` set is non-overlapping and tiles the globe.
 - `utilities.py` — shared helpers (re-exported widely via `from .utilities import *`).
 

examples/load_example_ECCO_grid.py

@@ -82,6 +82,47 @@ def download_ECCO_geometry(data_dir="../data"):
     return download_ECCO_data(ECCO_GEOMETRY_FILE, data_dir=data_dir)
 
 
+# The LLC90 domain's southern edge folds onto itself, along the great circle
+# through 65E and 115W under Antarctica. `face_connections` declares those four
+# tile edges as `None` -- walls -- because its schema maps a whole tile edge to a
+# whole tile edge, and this edge is glued to several partners over different index
+# ranges: tile 0's runs into tile 3's, tile 9's and tile 12's in turn, and two of
+# them fold onto themselves at the pivots. So one line of physical corners is
+# stored twice with nothing in the metadata to say the two copies are the same
+# point, and a region straddling it cannot close.
+#
+# Each block below is an affine map on the 'outer' corner lattice, `k -> const - k`
+# along the named edge. Together they merge 179 corner pairs; every merged pair is
+# bit-identically coincident in the archived coordinates (`validate_positions()`
+# returns 0.0 m), which is a check on the declaration rather than its source.
+LLC90_SOUTHERN_FOLD_BLOCKS = [
+    #  (face, edge,     range,        face, edge,     const)
+    (0,  "Y_low",  range(91),      3,  "Y_low",  116),
+    (3,  "Y_low",  range(27),      3,  "Y_low",   26),
+    (0,  "Y_low",  range(27),      9,  "X_high",  26),
+    (9,  "X_high", range(91),      12, "X_high", 116),
+    (12, "X_high", range(27),      12, "X_high",  26),
+]
+
+
+def llc90_southern_fold(n=91):
+    """The LLC90 southern boundary fold, as explicit corner identifications."""
+    from sectionate.topology import Identification
+
+    def slot(f, edge, k):
+        return (f, 0, k) if edge == "Y_low" else (f, k, n - 1)
+
+    out = []
+    for fa, ea, ks, fb, eb, const in LLC90_SOUTHERN_FOLD_BLOCKS:
+        ks = [k for k in ks if 0 <= const - k < n]
+        out.append(Identification(
+            [slot(fa, ea, k) for k in ks],
+            [slot(fb, eb, const - k) for k in ks],
+            name=f"LLC90 southern fold {fa}.{ea} <-> {fb}.{eb}",
+        ))
+    return out
+
+
 def load_ECCO_LLC90_grid(data_dir="../data"):
     """Load the ECCOv4r4 LLC90 grid as a native (MITgcm 'left'-staggered)
     ``xgcm.Grid``.
@@ -97,12 +138,16 @@ def load_ECCO_LLC90_grid(data_dir="../data"):
     ds = xr.open_dataset(path)
     ds = ds.rename({"XC": "geolon", "YC": "geolat",
                     "XG": "geolon_c", "YG": "geolat_c"})
-    return xgcm.Grid(
+    from sectionate.topology import declare_identifications
+
+    grid = xgcm.Grid(
         ds, padding="fill", autoparse_metadata=False,
         coords={"X": {"center": "i", "left": "i_g"},
                 "Y": {"center": "j", "left": "j_g"}},
         face_connections=LLC90_FACE_CONNECTIONS,
     )
+    declare_identifications(grid, llc90_southern_fold())
+    return grid
 
 
 # Natural Earth ocean basins that make up the Atlantic sector (open Atlantic plus

pyproject.toml

@@ -19,7 +19,10 @@ dependencies = [
     "xgcm >= 0.10.1",
     "regionmask >= 0.13",
     "contourpy",
-    "sectionate >= 0.4.0rc2"
+    # Needs the corner-topology API (`sectionate.topology`), which is not in any
+    # release yet: it is the PR this one is built on. Until it ships, install
+    # sectionate from that branch; the floor below is what to require once it does.
+    "sectionate >= 0.5.0"
 ]
 
 [project.urls]