Drop the seam hand-off corner in section-finding, not in transports

CLAUDE.md

@@ -122,9 +122,9 @@ The package is organized around a pipeline: define sections → map to grid →
 
 ### Core Modules
 
-- **`section.py`** — Section definition and grid path algorithms. `Section` holds named waypoint coordinates; `GriddedSection` extends it with grid index information. `grid_section()` is the main entry point that maps geographic waypoints to grid vorticity-point indices `(i_c, j_c)` by walking the grid between waypoints along the `curve` requested (`"great circle"`, the default geodesic, or `"latitude circle"`); each segment must span less than 180°. The walk is deterministic and direction-independent: it admits neighbors strictly closer to the endpoint plus any seam twin of the current cell (so periodic/fold-seam crossings do not depend on floating-point rounding), and breaks ties by index. Grid topology is inferred entirely from `xgcm.Grid` metadata — each axis' `boundary` (periodic wrap, fill/extend wall, or a single-tile bipolar north fold `{"Y": {"fold": ...}}`) and `face_connections` for multi-tile grids — so there is no `topology` keyword. The pathfinder consumes topology-aware neighbor maps built by `gridutils.build_neighbor_maps`.
+- **`section.py`** — Section definition and grid path algorithms. `Section` holds named waypoint coordinates; `GriddedSection` extends it with grid index information. `grid_section()` is the main entry point that maps geographic waypoints to grid vorticity-point indices `(i_c, j_c)` by walking the grid between waypoints along the `curve` requested (`"great circle"`, the default geodesic, or `"latitude circle"`); each segment must span less than 180°. The walk is deterministic and direction-independent: it admits neighbors strictly closer to the endpoint plus any seam twin of the current cell (so periodic/fold-seam crossings do not depend on floating-point rounding), and breaks ties by index. That seam-twin step visits the same physical corner twice, so it spans no grid cell and is not a velocity face; **`grid_section()` removes the redundant identity right after the walk, via `drop_repeated_corners()`**, keeping the later one (the frame the path continues in). Every section `grid_section()` returns therefore satisfies the invariant that **no two consecutive corners are the same physical point** — every consecutive pair is a real velocity face — and, on a multi-tile grid, that **every corner is the canonical native index of its corner node** (`_OuterTopology.node_native`). That is the only place the normalization happens: `transports.uvindices_from_qindices` *checks* the invariant and raises rather than re-establishing it, so hand-built index arrays are the caller's to fix (by calling `drop_repeated_corners` themselves). `_check_supported_topology()` additionally rejects, at the front door, a multi-tile grid that registers no cell-center dimension on both horizontal axes: the corner-node topology is rebuilt by fingerprinting the tracer *cells* around each corner, and velocities are staggered onto the center dimensions (U at (X-corner, Y-center), V at (X-center, Y-corner)), so such a grid can carry no velocities and its seam corners cannot be resolved to one identity. Only the center **dimensions** are needed — `_OuterTopology` numbers cells with `np.arange`, so no tracer longitude/latitude values are required. Single-tile grids are unaffected: their seam twins are detected geometrically from the corner coordinates alone, so a corner-only single-tile grid still traces (it just cannot produce transports, which `gridutils.coord_dict` reports by name). Grid topology is inferred entirely from `xgcm.Grid` metadata — each axis' `boundary` (periodic wrap, fill/extend wall, or a single-tile bipolar north fold `{"Y": {"fold": ...}}`) and `face_connections` for multi-tile grids — so there is no `topology` keyword. The pathfinder consumes topology-aware neighbor maps built by `gridutils.build_neighbor_maps`.
 
-- **`transports.py`** — Transport computation along sections. `uvindices_from_qindices()` converts vorticity-point indices to U/V velocity-point indices using a per-position corner offset (`gridutils.corner_offset`) covering all three C-grid staggerings: 'outer', 'right', and 'left' (see "Corner staggering" below). `convergent_transport()` is the main function: it lazily computes signed normal transports with configurable orientation (positive inward to the polygon defined by the section).
+- **`transports.py`** — Transport computation along sections. `uvindices_from_qindices()` converts vorticity-point indices to U/V velocity-point indices using a per-position corner offset (`gridutils.corner_offset`) covering all three C-grid staggerings: 'outer', 'right', and 'left' (see "Corner staggering" below). It emits exactly one face per consecutive pair of corners: it relies on the section-finding invariant above rather than filtering degenerate pairs itself, so an edge whose normal velocity is stored on neither the source nor the destination face raises instead of silently returning a zero-flux placeholder. `convergent_transport()` is the main function: it lazily computes signed normal transports with configurable orientation (positive inward to the polygon defined by the section).
 
 - **`tracers.py`** — `extract_tracer()` interpolates tracer data to U/V points along a section path for cross-section plotting.
 

sectionate/gridutils.py

@@ -142,9 +142,28 @@ def coord_dict(grid):
     dict
         Dictionary containing names of "X" and "Y" dimension variables, at both cell 'center'
         position and the corner position ('outer', 'right', or 'left'; see `corner_position`).
+
+    Raises
+    ------
+    ValueError
+        If either horizontal axis registers no 'center' position. Velocities are staggered
+        onto the center dimensions -- U lives at (X-corner, Y-center) and V at (X-center,
+        Y-corner) -- so a grid without them cannot locate any velocity, and anything that
+        derives faces, transports or tracers from a section needs this mapping.
     """
     corner_pos = corner_position(grid)
 
+    missing = [ax for ax in ("X", "Y") if "center" not in grid.axes[ax].coords]
+    if missing:
+        raise ValueError(
+            f"The {' and '.join(missing)} axis of this grid registers no 'center' "
+            "position, so sectionate cannot locate its velocities: U is stored at "
+            "(X-corner, Y-center) and V at (X-center, Y-corner). Declare the cell-center "
+            "dimensions, e.g. coords={'X': {'center': 'xh', 'outer': 'xq'}, "
+            "'Y': {'center': 'yh', 'outer': 'yq'}}. (Tracing a section with "
+            "`grid_section` needs only the corner coordinates and works without them.)"
+        )
+
     return {
         "X": {
             "center": grid.axes["X"].coords["center"],
@@ -256,15 +275,20 @@ def build_neighbor_maps(grid, geocorners):
         has_centers = all("center" in grid.axes[ax].coords for ax in ("X", "Y"))
         if has_centers:
             return outer_topology(grid).maps
-        # Without tracer-center coordinates the cell-identity construction is
-        # unavailable (and neither are velocities, so only walking is needed):
-        # fall back to reading xgcm's corner halos directly. Only shared-corner
-        # ('outer') tilings are safe here -- staggered corner arrays can pad one
-        # corner off across rotated/reversed seams.
+        # The cell-identity construction fingerprints each corner by the tracer cells
+        # around it, so it needs the cell-center *dimensions* (their names and lengths;
+        # no tracer longitude/latitude values are read). Without them, fall back to
+        # reading xgcm's corner halos directly. These maps step through each face's own
+        # copy of a shared seam corner, so a section traced on them would keep the
+        # zero-length hand-off step; `section.grid_section` therefore refuses a
+        # multi-tile grid with no center dimensions, and this fallback only ever serves
+        # a direct `build_neighbor_maps` call. Only shared-corner ('outer') tilings are
+        # safe here -- staggered corner arrays can pad one corner off across
+        # rotated/reversed seams.
         if corner_position(grid) != "outer":
             raise ValueError(
                 "Multi-tile grids with staggered ('left'/'right') corners require "
-                "tracer-center coordinates to derive their corner topology."
+                "cell-center dimensions to derive their corner topology."
             )
         return _multitile_padded_maps(grid, geocorners)
     da = geocorners["X"]
@@ -313,7 +337,7 @@ def pad(a):
 def _multitile_padded_maps(grid, geocorners):
     """
     Fallback multi-tile neighbor maps for shared-corner ('outer') tilings that
-    carry no tracer-center coordinates: pad index-valued corner arrays with
+    register no cell-center dimensions: pad index-valued corner arrays with
     xgcm's `face_connections` halos and read the neighbors off directly. A
     shared seam corner is stored on every face that touches it, so these maps
     step through each face's own copy (seam-twin semantics); `_validate_reciprocity`

sectionate/section.py

@@ -5,6 +5,7 @@
     get_geo_corners,
     get_facedim,
     build_neighbor_maps,
+    outer_topology,
 )
 
 # Two corner indices map to the same physical point on seams that fold or wrap (e.g.
@@ -284,18 +285,121 @@ def grid_section(grid, lons, lats, curve="great circle"):
     # grid's xgcm metadata by padding index arrays (see `build_neighbor_maps`).
     # Where a single physical corner carries two indices (a periodic seam, a shared
     # multi-tile boundary corner, or the fold seam), the walk simply steps through
-    # both; the resulting zero-length section edge carries no flux and is dropped when
-    # faces are derived (see `transports.uvindices_from_qindices`).
+    # both; that hand-off step spans no grid cell, so the extra identity is dropped
+    # right here, by `drop_repeated_corners`, before the section is returned.
     neighbor_maps = build_neighbor_maps(grid, geocorners)
 
-    return create_section_composite(
+    out = create_section_composite(
         geocorners["X"],
         geocorners["Y"],
         lons,
         lats,
         neighbor_maps=neighbor_maps,
         curve=curve,
     )
+    if facedim is not None:
+        i_c, j_c, f_c = out[0], out[1], out[2]
+    else:
+        i_c, j_c, f_c = out[0], out[1], None
+
+    i_c, j_c, f_c, lons_c, lats_c = drop_repeated_corners(grid, i_c, j_c, f_c)
+    if f_c is not None:
+        return i_c, j_c, f_c, lons_c, lats_c
+    return i_c, j_c, lons_c, lats_c
+
+
+def drop_repeated_corners(grid, i_c, j_c, f_c=None):
+    """
+    Normalise a traced grid path so that consecutive corners are always distinct
+    physical points.
+
+    A single physical vorticity point can be stored under more than one index: the
+    duplicated seam column (or row) of a symmetric ('outer') periodic grid, a corner
+    shared by several tiles of a multi-tile grid, or the two coincident lips of a grid
+    cut. `infer_grid_path_from_geo` deliberately admits the "seam twin" of the corner
+    it is standing on, so that a seam crossing is decided by the grid's topology rather
+    than by floating-point rounding. The resulting hand-off step spans no grid cell, so
+    that pair of corners does not define a velocity face.
+
+    `grid_section` calls this immediately after the path is traced, so that every section
+    it returns satisfies:
+
+    * no two consecutive corners are the same physical point -- every consecutive pair
+      is a real velocity face; and
+    * on a multi-tile grid, every corner is the canonical native index of the physical
+      corner point it denotes (`_OuterTopology.node_native`), which is what lets the
+      velocity-face attribution look a corner up in either the source or the
+      destination face's frame.
+
+    Of the two identities of a seam point the LATER one is kept -- the frame the path
+    continues in. The faces on either side of the hand-off then keep exactly the
+    indices they would have had if the walk had never changed frame: the face after
+    the seam is untouched, and the face before it is re-expressed as a wrapped step,
+    which the index arithmetic in `transports.uvindices_from_qindices` already handles.
+
+    `uvindices_from_qindices` requires this invariant and checks it, so hand-built index
+    arrays that were not produced by `grid_section` can be passed through here first.
+
+    PARAMETERS:
+    -----------
+    grid: xgcm.Grid
+        Grid the path was traced on.
+    i_c, j_c: array-like of int
+        Vorticity-point indices along the "X" and "Y" dimensions.
+    f_c: array-like of int or None
+        Face indices for multi-tile grids (`face_connections`); None for single-tile.
+
+    RETURNS:
+    --------
+    i_c, j_c, f_c, lons_c, lats_c: np.ndarray
+        The normalised path and the coordinates of its corners (`f_c` is None for
+        single-tile grids).
+    """
+    i_c = np.asarray(i_c, dtype=np.int64)
+    j_c = np.asarray(j_c, dtype=np.int64)
+    f_c = None if f_c is None else np.asarray(f_c, dtype=np.int64)
+
+    geocorners = get_geo_corners(grid)
+    glon = np.asarray(geocorners["X"].values)
+    glat = np.asarray(geocorners["Y"].values)
+
+    if f_c is not None:
+        # Multi-tile: identity is exact and purely topological -- two indices denote
+        # the same point iff they resolve to the same node of the corner graph.
+        ot = outer_topology(grid)
+        nodes = ot.node_id[f_c, j_c + ot.t, i_c + ot.t]
+        if (nodes < 0).any():
+            raise ValueError(
+                "Section contains indices that are not grid corners."
+            )
+        native = ot.node_native[nodes]
+        if (native[:, 0] < 0).any():
+            raise ValueError(
+                "Section passes through a corner point that is stored on no face of "
+                "the grid (e.g. a cube vertex or an unstored cap/cut corner), so it "
+                "has no native index."
+            )
+        f_c, j_c, i_c = native[:, 0].copy(), native[:, 1].copy(), native[:, 2].copy()
+        repeated = nodes[1:] == nodes[:-1]
+    else:
+        # Single tile: the seam twins of a periodic/fold boundary are separate indices
+        # with no shared identity to compare, so coincidence is measured physically.
+        # The tolerance sits far below any real grid spacing (see COINCIDENT_TOLERANCE_M).
+        lo, la = glon[j_c, i_c], glat[j_c, i_c]
+        repeated = distance_on_unit_sphere(
+            lo[:-1], la[:-1], lo[1:], la[1:]
+        ) < COINCIDENT_TOLERANCE_M
+
+    keep = np.ones(i_c.size, dtype=bool)
+    keep[:-1] = ~repeated                    # keep the LAST index of each run
+    i_c, j_c = i_c[keep], j_c[keep]
+    f_c = None if f_c is None else f_c[keep]
+
+    if f_c is not None:
+        lons_c, lats_c = glon[f_c, j_c, i_c], glat[f_c, j_c, i_c]
+    else:
+        lons_c, lats_c = glon[j_c, i_c], glat[j_c, i_c]
+    return i_c, j_c, f_c, np.asarray(lons_c), np.asarray(lats_c)
 
 
 def _check_supported_topology(grid):
@@ -306,6 +410,16 @@ def _check_supported_topology(grid):
     connects to itself along the "Y" axis) is not supported. Use the single-tile bipolar-fold
     padding instead -- ``padding={"X": "periodic", "Y": {"fold": "corner"}}`` (xgcm >= the
     bipolar-fold release) -- which sectionate handles natively via xgcm's fold padding.
+
+    A multi-tile grid must also register a cell-center dimension on both horizontal axes.
+    That is what the corner-node topology is reconstructed from (`gridutils.outer_topology`
+    fingerprints each corner by the tracer cells around it, because cells -- unlike corner
+    arrays -- pad reliably across any seam), and it is also where half of each velocity's
+    dimensions live: U sits at (X-corner, Y-center) and V at (X-center, Y-corner). Only the
+    *dimensions* are needed, not tracer longitudes/latitudes. Without them a multi-tile grid
+    can carry no velocities and its seam corners cannot be resolved to a single identity, so
+    a section across a seam could not be traced correctly; refuse it at the front door rather
+    than return a path whose seam steps are wrong.
     """
     facedim = grid._facedim
     connections = (getattr(grid, "_face_connections", None) or {}).get(facedim, {})
@@ -323,6 +437,20 @@ def _check_supported_topology(grid):
                         "padding instead: padding={'X':'periodic','Y':{'fold':'corner'}}."
                     )
 
+    missing = [ax for ax in ("X", "Y") if "center" not in grid.axes[ax].coords]
+    if missing:
+        raise ValueError(
+            "Multi-tile grids (`face_connections`) must register a cell-center dimension on "
+            f"both horizontal axes; {' and '.join(missing)} declare none. Sectionate rebuilds "
+            "the shared-corner topology from the tracer cells around each corner, and the "
+            "velocities themselves are staggered onto the center dimensions (U at (X-corner, "
+            "Y-center), V at (X-center, Y-corner)), so a grid without them carries no "
+            "velocities and its seam corners cannot be identified. Declare them, e.g. "
+            "coords={'X': {'center': 'xh', 'outer': 'xq'}, "
+            "'Y': {'center': 'yh', 'outer': 'yq'}} -- the center entries need only name "
+            "dimensions of the dataset, they need not carry coordinate values."
+        )
+
 def create_section_composite(
     gridlon,
     gridlat,
@@ -359,6 +487,12 @@ def create_section_composite(
         (i_c, j_c[, f_c]) correspond to indices of vorticity points that define velocity faces;
         the face index f_c is only returned for multi-tile grids.
         (lons_c, lats_c) are the corresponding longitude and latitudes.
+
+        This is the raw walk: where a section crosses a seam whose corner is stored under
+        two indices, both are still present, and the step between them is not a velocity
+        face. `grid_section` (which has the `xgcm.Grid` this needs) drops the redundant
+        one via `drop_repeated_corners`; call that yourself if you use this routine
+        directly and intend to derive faces from the result.
     """
 
     # A face dimension (3-D corner arrays) marks a multi-tile grid, whose sections

sectionate/tests/test_ecco_llc90.py

@@ -83,6 +83,37 @@ def test_llc90_polar_cut_lips_share_nodes():
     assert shared, "no node shares native storage across tiles 0 and 3"
 
 
+def test_llc90_sections_contain_no_repeated_corner():
+    """Sections traced on the real LLC90 grid satisfy the public invariant: consecutive
+    corners are always distinct physical points -- so every consecutive pair is a real
+    velocity face -- and every corner is the canonical native index of its corner node,
+    which is what lets the face attribution read a corner in either face's frame."""
+    from sectionate.gridutils import outer_topology
+    from sectionate.section import grid_section
+    from sectionate.transports import uvindices_from_qindices
+
+    grid = _load_grid()
+    ot = outer_topology(grid)
+
+    for lat in (-70., 10., 70.):
+        i_c, j_c, f_c, _, _ = grid_section(
+            grid, [0., 90., 180., 270., 360.], [lat] * 5, curve="latitude circle"
+        )
+        nodes = ot.node_id[f_c, j_c + ot.t, i_c + ot.t]
+        assert (nodes >= 0).all()
+        assert np.all(nodes[1:] != nodes[:-1]), f"repeated corner at lat {lat}"
+
+        native = ot.node_native[nodes]
+        assert np.array_equal(native[:, 0], f_c)
+        assert np.array_equal(native[:, 1], j_c)
+        assert np.array_equal(native[:, 2], i_c)
+
+        # every consecutive pair yields a face, and all of them are real velocities
+        uv = uvindices_from_qindices(grid, i_c, j_c, f_c=f_c)
+        assert uv["var"].size == i_c.size - 1
+        assert set(np.unique(uv["var"]).tolist()) <= {"U", "V"}
+
+
 def test_llc90_face_corners_resolve_topologically():
     """Every native corner appears in the neighbour maps with the correct arity:
     building the maps requires resolving the four-tile-junction face corners via