Topology-driven section-finding + outer-lattice multi-tile transports (from xgcm metadata)

[hdrake]

Consolidated PR. #48 (outer-lattice-topology) built directly on this branch; its work has been folded in here by fast-forward so review and merge happen in one place, and #48 is now closed as superseded. The outer-lattice topology + exact padded-transport work from #48 is summarized under Outer-lattice multi-tile topology below.

What & why

Section-finding used to hard-code a handful of grid topologies behind a topology= string ("latlon"/"cartesian"/"MOM-tripolar"), and the transport/tracer pipeline assumed a single logically-rectangular (Y, X) tile. This PR drops the string and instead reads topology from the xgcm.Grid the caller already supplies — its per-axis boundary and its face_connections — so periodicity, the bipolar north fold, and inter-tile connectivity are inferred rather than hand-coded. That single change extends both section-finding and transports/tracers to multi-tile grids (lat-lon-cap, cubed-sphere) and to all three native C-grid vorticity staggerings, with real-world example notebooks that double as validation.

Changes

Topology from xgcm metadata — no more topology=. Neighbors follow from each axis' padding: "periodic" → wrap, fill/extend wall → clip, and a single-tile bipolar north fold {"Y": {"fold": ...}} → crossed natively via xgcm's fold padding. The topology= argument and the hard-coded north-fold index formula are gone. gridutils.build_neighbor_maps builds the topology-aware neighbor maps the pathfinder consumes.

Multi-tile geometry via face_connections. build_neighbor_maps builds index-valued corner arrays and reads their cross-tile neighbors with any axis rotation/reversal already applied. A reciprocity check raises rather than returning silently-wrong neighbors. Section state now carries a face index, which grid_section returns for multi-tile grids.

All three C-grid staggerings, natively. Vorticity may sit at "outer" (MOM6 symmetric, M+1×N+1), "right" (MOM6 non-symmetric, M×N), or "left" (MITgcm/ECCO incl. lat-lon-cap, M×N). They are distinguished only by a per-position velocity index offset (gridutils.corner_position/corner_offset); no up-front conversion to a symmetric grid is needed.

Deterministic walk + curve option. The section walk is deterministic and direction-independent — it admits neighbors strictly closer to the endpoint plus any seam twin of the current cell, breaking ties by index — so periodic/fold-seam crossings no longer depend on floating-point rounding. grid_section gains curve=: "great circle" (default geodesic) or "latitude circle" (constant-latitude paths that march in longitude).

Transports / tracers across seams. Each section edge's normal velocity is read in a single face's frame (source-first dual-anchor), so no velocity is ever rotated across a seam. The per-edge sign is computed geometrically — "does the stored velocity's positive direction point to the left of the section's direction of travel?" — which stays consistent across rotated faces without any global orientation factor. convergent_transport/extract_tracer pick the contributing face pointwise.

Outer-lattice multi-tile topology (folded in from #48)

Rebuilds all multi-tile corner topology on the outer (shared-corner) lattice so rotated/reversed seams, cube-vertex junctions, boundary folds, and grid cuts become ordinary graph nodes — each with one canonical native (face, j, i) — and adds an exact, xgcm-independent transport-padding path.

• Outer-lattice corner graph (_OuterTopology). Replaces the earlier _correct_seam_neighbors (which walled off twisted 4-face cube-vertex junctions): reconstruct each face's full outer lattice, resolve every extended slot to the native corner storing its physical point (exact cell-identity matching, plus extrapolate+snap at face corners/cuts), merge slots into physical nodes, and project the global corner graph back onto native-frame neighbor maps. The walker and transport attribution consume these maps unchanged, preserving the public native (i_c, j_c, f_c) contract. On ECCO LLC90 the 8 twisted cube-vertices become ordinary 4-valent corners, the 3 Arctic-cap vertices 3-valent, and the south-pole fold + lon=−115 cut unify into single well-defined nodes.
• Exact padded transports. _OuterTopology.padded_transports extends native staggered transports to the outer lattice by resolving each missing edge slot to the stored velocity of that physical edge — a topology-exact, xgcm-independent alternative to xgcm.pad(..., other_component=...). On ECCO LLC90 the cell convergence sums to exactly zero globally. This remains the exact path even with the correct vector pad (xgcm#749), because no halo can supply edges stored on no face.
• Points stored on no face raise an error. A cubed-sphere stores 6N² corners for 6N²+2 points; LLC90 adds its 4th Arctic-cap vertex and the two sides of the lon=−115 cut. These become edge-less nodes — sections/traced boundaries through them raise an error rather than return invalid indices.
• Geometry-first 'left' cubed-sphere fixture. A gnomonic cubed-sphere in native 'left' staggering with face_connections from corner coincidence; transports from a corner streamfunction are exactly non-divergent, giving cell-by-cell convergence that sums to zero, section transport = endpoint streamfunction difference, and exactly zero net transport around a stored cube vertex.

Examples

• examples/4_sections_on_global_tripolar_grid.ipynb — sections traced across the bipolar north fold on a global single-tile tripolar grid, including a curve="latitude circle" global section.
• examples/5_MOC_transports_ECCOv4r4.ipynb — end-to-end depth-space global meridional overturning streamfunction from ECCOv4r4 on its native 13-tile lat-lon-cap (LLC90) grid: rotated multi-tile seam crossings + native "left" staggering + latitude-circle sections, on real model output. Its small ECCO data subset is fetched from Zenodo (10.5281/zenodo.21051424) with MD5 checks, so it runs with no NASA Earthdata login.

Verification

• Rotated multi-tile seams. Streamfunction-difference test on a genuinely rotated 2-face grid: transport across the rotated seam must equal psi(P2) − psi(P1), a known, non-zero, rotation-independent answer needing no model data. The ECCOv4r4 example is the real-grid counterpart.
• Bipolar north fold. Dedicated single-tile tests in test_section_fold.py: fold neighbor connectivity and a real-grid Arctic-crossing latitude circle that only completes once the grid declares its {"Y": {"fold": ...}} boundary.
• Split-grid equivalence. A 2-face grid is a flat grid cut in half; transport matches the uncut grid to rtol 1e-12 for both symmetric and non-symmetric axis-preserving seams.
• Existing single-tile tests pass unchanged; the geographic sign reproduces the previous Xinc/Yinc results to machine precision.
• Reversibility, "left"-staggering, and corner-case coverage across the new suites.
• Full pytest: 61 passed, 0 skipped against released xgcm 0.10.1, deterministic across PYTHONHASHSEED. All five example notebooks re-executed cleanly.

Dependency / CI note

The fold + multi-tile paths need xgcm >= 0.10.1 (now released on PyPI), which introduces the bipolar fold xgcm#711, the face-connection padding fix #712, and the bare-DataArray other_component vector-pad fix xgcm#749. padded_transports is independent of #749; the vector Grid.diff/Grid.interp used elsewhere benefits from it. Because 0.10.1 is a plain PyPI release, pip install -e . pulls it in — the previous --force-reinstall --no-deps install of the unreleased hdrake/xgcm@dev-v1.0.0 dev build is gone from CI. Fold/multi-tile tests skip on xgcm < 0.10.1; everything else runs on the same floor.

Note: the C-grid boundary-condition API rename lands in the 0.10.x line, not a later v1.0.0: Grid(boundary=) → padding=, Axis.boundary → .padding, pad(boundary=) → padding=, and periodic= is removed (periodic=False → padding="fill"). Sectionate source, tests, example loaders, and notebooks target this API; values and semantics are unchanged. (The remaining xgcm v1.0.0-line work not yet in 0.10.1 is not relevant to this PR.)

Merged with master

master has been merged into this branch (a merge, not a rebase, so the review history is preserved). Three conflicts, resolved as follows:

• pyproject.toml — kept this branch's xgcm >= 0.10.1 over master's xgcm == 0.9.*. This is the one dependency change the PR still carries.
• sectionate/version.py — took master's dynamic-versioning shim (Derive the version from the git tag instead of a checked-in file #50). The version now comes from the git tag via hatch-vcs, so the __version__ literal this branch carried is gone and no version string remains in the tree.
• .gitignore — union of both sides: this branch's data/*.nc, *.tar.gz, .DS_Store, plus master's generated sectionate/_version.py.

Because the versioning work is now merged in rather than diverged from, it no longer appears in this PR's diff — pyproject.toml here is a one-line dependency change. conda/meta.yaml and .github/workflows/ci.yml merged cleanly, keeping master's hatch-vcs host requirement and fetch-depth: 0 alongside this branch's xgcm floor.

CLAUDE.md (new in this PR) now also specifies one conda environment per branch/worktree, named docs_env_sectionate_<branch-or-worktree-name>: branches differ in their xgcm floor, so a single shared environment silently tests the wrong versions. Its Definition of Done applies before pushing as well as committing.

Known limitations (guarded, not silent)

• A tripolar/bipolar fold expressed as a face self-connection in face_connections raises NotImplementedError — build the grid as a single tile with the fold as a boundary ({"X": "periodic", "Y": {"fold": "corner"}}), which is supported natively.

---

AI-assisted PR, disclosed per the project AI Usage Policy: commits carry Co-Authored-By: trailers; the author has reviewed and validated the full diff. 🤖 Drafted with Claude Code.

[geoffstanley]

This looks great! I poked around and have a few small comments or questions:

• "Bipolar-fold and reversibility tests; the ECCO MOC example reproduces the expected global overturning structure on a real lat-lon-cap grid." this seems like a typo here, as the lat-lon-cap grid of ECCO isn't testing the bipolar north fold. (Are more bipolar north fold unit tests needed?)
• Thanks for implementing curve="latitude circle"!
• _left_sign tests the sign of a cross product of two vectors in lat-lon space. I worry this won't work near the geographic poles, or at the dateline (though your comments suggest it would; you may have thought about this more than me). The approach I quickly implemented was to pick the sign of U or V based on the F point grid indices alone: If the current F point indices (ic,jc) and previous F point indices (ip, jp) have ic==ip then we want U not V, and if jc = jp + 1 then take +U from (ic,jc), otherwise take -U from (ip,jp); similarly for taking V. This requires inserting duplicate F points in the section's list of indices where it crosses a boundary like the north fold or dateline, so the same geographic point but with different F grid indices is represented twice. Food for thought.

[hdrake]

These comments need to be addressed before this PR can be merged. Some of the comments are minor suggestions. Some are questions, which should be answered directly as a response to the comment, which @hdrake will either resolve or follow up on manually. Any questions that require large changes (e.g. changing at which step in the sectionate workflow the artificial cross-tile faces with zero physical length get dropped) should be implemented as separate PRs back into this PR.

[hdrake]

This is now outdated, should use v3.0.0 10.5281/zenodo.21479854; comments like this should not hard-core things like this, since the place where it should be tracked is the ZENODO_RECORD_ID below.

[hdrake]

Updated to v3.0.0 and rewrote the paragraph so it no longer carries a version DOI at all — it now cites only the version-independent concept DOI and says explicitly that the version being downloaded is pinned by ZENODO_RECORD_ID below, which is the single place that version is recorded.

Worth flagging that the bump was more than a number: v3.0.0 concatenates each collection along time into one file per collection, so the twelve monthly volume-flux entries in ECCO_FILE_MD5 collapse to a single 2010 file with a new checksum. I downloaded it and checked that its MD5 matches the published one, that it carries the same UVELMASS/VVELMASS on the same dims with the same twelve monthly timestamps, and that load_ECCO_MOC_grid builds utr/vtr with identical dims and shapes. The geometry file's checksum is unchanged.

In hdrake#9.

Written by Claude Code.

[hdrake]

Up to v3.0.0

[hdrake]

Done, and I removed the version marker rather than bumping it, for the same reason as the docstring above: ZENODO_RECORD_ID should be the only thing that records which version we pull. The comment now explains what the two constants are for and notes that changing the record id means re-checking the checksum table, since each version publishes its own file list — which is exactly what bit here.

In hdrake#9.

Written by Claude Code.

[hdrake]

Avoid weird jargon like " -- reproducing clip-to-edge". Either explain it in sectionate's own terms or in plain language.

[hdrake]

Dropped, in both places it appeared. It now says the thing plainly instead: at a closed boundary a corner is recorded as its own neighbour, so a step in that direction stays put, never gets the walk closer to the target, and is therefore never taken — the section stops at the boundary rather than running off the array.

In hdrake#9.

Written by Claude Code.

[hdrake]

Add a comment clarifying that this is the single tile case.

[hdrake]

Added, right after the multi-tile early returns. It notes that a single-tile grid's whole topology is its axis boundary conditions — periodic wrap, closed edge, north fold — that xgcm's halo padding already encodes all three, and that these maps come back with fmap = None.

In hdrake#9.

Written by Claude Code.

[hdrake]

The grid cut terminology is a bit confusing. It should be defined more clearly. The other cases are more obvious.

[hdrake]

Defined it in the class docstring and used it consistently from there. A grid cut is a seam along which the mesh has been slit open, so one line of physical corner points appears twice in the lattice, as the two coincident "lips" of the cut.

The part that makes it different from a face seam is now stated up front, because that is the bit that was doing the unexplained work: a cut is not declared in face_connections, so neither lip enters the other's halo and nothing combinatorial reveals that they are the same points — only their coordinates do. That is why the cut is the one place a coordinate fallback is needed at all. LLC90's 65°E/115°W slit under Antarctica is given as the example.

In hdrake#9.

Written by Claude Code.

[hdrake]

As mentioned in https://github.com/MOM6-community/sectionate/pull/47/changes#r3753380399, I think these should be dropped way earlier so that the public API never surfaces these artificial faces/index pairs.

[hdrake]

Moved into section-finding. The two halves were doing different jobs, and both belong there: resolving each corner to its canonical native (face, j, i) is what lets the face attribution look a corner up in either face's frame, and dropping consecutive duplicates is the fix you are asking for.

Both now happen in section.drop_repeated_corners, which needs the xgcm.Grid — so it is called from grid_section, from GriddedSection.__init__ when indices are passed in directly (covering .copy() and hand-built sections), and from load_gridded_section. That last one is what makes removing the old filter safe: sections saved before this change still contain both identities, and are normalised on load rather than silently mis-derived.

In hdrake#11.

Written by Claude Code.

[hdrake]

As mentioned in https://github.com/MOM6-community/sectionate/pull/47/changes#r3753380399, I think these faces should be dropped immediately after they are used in the section-finding algorithm.

[hdrake]

Done — it now happens immediately after the walk instead of at the end of the face derivation. The reason it lived in uvindices_from_qindices was to cover reloaded sections too; that is now handled by normalising on load, so this filter is unreachable and removed.

Which of the two identities to keep turned out to matter. Keeping the later one — the frame the path continues in — makes both flanking faces come out with exactly the indices they had before, so no transport number moves. I verified that rather than assuming it: reference values captured from the unmodified base checkout, compared against the branch across six ECCO LLC90 latitude circles and three MOM6 tripolar-fold sections. Exactly one recorded value changed, n_corners 308 -> 307 on the zonal periodic section, which is the point of the change. Face counts, index sums, index heads, U/V counts and the overturning streamfunction are all identical.

In hdrake#11.

Written by Claude Code.

[hdrake]

Remove "(e.g. MOM6)"

[hdrake]

Change "Only structured grids are supported" to supports any structured model whose grid can be defined by a xgcm.Grid object.

[hdrake]

Removed.

In hdrake#9.

Written by Claude Code.

[hdrake]

Done — it now reads that sectionate supports any structured model whose grid can be described by an xgcm.Grid object.

In hdrake#9.

Written by Claude Code.

[hdrake]

Remove: "Branches diverge in their dependencies — the xgcm floor in particular — so a single
shared environment silently tests the wrong versions."

[hdrake]

Removed.

In hdrake#9.

Written by Claude Code.

[hdrake]

Remove all of this, which seems unnecessary.

[hdrake]

Removed, all of it. One dangling reference was left behind in the following paragraph, which cited xgcm < 0.10.1 as a skip cause; that now refers to the floor in pyproject.toml instead of restating a version. Nothing else in CLAUDE.md, README.md or docs/ depended on the removed text.

In hdrake#9.

Written by Claude Code.

[hdrake]

Disclosure, per the project's AI usage policy: the replies I am posting on this review's comment threads were drafted with Claude Code, as were the three follow-up PRs linked from them. Every claim in them was checked against a run rather than asserted, and where a measurement contradicted what the PR already said, the PR is being corrected rather than the claim defended (see the padded_transports thread).

The review's comments are being addressed as three PRs back into topology-driven-neighbors, as you asked:

• Review follow-up: wording, docstrings, ECCO record v3.0.0, generic topology check hdrake/sectionate#9 — wording, docstrings, CLAUDE.md, the ECCO Zenodo record (now v3.0.0), and the generic topology check.
• Always take the shortest path between section waypoints hdrake/sectionate#10 — always take the shortest path between waypoints; one ambiguity error for both curve types.
• A third PR for dropping the artificial zero-length cross-tile faces immediately after section-finding, which is the large change you singled out.

Each runs the full suite with 0 skips on its own branch. The example notebooks are re-executed once, on the integrated branch, after all three land.

Written by Claude Code.

[hdrake]

• The approach I quickly implemented was to pick the sign of U or V based on the F point grid indices alone: If the current F point indices (ic,jc) and previous F point indices (ip, jp) have ic==ip then we want U not V, and if jc = jp + 1 then take +U from (ic,jc), otherwise take -U from (ip,jp); similarly for taking V. This requires inserting duplicate F points in the section's list of indices where it crosses a boundary like the north fold or dateline, so the same geographic point but with different F grid indices is represented twice. Food for thought.

@geoffstanley, I looked into your approach but it got a bit messy when it came to multi-tile grids like ECCO's LLC because the meaning of U and V flip across the boundaries of some tiles. I like that with the geometric approach you don't need to worry about such things because it all comes for free from the coordinates.

[geoffstanley]

@geoffstanley, I looked into your approach but it got a bit messy when it came to multi-tile grids like ECCO's LLC because the meaning of U and V flip across the boundaries of some tiles. I like that with the geometric approach you don't need to worry about such things because it all comes for free from the coordinates.

I've continued a bit further with my approach and it actually seems to be working very well. I'll have some tests done soon to confirm either way though.

When crossing the north fold (or perhaps other tile boundaries), how do you handle the sign flip of the velocity components? In the ORCA north fold for example on a T-point pivot, the right half of the uppermost U grid is duplicate and carries a sign flip, while the entire uppermost V grid is duplicate and carries a sign flip. For an F-point pivot, the U grid has no duplicates but the V grid uppermost right half is duplicate with a sign flip. See Appendix E of the NEMO 4.2 manual.

If your list of F-point indices for the section does not include duplicates (i.e. different integer indices that represent the same geographic point), then how do you avoid double counting, and how do you capture the sign flip for certain grids if your transect happens to be on that side (the right half)? It will require considerable knowledge of the grid, which will get hard to generalize to many grids (LLC, cube-sphere, bifold, etc).

What I have done is to purposely include duplicate (i,j) points on the F grid, and then when calculating transport. Specifically, I run the GriddedSection algorithm twice: once forwards; and again with the section reversed (using Section.reverse()) followed by reversing the list of (i,j) F grid indices. If the section does not cross any duplicate points (on the F grid), then the two lists are identical and nothing needs to be done. But when the two lists disagree, I merge the two lists together, taking all the points that agree as they are, but where they differ, I take all the points from the first list first, then all the points from the second list second. The result is a list of indices that are all on the same "face" ending with points that are in the duplicate zone, followed by a jump to a list of indices that begin with points in the duplicate zone and are all on the other face.
(In practice, the two cases are handled identically; no check needs to be done on whether the two lists are identical.)

Then, accumulating transport across the section happens naturally with the above algorithm, with the simple modification that I mask out the contributions from duplicated U or V points. In NEMO, those masks are given by 2D fields called umaskutil and vmaskutil (not the standard 3D masks umask and vmask). Note, in NEMO, the non-duplicated points always carry the normal (positive) sign.

My only assumption is that the "jump" involves two points that are not adjacent to each other in (i,j) space (ie same i and j differs by 1, or same j and i differs by 1). If it did, then transport could be double counted there. With the bifold, this should not happen as the pivot points are in land -- so even if you made a section that went through the pivot points, you'd be double counting zero transport, so no problem. For the LLC or cube sphere or anything else involving multiple tiles, the logical coordinates are (tile, i, j) and so the tile coordinate would be different, so again no double counting would happen.