Accept either seam representation, and map velocity faces back to corner steps

CLAUDE.md

@@ -139,7 +139,7 @@ The package is organized around a pipeline: define sections → map to grid →
 
 - **`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. There are three: `"great circle"` (the default — every segment follows the geodesic), `"latitude circle"` (every segment follows a parallel, and a segment whose endpoints do not share a latitude — meridional ones included — raises), and `"latitude and great circle"` (decided per segment: constant-latitude segments follow the parallel, all others the geodesic). Under all three, every segment follows the **shortest** path between its two waypoints — raw longitudes are never read as a request to go the long way round, so a `"latitude circle"` segment written `0 -> 270` runs 90° *west* — which is why encircling the globe takes intermediate waypoints (e.g. `0 -> 120 -> 240 -> 360`). `_check_segment_span` resolves each segment's curve (via `_is_constant_latitude`, the one classification the metric choice and the legality check share) and raises for ill-posed ones: endpoints written exactly half a circle apart (neither way round is shorter — one error shared by all curves, quoting degrees of longitude along the parallel or degrees of arc, as applicable), or, under `"latitude circle"`, endpoints that do not share a latitude. 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`.
 
-- **`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). Before doing so it resolves **corner identity**, since a seam gives one physical corner two indices: multi-tile grids canonicalize each corner through `_OuterTopology`, and single-tile grids run `_insert_seam_twins()`, which splices the seam twin into the corner list wherever a step skips it (see "Seam twins" 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).
 
 - **`tracers.py`** — `extract_tracer()` interpolates tracer data to U/V points along a section path for cross-section plotting.
 
@@ -151,6 +151,8 @@ The package is organized around a pipeline: define sections → map to grid →
 
 - **Vorticity points (q-points):** Sections are defined as paths through vorticity-point indices `(i_c, j_c)`. Consecutive q-points define velocity faces (either U or V).
 - **Corner staggering (three positions):** Vorticity sits at one of three xgcm positions: `"outer"` (MOM6 symmetric, M+1×N+1), `"right"` (MOM6 non-symmetric, M×N), or `"left"` (MITgcm/ECCO, incl. the lat-lon-cap, M×N). All three are native; they differ only by a per-position velocity index offset (`gridutils.corner_offset`: outer→0, right→+1, left→0). `"left"` indexes like `"outer"`; it differs only in array length and in that the *high* corner row/column is absent (so a section exactly on the north/east domain wall clips one row inside).
+- **Seam twins (corner identity):** A seam gives one physical corner more than one index — the periodic wrap's first/last column, a shared multi-tile boundary corner, and the bipolar fold's mirrored seam-row columns. Transport attribution therefore has to know *which* indices denote the same point, or it will name the mirror image of the face a section crossed (with the sign that goes with it) and can count a face twice. Multi-tile grids get this from `_OuterTopology`, which resolves every corner to one canonical native `(face, j, i)`. Single-tile grids get it from `transports._insert_seam_twins()`: wherever consecutive corners are a physical edge but *not* index-adjacent — which is how the walk crosses a fold seam — it splices in the twin, so the crossing becomes one ordinary edge plus one zero-length twin edge, and zero-length edges emit no velocity face. This is the same construction as the `umaskutil`/`vmaskutil` masking of duplicated points in NEMO. The **sign flip** of an ORCA-style duplicated seam row (NEMO 4.2 manual, Appendix E) needs no handling of its own: the traversal direction and the stored velocity are read in the same index frame (geographically on both sides, for multi-tile), so a frame reversal cancels. A declared fold whose corner coordinates do not carry the fold symmetry has no twins to splice, and a crossing raises rather than silently mis-attributing.
+- **Seam representation contract (two axes, and `q`):** A section crossing a seam can be written *seam-explicit* (both indices of the shared corner, the step between them spanning no cell) or *twin-free* (one of them). **Input is liberal:** `uvindices_from_qindices()` accepts either and normalizes internally — physical coincidence is tested *before* index adjacency, because the two indices of a fold-seam corner sit far apart in the index lattice (column `i` against column `nx-i`). This is what lets a section arrive from `grid_section`, from a mask traced on the grid (`regionate.boundaries` deliberately keeps both coincident corners at a seam junction, and closes its loops with a coincident edge), or from saved indices reloaded off disk. **Output is strict:** the returned faces are all real velocity faces. So the two axes have different lengths, on purpose: the **corner** axis holds `i_c`/`j_c`/`f_c`/`lons_c`/`lats_c`, all 1:1 with each other (a repeated corner is *kept* — what is dropped is the face between the pair, never the corner); the **`sect`** axis holds the velocity faces and their own `lon`/`lat`. Relating them is what `uvindices["q"]` is for: `q[k]` is the index into the caller's corner arrays that face `k` starts at, so face `k` spans corners `q[k]` and `q[k]+1`. It is strictly increasing, skips exactly the steps that carried no face, and is carried as a coordinate on `convergent_transport()` and `extract_tracer()` output.
 - **Sign conventions:** For a *closed* section, `convergent_transport()` determines orientation (clockwise/counterclockwise) using stereographic projection and signed polygon area, then applies sign corrections so positive transport means "inward" (toward the enclosed polygon). For an *open* section there is no enclosing polygon, so `positive_in` is undefined; it instead uses the **left-of-transect** convention — `positive_in=True` makes positive transport point to the left of the section as traversed from the first to the last waypoint — and emits a `UserWarning`. `is_section_counterclockwise()` is only consulted in the closed case.
 - **xgcm.Grid dependency:** The package relies heavily on `xgcm.Grid` for grid metadata (axis boundaries, coordinate positions, dataset access via `grid._ds`).
 

docs/source/algorithm.md

@@ -226,7 +226,8 @@ from sectionate.transports import uvindices_from_qindices
 
 uv = uvindices_from_qindices(grid, i_c, j_c)
 # uv["var"] is "U" or "V" per face; uv["i"], uv["j"] index that velocity;
-# uv["Xinc"], uv["Yinc"] record the direction of travel through each face.
+# uv["Xinc"], uv["Yinc"] record the direction of travel through each face;
+# uv["q"] says which step of (i_c, j_c) each face came from.
 ```
 
 Whether a face is a `U` or a `V` point falls straight out of which way the step went, though
@@ -239,6 +240,29 @@ staggering (`outer`, `right` or `left`) and is read from the grid rather than as
 "Up to" `N-1`, because a consecutive pair that resolves to the same physical point — the twin
 corners of a seam again — spans no cell and carries no flux, so it emits no face at all.
 
+Twins matter here for a second reason. Naming the face between two corners assumes they are
+neighbours *in index space*, and across a bipolar fold that can fail: the walk may leave one
+index representation of a seam corner and land next to the *other*, a real physical edge whose
+endpoints are far apart in `i`. Attributing that step from the source corner alone would name
+the mirrored column's face, with the sign that goes with it. So before converting the chain,
+`uvindices_from_qindices` resolves corner identity — multi-tile grids canonicalise every corner
+through the outer-lattice topology, single-tile grids splice the skipped twin back into the
+chain — after which every remaining step is an ordinary adjacent one, and the zero-length twin
+edge it introduces drops out by the rule above rather than being counted twice.
+
+Which representation the chain arrives in does not matter. A section may already carry both
+indices of a seam corner — that is what the walk returns across a periodic wrap, and what a
+boundary traced on a cell mask returns at every seam junction — or only one, which is what it
+returns across a fold. Either way the faces come out the same, because coincidence is tested
+before adjacency: a pair that is *already* the zero-length twin edge is passed straight to the
+drop rather than being mistaken for a step needing repair.
+
+The consequence is that faces and corner steps are not one-to-one, so `uv["q"]` says which is
+which: face `k` spans corners `q[k]` and `q[k]+1`, and the steps `q` skips are exactly the ones
+that carried no flux. It rides along on `convergent_transport` and `extract_tracer` output, so
+a transport can be put back onto the section's corners — to plot it against them, to slice out
+the part belonging to one stretch of the path, or to find a child section inside its parent.
+
 `transports.convergent_transport` then accumulates the signed normal transport through those
 faces. For a **closed** section it works out the traversal orientation and signs everything so
 that positive means *into* the enclosed region; for an **open** section there is no inside, so