From 7b93104f9d9c1374c4e80a282b8ddee734f98907 Mon Sep 17 00:00:00 2001 From: Henri Drake Date: Mon, 10 Aug 2026 15:27:28 -0700 Subject: [PATCH 1/4] Always take the shortest path between section waypoints Sectionate now resolves every section segment to the shortest path between its two waypoints, for both `curve="great circle"` and `curve="latitude circle"`. Raw longitudes are no longer read as a request to travel the long way round: a latitude-circle segment written 0 -> 270 runs 90 degrees west, matching what a great circle already did. The old "each segment must span less than 180 degrees" ValueError is gone. `_check_segment_span` now raises a single, shared ambiguity error for both curve types, and only for genuinely ill-posed segments: - endpoints exactly half a circle apart, where the two ways round are equally short so the shortest path is not unique (antipodal for a great circle, +/-180 degrees of longitude for a latitude circle); - for "latitude circle" only, oblique endpoints that differ in both latitude and longitude, since no circle of constant latitude passes through them. Latitude-circle segments must be zonal (shared latitude) or meridional (shared longitude); a meridional segment is the connector used to join zonal arcs at different latitudes and walks straight down the meridian. Previously an oblique latitude-circle segment was traced as an L-shaped path -- zonal along the starting latitude, then meridional -- which was direction-dependent and, when it crossed the periodic seam, degenerated into a staircase of zero-length seam steps. Co-Authored-By: Claude Opus 5 (1M context) --- CLAUDE.md | 2 +- sectionate/section.py | 132 ++++++++++++++----- sectionate/tests/test_section_cornercases.py | 94 +++++++++++-- 3 files changed, 185 insertions(+), 43 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 3d90752..868126a 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -122,7 +122,7 @@ 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"`). 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` raises only for ill-posed segments: endpoints exactly half a circle apart (no unique shortest path, one error shared by both curves), or, for `"latitude circle"`, oblique endpoints differing in both latitude and longitude (no constant-latitude circle passes through them; such segments must be zonal or meridional). 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). diff --git a/sectionate/section.py b/sectionate/section.py index b730346..9f9040a 100644 --- a/sectionate/section.py +++ b/sectionate/section.py @@ -17,6 +17,11 @@ # index, so the path is independent of platform and of travel direction. WALK_DEVIATION_ATOL = 1.e-9 +# Angular tolerance, in degrees, for the well-posedness checks a section segment must +# pass (`_check_segment_span`). It exists only to absorb floating-point round-off in the +# waypoint coordinates: endpoints further apart than this really are distinct. +SEGMENT_ATOL_DEG = 1.e-9 + class Section(): """A named hydrographic section""" def __init__(self, name, coords, children = {}, parent = None): @@ -263,8 +268,18 @@ def grid_section(grid, lons, lats, curve="great circle"): Latitudes, in degrees (in range [-90, 90]), of consecutive vertices defining a piece-wise section. curve: str Curve followed between consecutive vertices: "great circle" (default, the geodesic) or - "latitude circle" (constant latitude, marching in longitude). Each segment must span - less than 180 degrees, otherwise the direction is ambiguous and a ValueError is raised. + "latitude circle" (constant latitude, marching in longitude). + + Each segment follows the **shortest** path between its two vertices; raw longitudes + are never taken as a request to go the long way round, so a "latitude circle" + segment from 0 to 270 degrees runs 90 degrees *west*. Encircle the globe by giving + intermediate vertices (e.g. 0 -> 120 -> 240 -> 360), which is also what says which + way round it goes. A ValueError is raised only if a segment is ill posed: its + endpoints are exactly half a circle apart (so the shortest path is not unique), or + -- for "latitude circle" -- they are oblique, differing in both latitude and + longitude, so that no circle of constant latitude passes through them. Latitude- + circle segments must therefore be zonal (shared latitude) or meridional (shared + longitude); see `_check_segment_span`. Returns ------- @@ -351,6 +366,13 @@ def create_section_composite( Topology-aware neighbor maps from `sectionate.gridutils.build_neighbor_maps` (single- or multi-tile). Sections are always built from an `xgcm.Grid`, so these are always supplied; the usual entry point is `sectionate.grid_section`. + curve: str + Curve followed between consecutive vertices: "great circle" (default, the geodesic) + or "latitude circle" (constant latitude, marching in longitude). Every segment + follows the shortest path between its two vertices, so a piece-wise section that + encircles the globe needs intermediate vertices (e.g. 0 -> 120 -> 240 -> 360) -- + both to go the long way round and to say which way round. Each segment is checked + by `_check_segment_span`, whose docstring gives the full contract. RETURNS: ------- @@ -452,43 +474,68 @@ def create_section(gridlon, gridlat, lonstart, latstart, lonend, latend, neighbo curve=curve, ) +def _wrapped_dlon(lon1, lon2): + """Signed longitude change from `lon1` to `lon2`, wrapped into [-180, 180). + + This is the change along the *shortest* way round: a raw change of +270 degrees comes + back as -90. Endpoints that coincide modulo 360 degrees (e.g. a 360 -> 0 loop + closure) give 0. + """ + return (lon2 - lon1 + 180.) % 360. - 180. + + def _check_segment_span(lon1, lat1, lon2, lat2, curve): - """Raise if a section segment's direction between its endpoints is ambiguous. - - A unique *directed* shortest path needs the two endpoints to be less than half a - circle apart: - - - "latitude circle": the longitude change must be less than 180 degrees; at or beyond - that the east/west direction is equally far either way (and a full circle is - degenerate). Longitudes are taken as given, so write a >180-degree arc with one or - more intermediate waypoints (e.g. split 0 -> 270 into 0 -> 135 -> 270). Endpoints - that coincide modulo 360 degrees (e.g. a 360 -> 0 loop closure) describe a - zero-length segment, not a full circle, and are allowed. - - "great circle": the endpoints must not be (near-)antipodal, where infinitely many - geodesics connect them. + """Raise if the segment between two consecutive waypoints is not well posed. + + Sectionate always follows the **shortest** path between consecutive waypoints, for + both curve types. Raw coordinates are never read as an instruction to take the long + way round: a ``curve="latitude circle"`` segment written 0 -> 270 degrees travels 90 + degrees *west*, not 270 degrees east, exactly as a great-circle segment would. To go + the long way round, say so with intermediate waypoints (e.g. 0 -> 120 -> 240 -> 360). + + Two things can still leave a segment ill posed: + + 1. **The shortest path is not unique.** This happens only when the endpoints are + exactly half a circle apart -- antipodal under "great circle" (infinitely many + geodesics join them), or exactly +/-180 degrees of longitude apart under + "latitude circle" (east and west are equally far). A single ambiguity error + covers both curves. Split such a segment with an intermediate waypoint. + 2. **The requested curve does not pass through both endpoints.** For "latitude + circle" only: a circle of constant latitude through both endpoints exists only if + they share a latitude. Segments that instead share a longitude are also accepted + -- these are the meridional connectors that join zonal arcs at different + latitudes, and a meridian is unambiguous. An *oblique* segment, differing in both + latitude and longitude, lies on no latitude circle; use ``curve="great circle"`` + for it, or split it into a zonal leg and a meridional leg. """ if curve == "latitude circle": - dlon = abs(lon2 - lon1) - # Wrapped longitude change in (-180, 180]. A magnitude near 0 means the endpoints - # coincide modulo 360 degrees -- a zero-length segment (e.g. a 360 -> 0 loop - # closure) with no east/west ambiguity -- so it is allowed even though the raw - # dlon is a multiple of 360. Genuine arcs keep the "taken as given" rule below. - wrapped = (lon2 - lon1 + 180.) % 360. - 180. - if abs(wrapped) > 1.e-9 and dlon >= 180.: + # Signed shortest-way-round longitude change. It is also the separation that + # decides ambiguity, because a latitude-circle segment travels only in longitude. + dlon = _wrapped_dlon(lon1, lon2) + zonal = abs(lat2 - lat1) <= SEGMENT_ATOL_DEG + meridional = abs(dlon) <= SEGMENT_ATOL_DEG + if not (zonal or meridional): raise ValueError( - f"Latitude-circle segment from lon={lon1} to lon={lon2} spans {dlon} " - "degrees of longitude; each segment must span less than 180 degrees, " - "otherwise the east/west direction is ambiguous. Add intermediate " - "waypoints to subdivide longer arcs." + f"Latitude-circle segment from (lon={lon1}, lat={lat1}) to " + f"(lon={lon2}, lat={lat2}) is oblique: its endpoints differ in both " + "latitude and longitude, so no circle of constant latitude passes " + "through them. Use curve='great circle' for oblique segments, or split " + "this one into a zonal leg (same latitude) and a meridional leg (same " + "longitude)." ) + # A meridional segment runs along a single meridian, so `dlon` is 0 and the + # ambiguity test below is trivially passed: there is no east/west choice to make. + sep = abs(dlon) else: sep = np.rad2deg(distance_on_unit_sphere(lon1, lat1, lon2, lat2, R=1.)) - if sep >= 180. - 1.e-9: - raise ValueError( - f"Great-circle segment from ({lon1}, {lat1}) to ({lon2}, {lat2}) is " - f"(near-)antipodal (separation ~{sep:.4f} degrees); the geodesic " - "direction is ambiguous. Add an intermediate waypoint to disambiguate." - ) + + if sep >= 180. - SEGMENT_ATOL_DEG: + raise ValueError( + f"Segment from (lon={lon1}, lat={lat1}) to (lon={lon2}, lat={lat2}) has " + f"endpoints half a circle apart (separation ~{sep:.4f} degrees), so the two " + "ways round are equally short and the shortest path between them is not " + "unique. Add an intermediate waypoint to say which way the section goes." + ) def infer_grid_path_from_geo(lonstart, latstart, lonend, latend, gridlon, gridlat, neighbor_maps, curve="great circle"): @@ -638,11 +685,28 @@ def deviation(lon, lat): return (spherical_angle(lon2, lat2, lon1, lat1, lon, lat) + spherical_angle(lon1, lat1, lon2, lat2, lon, lat)) elif curve == "latitude circle": + # These two metrics measure progress purely in longitude and deviation purely in + # latitude, which is exactly right for the two segment shapes `_check_segment_span` + # admits for this curve, and for nothing else: + # - a *zonal* segment (endpoints on one parallel) marches in longitude at constant + # latitude, which is what `progress` drives and `deviation` holds it to; + # - a *meridional* segment (endpoints on one meridian) has lon2 == lon1, so + # `progress` is flat and admits no neighbor. The walk then falls through to the + # "no admissible forward move" branch below, which steps to the neighbor + # geodesically closest to the endpoint -- i.e. straight down the meridian. + # An oblique segment would instead walk zonally along the *start* latitude and only + # then close the latitude gap, an L-shaped path that depends on which end you start + # from; `_check_segment_span` rejects those up front rather than tracing one. def progress(lon, lat): - # monotonic in |delta-lon| over each (sub-180-degree) segment; direction-symmetric. + # sin^2(delta-lon/2) is the haversine of the longitude gap: periodic in 360 + # degrees and monotonic in |delta-lon| up to 180, so it measures the *shortest* + # way round regardless of how the endpoint longitudes were written (0 -> 270 + # descends westward just as 0 -> -90 does). Direction-symmetric. return np.sin(np.deg2rad((lon - lon2) / 2.)) ** 2 def deviation(lon, lat): - # angular distance off the constant-latitude curve through the endpoints (radians). + # angular distance off the constant-latitude curve through the endpoints + # (radians). Flat between the two endpoint latitudes, which absorbs the + # sub-cell mismatch left when each endpoint snaps to its nearest grid corner. return np.deg2rad(abs(lat - lat1)) + np.deg2rad(abs(lat - lat2)) else: raise ValueError( diff --git a/sectionate/tests/test_section_cornercases.py b/sectionate/tests/test_section_cornercases.py index 8299ba9..2ab51e2 100644 --- a/sectionate/tests/test_section_cornercases.py +++ b/sectionate/tests/test_section_cornercases.py @@ -1,4 +1,5 @@ import numpy as np +import pytest import xarray as xr import xgcm @@ -87,11 +88,88 @@ def test_latitude_circle_zero_length_closure(): ]) -def test_latitude_circle_long_arc_still_rejected(): - """The zero-length exemption must not relax the genuine-ambiguity guard: a real - arc spanning >= 180 degrees of longitude (taken as given) is still rejected.""" - from sectionate.section import _check_segment_span - import pytest - for lon2 in (180., 270., 200.): # genuine arcs, endpoints do NOT coincide - with pytest.raises(ValueError, match="spans"): - _check_segment_span(0., 0., lon2, 0., "latitude circle") +def test_latitude_circle_takes_shortest_path_west(): + """Sectionate always walks the *shortest* path between two waypoints, so a + latitude-circle segment written 0 -> 270 travels 90 degrees WEST rather than 270 + degrees east. Raw longitudes carry no intent to go the long way round.""" + from sectionate.section import grid_section, _check_segment_span + + # Unit level: a raw 270-degree change is a 90-degree westward one; not an error. + _check_segment_span(0., 0., 270., 0., "latitude circle") + + # End to end: 0 -> 270 leaves lon=0 westward across the periodic seam (the seam + # vertex carries both index 6 (lon 360) and index 0 (lon 0)); the endpoint snaps to + # the nearest corner at lon=240. The path never visits the eastern half. + i, j, lons, lats = grid_section(grid, [0., 270.], [0., 0.], curve="latitude circle") + assert np.all([ + modequal(i, np.array([0, 6, 5, 4])), + modequal(j, np.array([2, 2, 2, 2])), + modequal(lons, np.array([0., 360., 300., 240.])), + modequal(lats, np.array([0., 0., 0., 0.])), + ]) + assert not np.any(np.isin(np.mod(lons, 360.), [60., 120., 180.])) + + # The default great circle already behaved this way; both curves now agree. + i_gc, j_gc, lons_gc, lats_gc = grid_section(grid, [0., 270.], [0., 0.]) + assert np.all([modequal(i_gc, i), modequal(j_gc, j)]) + + +def test_half_circle_separation_is_the_only_ambiguity_error(): + """The single remaining ambiguity error covers both curve types: endpoints exactly + half a circle apart have two equally short paths between them, so the shortest path + is not unique. Anything short of that resolves silently to the shortest path.""" + from sectionate.section import grid_section, _check_segment_span + + ambiguous = [ + # (lon1, lat1, lon2, lat2, curve) + (0., 0., 180., 0., "latitude circle"), # exactly +180 degrees of longitude + (0., 0., -180., 0., "latitude circle"), # exactly -180 degrees of longitude + (10., 30., 190., 30., "latitude circle"), # 180 degrees away from lon=10 + (0., 0., 180., 0., "great circle"), # antipodal on the equator + (0., 30., 180., -30., "great circle"), # antipodal off the equator + ] + for lon1, lat1, lon2, lat2, curve in ambiguous: + with pytest.raises(ValueError, match="half a circle apart"): + _check_segment_span(lon1, lat1, lon2, lat2, curve) + + # Just short of half a circle is unambiguous for both curves. + _check_segment_span(0., 0., 179.9, 0., "latitude circle") + _check_segment_span(0., 0., 179.9, 0., "great circle") + + # The same single error reaches the user through the public API. + with pytest.raises(ValueError, match="half a circle apart"): + grid_section(grid, [0., 180.], [0., 0.], curve="latitude circle") + with pytest.raises(ValueError, match="half a circle apart"): + grid_section(grid, [0., 180.], [0., 0.]) + + +def test_latitude_circle_meridional_segment_walks_the_meridian(): + """A latitude-circle segment whose endpoints share a longitude is the meridional + connector used to join zonal arcs at different latitudes. It is unambiguous, so it + is allowed, and it walks straight down the meridian.""" + from sectionate.section import grid_section + + i, j, lons, lats = grid_section(grid, [60., 60.], [-80., 80.], curve="latitude circle") + assert np.all([ + modequal(i, np.array([1, 1, 1, 1, 1])), + modequal(j, np.array([0, 1, 2, 3, 4])), + modequal(lons, np.array([60., 60., 60., 60., 60.])), + modequal(lats, np.array([-80., -40., 0., 40., 80.])), + ]) + + +def test_latitude_circle_oblique_segment_rejected(): + """No circle of constant latitude passes through two points that differ in both + latitude and longitude, so such a segment is rejected rather than traced as an + L-shaped, direction-dependent path. Great circle handles it instead.""" + from sectionate.section import grid_section, _check_segment_span + + with pytest.raises(ValueError, match="oblique"): + _check_segment_span(0., 0., 120., 40., "latitude circle") + with pytest.raises(ValueError, match="oblique"): + grid_section(grid, [0., 120.], [0., 40.], curve="latitude circle") + + # The same waypoints are fine as a great circle, and as a zonal leg plus a + # meridional leg under latitude circle. + grid_section(grid, [0., 120.], [0., 40.]) + grid_section(grid, [0., 120., 120.], [0., 0., 40.], curve="latitude circle") From d0ef93b4b7f15ed8927e2b908dc24140d35390d9 Mon Sep 17 00:00:00 2001 From: Henri Drake Date: Tue, 11 Aug 2026 10:28:20 -0700 Subject: [PATCH 2/4] Decide the curve per segment, and make "latitude circle" mean it `curve` now takes three values instead of two: - "great circle" (default) -- unchanged, every segment follows the geodesic. - "latitude circle" -- every segment follows a parallel. A segment whose endpoints do not share a latitude lies on no circle of constant latitude, so it now raises. Previously meridional segments were quietly accepted here on the grounds that a meridian is unambiguous; they are not latitude circles, and accepting them hid a real failure (below). - "latitude and great circle" (new) -- decided segment by segment: segments whose endpoints share a latitude follow the parallel, all others follow the geodesic. This is what a section that is zonal in places and joined up by meridional or slanted legs elsewhere actually wants. This fixes a bug. The constant-latitude `progress` metric measures progress purely in longitude, so along a meridian it is flat: the walk admits no neighbour, falls through to the fallback branch, and on a real grid can fail to converge. `grid_section(grid, [0., 0.], [80., 60.], curve="latitude circle")` on ECCO LLC90 raised `RuntimeError: Should have reached the endpoint by now.` in that direction while 60->80 traced fine -- a direction-dependent failure on a section that never should have been routed to those metrics. Under "latitude circle" the segment is now refused with a clear message; under "latitude and great circle" it routes to the geodesic and traces identically in both directions. The per-segment choice is made once, in `infer_grid_path_from_geo` (which is called once per segment), from the requested waypoints rather than from the grid corners they snap to, and threaded down to the metric selection. `_is_constant_latitude` is the single classification behind both the legality check and the metric choice, so the two cannot disagree. Also, from review of the shortest-path change: - The segment tolerance was 1e-9 degrees, 0.11 mm. It is now a *classification* tolerance, and at that size a latitude that has been through float32 -- which is how ECCO and many models store corner coordinates -- reads as a different latitude and turns a zonal segment oblique. Split into two named constants: `CONSTANT_LATITUDE_ATOL_DEG` (1e-6 deg, ~11 cm, five orders of magnitude below any grid cell) and `HALF_CIRCLE_ATOL_DEG` (1e-9 deg), which measure different things. - The shared half-a-circle error quoted a bare "separation ~180.0000 degrees" for two points 222 km apart at 89N, because along a parallel the number is degrees of longitude and not degrees of arc. Still one error, but it now names the measure it is reporting. - The `180 - atol` guard is described for what it does: it catches endpoints written exactly half a circle apart. A hair short of that (0 -> 179.9999999 east, 0 -> 180.0000001 west) resolves by round-off with no error, so the docstring says to subdivide rather than to rely on the check. - The contract was restated in full in five places. It is now stated once, in `grid_section`, and cross-referenced from the rest; nothing points a user at a private name, since there is no rendered API page for one to land on. - Dropped the comment describing the L-shaped oblique latitude-circle path, which no longer reaches that code. Tests: all three options over zonal, meridional, oblique and exactly-half-a- circle segments; a 2-degree global fixture where the two metrics visibly differ (0->120 at 40N holds 40N for 61 points, while the geodesic bows to 60N over 83) -- the existing 7x5 fixture has corner rows 40 degrees apart and cannot show this; float32-derived waypoints classifying as constant-latitude; and the ECCO LLC90 meridional segment tracing identically in both directions. `test_oblique_latitude_circle_segments_reversible` is renamed: its legs are zonal and meridional, so under these rules it was never oblique. Co-Authored-By: Claude Opus 5 (1M context) --- CLAUDE.md | 2 +- sectionate/section.py | 213 +++++++++++------- sectionate/tests/test_ecco_llc90.py | 43 ++++ sectionate/tests/test_section_cornercases.py | 171 ++++++++++++-- .../tests/test_section_reversibility.py | 14 +- 5 files changed, 328 insertions(+), 115 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 868126a..11cbcf7 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -122,7 +122,7 @@ 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"`). 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` raises only for ill-posed segments: endpoints exactly half a circle apart (no unique shortest path, one error shared by both curves), or, for `"latitude circle"`, oblique endpoints differing in both latitude and longitude (no constant-latitude circle passes through them; such segments must be zonal or meridional). 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. 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). diff --git a/sectionate/section.py b/sectionate/section.py index 9f9040a..238d5b4 100644 --- a/sectionate/section.py +++ b/sectionate/section.py @@ -17,10 +17,27 @@ # index, so the path is independent of platform and of travel direction. WALK_DEVIATION_ATOL = 1.e-9 -# Angular tolerance, in degrees, for the well-posedness checks a section segment must -# pass (`_check_segment_span`). It exists only to absorb floating-point round-off in the -# waypoint coordinates: endpoints further apart than this really are distinct. -SEGMENT_ATOL_DEG = 1.e-9 +# The three curves a section can be asked to follow. "latitude and great circle" is not a +# curve in its own right: it resolves, segment by segment, to one of the other two (see +# `_segment_curve`). `grid_section` documents what each one means. +CURVES = ("great circle", "latitude circle", "latitude and great circle") + +# Angular tolerance, in degrees, for deciding whether a segment's two endpoints share a +# latitude -- the classification that says whether a segment follows a parallel. About +# 11 cm on Earth: five orders of magnitude finer than any grid cell, yet coarse enough to +# absorb round-off in waypoints that have been through single precision (a float32 +# latitude is off by ~4.e-7 degrees, and model corner coordinates are commonly stored as +# float32). +CONSTANT_LATITUDE_ATOL_DEG = 1.e-6 + +# Angular tolerance, in degrees, on the separation at which the two ways round a segment +# are equally long. This is a different quantity from the classification tolerance above, +# and it is deliberately tight: it catches endpoints *written* exactly half a circle +# apart (lon 0 -> 180) and nothing else. A segment even slightly short of half a circle is +# not ambiguous to the code -- which way round it runs is settled by the sign of the +# round-off, so 0 -> 179.9999999 goes east and 0 -> 180.0000001 goes west -- so subdivide +# such segments with an intermediate waypoint rather than relying on this check. +HALF_CIRCLE_ATOL_DEG = 1.e-9 class Section(): """A named hydrographic section""" @@ -267,19 +284,31 @@ def grid_section(grid, lons, lats, curve="great circle"): lats: list or np.ndarray Latitudes, in degrees (in range [-90, 90]), of consecutive vertices defining a piece-wise section. curve: str - Curve followed between consecutive vertices: "great circle" (default, the geodesic) or - "latitude circle" (constant latitude, marching in longitude). - - Each segment follows the **shortest** path between its two vertices; raw longitudes - are never taken as a request to go the long way round, so a "latitude circle" - segment from 0 to 270 degrees runs 90 degrees *west*. Encircle the globe by giving - intermediate vertices (e.g. 0 -> 120 -> 240 -> 360), which is also what says which - way round it goes. A ValueError is raised only if a segment is ill posed: its - endpoints are exactly half a circle apart (so the shortest path is not unique), or - -- for "latitude circle" -- they are oblique, differing in both latitude and - longitude, so that no circle of constant latitude passes through them. Latitude- - circle segments must therefore be zonal (shared latitude) or meridional (shared - longitude); see `_check_segment_span`. + Curve followed between consecutive vertices. One of: + + - "great circle" (default): every segment follows the geodesic. + - "latitude circle": every segment follows a circle of constant latitude, marching + in longitude. A segment whose endpoints do not share a latitude (to within + 1.e-6 degrees) lies on no such circle, so it raises a ValueError -- meridional + segments included. + - "latitude and great circle": decided per segment. A segment whose endpoints + share a latitude follows the parallel; every other segment follows the geodesic. + This is the option for a section that is zonal in places and joined up by + meridional or slanted legs elsewhere. + + Under every option each segment takes the **shortest** path between its two + vertices. Raw longitudes are never read as a request to go the long way round, so + a segment written 0 -> 270 along a parallel runs 90 degrees *west*. Encircle the + globe by giving intermediate vertices (e.g. 0 -> 120 -> 240 -> 360), which is also + what says which way round it goes. + + A ValueError is also raised for a segment whose endpoints are written exactly half + a circle apart -- antipodal for a geodesic segment, or exactly +/-180 degrees of + longitude apart for one along a parallel -- because then neither way round is the + shorter. That check only catches endpoints typed as exactly half a circle apart: + a hair short of it, which way the segment runs is decided by the sign of the + round-off, so subdivide near-half-circle segments with an intermediate vertex + rather than expecting to be warned about them. Returns ------- @@ -367,12 +396,10 @@ def create_section_composite( (single- or multi-tile). Sections are always built from an `xgcm.Grid`, so these are always supplied; the usual entry point is `sectionate.grid_section`. curve: str - Curve followed between consecutive vertices: "great circle" (default, the geodesic) - or "latitude circle" (constant latitude, marching in longitude). Every segment - follows the shortest path between its two vertices, so a piece-wise section that - encircles the globe needs intermediate vertices (e.g. 0 -> 120 -> 240 -> 360) -- - both to go the long way round and to say which way round. Each segment is checked - by `_check_segment_span`, whose docstring gives the full contract. + Curve followed between consecutive vertices: "great circle" (default), + "latitude circle", or "latitude and great circle". Each segment is resolved and + checked independently; `sectionate.grid_section` documents what the three options + mean and which segments they reject. RETURNS: ------- @@ -484,59 +511,70 @@ def _wrapped_dlon(lon1, lon2): return (lon2 - lon1 + 180.) % 360. - 180. +def _is_constant_latitude(lat1, lat2): + """Whether a segment's two endpoint latitudes agree, and so lie on one parallel. + + The single classification used both to accept or reject a "latitude circle" segment + and to pick the metrics a segment is walked with, so those two can never disagree. + """ + return abs(lat2 - lat1) <= CONSTANT_LATITUDE_ATOL_DEG + + +def _segment_curve(lat1, lat2, curve): + """The curve one segment actually follows: always "great circle" or "latitude circle". + + Resolves the section-wide `curve` request (one of `CURVES`) for a single segment, and + raises ValueError if it is not a recognized request. Passing an already-resolved value + returns it unchanged, so this is safe to apply more than once along a call chain. + """ + if curve in ("great circle", "latitude circle"): + return curve + if curve == "latitude and great circle": + return "latitude circle" if _is_constant_latitude(lat1, lat2) else "great circle" + raise ValueError( + f"curve must be one of {', '.join(repr(c) for c in CURVES)}; got {curve!r}." + ) + + def _check_segment_span(lon1, lat1, lon2, lat2, curve): - """Raise if the segment between two consecutive waypoints is not well posed. - - Sectionate always follows the **shortest** path between consecutive waypoints, for - both curve types. Raw coordinates are never read as an instruction to take the long - way round: a ``curve="latitude circle"`` segment written 0 -> 270 degrees travels 90 - degrees *west*, not 270 degrees east, exactly as a great-circle segment would. To go - the long way round, say so with intermediate waypoints (e.g. 0 -> 120 -> 240 -> 360). - - Two things can still leave a segment ill posed: - - 1. **The shortest path is not unique.** This happens only when the endpoints are - exactly half a circle apart -- antipodal under "great circle" (infinitely many - geodesics join them), or exactly +/-180 degrees of longitude apart under - "latitude circle" (east and west are equally far). A single ambiguity error - covers both curves. Split such a segment with an intermediate waypoint. - 2. **The requested curve does not pass through both endpoints.** For "latitude - circle" only: a circle of constant latitude through both endpoints exists only if - they share a latitude. Segments that instead share a longitude are also accepted - -- these are the meridional connectors that join zonal arcs at different - latitudes, and a meridian is unambiguous. An *oblique* segment, differing in both - latitude and longitude, lies on no latitude circle; use ``curve="great circle"`` - for it, or split it into a zonal leg and a meridional leg. + """Validate one section segment and return the curve it follows. + + `curve` is the section-wide request; the return value is what it resolves to for this + segment, either "great circle" or "latitude circle". Raises ValueError if the segment + is ill posed under `curve` -- see `sectionate.grid_section` for the rules and for what + to write instead. """ - if curve == "latitude circle": - # Signed shortest-way-round longitude change. It is also the separation that - # decides ambiguity, because a latitude-circle segment travels only in longitude. - dlon = _wrapped_dlon(lon1, lon2) - zonal = abs(lat2 - lat1) <= SEGMENT_ATOL_DEG - meridional = abs(dlon) <= SEGMENT_ATOL_DEG - if not (zonal or meridional): - raise ValueError( - f"Latitude-circle segment from (lon={lon1}, lat={lat1}) to " - f"(lon={lon2}, lat={lat2}) is oblique: its endpoints differ in both " - "latitude and longitude, so no circle of constant latitude passes " - "through them. Use curve='great circle' for oblique segments, or split " - "this one into a zonal leg (same latitude) and a meridional leg (same " - "longitude)." - ) - # A meridional segment runs along a single meridian, so `dlon` is 0 and the - # ambiguity test below is trivially passed: there is no east/west choice to make. - sep = abs(dlon) + segment_curve = _segment_curve(lat1, lat2, curve) + + if curve == "latitude circle" and not _is_constant_latitude(lat1, lat2): + raise ValueError( + f"Segment from (lon={lon1}, lat={lat1}) to (lon={lon2}, lat={lat2}) does not " + f"follow a circle of constant latitude: its endpoints differ in latitude by " + f"{abs(lat2 - lat1)} degrees. Use curve='latitude and great circle' to follow " + "the parallel where the endpoints do share a latitude and the geodesic " + "everywhere else, or curve='great circle' throughout." + ) + + if segment_curve == "latitude circle": + # A segment along a parallel travels only in longitude, so the separation that + # decides ambiguity is the shortest-way-round longitude change -- degrees of + # longitude, which near the poles is a far larger number than the arc it spans. + sep = abs(_wrapped_dlon(lon1, lon2)) + measure = "degrees of longitude along the parallel" else: sep = np.rad2deg(distance_on_unit_sphere(lon1, lat1, lon2, lat2, R=1.)) + measure = "degrees of arc" - if sep >= 180. - SEGMENT_ATOL_DEG: + if sep >= 180. - HALF_CIRCLE_ATOL_DEG: raise ValueError( f"Segment from (lon={lon1}, lat={lat1}) to (lon={lon2}, lat={lat2}) has " - f"endpoints half a circle apart (separation ~{sep:.4f} degrees), so the two " - "ways round are equally short and the shortest path between them is not " - "unique. Add an intermediate waypoint to say which way the section goes." + f"endpoints half a circle apart ({sep:.4f} {measure}), so neither way round " + "is the shorter and there is no shortest path to take. Add an intermediate " + "waypoint to say which way the section goes." ) + return segment_curve + def infer_grid_path_from_geo(lonstart, latstart, lonend, latend, gridlon, gridlat, neighbor_maps, curve="great circle"): """ @@ -572,7 +610,12 @@ def infer_grid_path_from_geo(lonstart, latstart, lonend, latend, gridlon, gridla (lons_c, lats_c) are the corresponding longitude and latitudes. """ - _check_segment_span(lonstart, latstart, lonend, latend, curve) + # This function is called once per segment, so it is where the section-wide `curve` + # request becomes the one curve this segment follows. Resolve it here, from the + # *requested* waypoints, rather than leaving it to `infer_grid_path`, which sees only + # the grid corners the waypoints snap to -- and whose latitudes can differ from the + # requested ones by up to half a cell. + segment_curve = _check_segment_span(lonstart, latstart, lonend, latend, curve) multitile = np.ndim(gridlon) == 3 if multitile: @@ -591,7 +634,7 @@ def infer_grid_path_from_geo(lonstart, latstart, lonend, latend, gridlon, gridla gridlon, gridlat, neighbor_maps=neighbor_maps, - curve=curve, + curve=segment_curve, f1=fstart, f2=fend, ) @@ -625,6 +668,10 @@ def infer_grid_path(i1, j1, i2, j2, gridlon, gridlat, neighbor_maps, f1=None, f2 them -- typically via `sectionate.grid_section`. f1, f2: integer or None Face indices of the starting and ending points (multi-tile grids only); None otherwise. + curve: str + Curve this segment follows; see `sectionate.grid_section`. "latitude and great + circle" resolves here from the two endpoint corners' latitudes, since this entry + point is given indices rather than requested waypoints. RETURNS: ------- @@ -669,6 +716,11 @@ def neighbor(direction, f, j, i): lon1, lat1 = coord(gridlon, f1, j1, i1), coord(gridlat, f1, j1, i1) lon2, lat2 = coord(gridlon, f2, j2, i2), coord(gridlat, f2, j2, i2) + # Which curve this segment follows. `infer_grid_path_from_geo` has normally resolved + # it already, from the requested waypoints; resolving again here is a no-op for it and + # makes all three `CURVES` usable when this function is called directly with indices. + segment_curve = _segment_curve(lat1, lat2, curve) + # Per-curve metrics used by the deterministic neighbor selection below. # - progress(lon, lat): remaining distance to the segment endpoint (smaller = nearer); # admits only neighbors that do not move away from the endpoint. @@ -678,25 +730,16 @@ def neighbor(direction, f, j, i): # for both curve types so WALK_DEVIATION_ATOL is meaningful for both. # Physical coincidence with the endpoint (the seam-twin stop) always uses true geodesic # distance, independent of `curve`. - if curve == "great circle": + if segment_curve == "great circle": def progress(lon, lat): return distance_on_unit_sphere(lon, lat, lon2, lat2) def deviation(lon, lat): return (spherical_angle(lon2, lat2, lon1, lat1, lon, lat) + spherical_angle(lon1, lat1, lon2, lat2, lon, lat)) - elif curve == "latitude circle": - # These two metrics measure progress purely in longitude and deviation purely in - # latitude, which is exactly right for the two segment shapes `_check_segment_span` - # admits for this curve, and for nothing else: - # - a *zonal* segment (endpoints on one parallel) marches in longitude at constant - # latitude, which is what `progress` drives and `deviation` holds it to; - # - a *meridional* segment (endpoints on one meridian) has lon2 == lon1, so - # `progress` is flat and admits no neighbor. The walk then falls through to the - # "no admissible forward move" branch below, which steps to the neighbor - # geodesically closest to the endpoint -- i.e. straight down the meridian. - # An oblique segment would instead walk zonally along the *start* latitude and only - # then close the latitude gap, an L-shaped path that depends on which end you start - # from; `_check_segment_span` rejects those up front rather than tracing one. + else: # "latitude circle" -- the only other value `_segment_curve` returns + # Progress purely in longitude, deviation purely in latitude: the metrics of a + # march along a parallel, and meaningful only for a segment whose endpoints share + # a latitude. Segments that do not are never resolved to this curve. def progress(lon, lat): # sin^2(delta-lon/2) is the haversine of the longitude gap: periodic in 360 # degrees and monotonic in |delta-lon| up to 180, so it measures the *shortest* @@ -708,10 +751,6 @@ def deviation(lon, lat): # (radians). Flat between the two endpoint latitudes, which absorbs the # sub-cell mismatch left when each endpoint snaps to its nearest grid corner. return np.deg2rad(abs(lat - lat1)) + np.deg2rad(abs(lat - lat2)) - else: - raise ValueError( - f"curve must be 'great circle' or 'latitude circle'; got {curve!r}." - ) def order_key(pt): _f, _j, _i = pt diff --git a/sectionate/tests/test_ecco_llc90.py b/sectionate/tests/test_ecco_llc90.py index 7a08a25..27a764c 100644 --- a/sectionate/tests/test_ecco_llc90.py +++ b/sectionate/tests/test_ecco_llc90.py @@ -97,3 +97,46 @@ def test_llc90_face_corners_resolve_topologically(): assert fm.shape == jm.shape == im.shape assert (fm >= 0).all() and (fm < nf).all() assert (jm >= 0).all() and (im >= 0).all() + + +def test_llc90_meridional_segment_traces_both_directions(): + """A meridional segment on a real multi-tile grid, under + ``curve="latitude and great circle"``. Its endpoints do not share a latitude, so the + combined option routes it to the geodesic -- which is what makes it traceable at all: + the constant-latitude metrics measure progress purely in longitude, so along a + meridian they are flat, the walk never converges, and it eventually gives up with + "Should have reached the endpoint by now." Being direction-independent, it must also + trace identically whichever end it starts from.""" + from sectionate.section import grid_section + + grid = _load_grid() + fwd = grid_section(grid, [0., 0.], [60., 80.], curve="latitude and great circle") + rev = grid_section(grid, [0., 0.], [80., 60.], curve="latitude and great circle") + + i, j, f, lons, lats = fwd + assert len(i) == 49 + assert lats[0] < lats[-1] # it really does head north + for a, b in zip(fwd, rev): + assert np.array_equal(a, b[::-1]) + + # curve="latitude circle" refuses the same segment outright rather than walking it. + with pytest.raises(ValueError, match="constant latitude"): + grid_section(grid, [0., 0.], [60., 80.], curve="latitude circle") + + +def test_llc90_zonal_segment_holds_its_parallel(): + """The complement of the test above: a segment whose endpoints do share a latitude is + routed to the parallel by the combined option, giving the same path as an explicit + ``curve="latitude circle"``. The grid stores its corner latitudes in float32, so this + also exercises the classification tolerance on real single-precision coordinates.""" + from sectionate.section import grid_section + + grid = _load_grid() + combined = grid_section(grid, [0., 60.], [20., 20.], + curve="latitude and great circle") + parallel = grid_section(grid, [0., 60.], [20., 20.], curve="latitude circle") + for a, b in zip(combined, parallel): + assert np.array_equal(a, b) + + lats = combined[4] + assert np.ptp(lats) < 2. # stays within a cell of 20N diff --git a/sectionate/tests/test_section_cornercases.py b/sectionate/tests/test_section_cornercases.py index 2ab51e2..e3e658c 100644 --- a/sectionate/tests/test_section_cornercases.py +++ b/sectionate/tests/test_section_cornercases.py @@ -115,9 +115,9 @@ def test_latitude_circle_takes_shortest_path_west(): def test_half_circle_separation_is_the_only_ambiguity_error(): - """The single remaining ambiguity error covers both curve types: endpoints exactly - half a circle apart have two equally short paths between them, so the shortest path - is not unique. Anything short of that resolves silently to the shortest path.""" + """The single remaining ambiguity error covers every curve: endpoints written exactly + half a circle apart have two equally long ways round, so neither is the shorter. + Anything short of that resolves silently to the shortest path.""" from sectionate.section import grid_section, _check_segment_span ambiguous = [ @@ -127,14 +127,17 @@ def test_half_circle_separation_is_the_only_ambiguity_error(): (10., 30., 190., 30., "latitude circle"), # 180 degrees away from lon=10 (0., 0., 180., 0., "great circle"), # antipodal on the equator (0., 30., 180., -30., "great circle"), # antipodal off the equator + (0., 0., 180., 0., "latitude and great circle"), # resolves to the parallel + (0., 30., 180., -30., "latitude and great circle"), # resolves to the geodesic ] for lon1, lat1, lon2, lat2, curve in ambiguous: with pytest.raises(ValueError, match="half a circle apart"): _check_segment_span(lon1, lat1, lon2, lat2, curve) - # Just short of half a circle is unambiguous for both curves. + # Just short of half a circle is unambiguous for every curve. _check_segment_span(0., 0., 179.9, 0., "latitude circle") _check_segment_span(0., 0., 179.9, 0., "great circle") + _check_segment_span(0., 0., 179.9, 0., "latitude and great circle") # The same single error reaches the user through the public API. with pytest.raises(ValueError, match="half a circle apart"): @@ -143,13 +146,54 @@ def test_half_circle_separation_is_the_only_ambiguity_error(): grid_section(grid, [0., 180.], [0., 0.]) -def test_latitude_circle_meridional_segment_walks_the_meridian(): - """A latitude-circle segment whose endpoints share a longitude is the meridional - connector used to join zonal arcs at different latitudes. It is unambiguous, so it - is allowed, and it walks straight down the meridian.""" - from sectionate.section import grid_section +def test_half_circle_error_names_the_measure_it_reports(): + """The separation the error quotes is degrees of LONGITUDE for a segment along a + parallel and degrees of ARC for a geodesic one -- two different quantities that + coincide only on the equator. At 89N, 180 degrees of longitude is a 222 km hop, so + the number is only meaningful if the message says which measure it is.""" + from sectionate.section import _check_segment_span + + with pytest.raises(ValueError, match="degrees of longitude along the parallel"): + _check_segment_span(0., 89., 180., 89., "latitude circle") + with pytest.raises(ValueError, match="degrees of arc"): + _check_segment_span(0., 89., 180., -89., "great circle") + + +def test_latitude_circle_requires_constant_latitude(): + """Under curve="latitude circle" every segment must lie on a circle of constant + latitude. Meridional and oblique segments alike lie on none, so both raise.""" + from sectionate.section import grid_section, _check_segment_span + + not_constant_latitude = [ + (60., -80., 60., 80.), # meridional + (0., 0., 120., 40.), # oblique + ] + for lon1, lat1, lon2, lat2 in not_constant_latitude: + with pytest.raises(ValueError, match="constant latitude"): + _check_segment_span(lon1, lat1, lon2, lat2, "latitude circle") + + with pytest.raises(ValueError, match="constant latitude"): + grid_section(grid, [60., 60.], [-80., 80.], curve="latitude circle") + with pytest.raises(ValueError, match="constant latitude"): + grid_section(grid, [0., 120.], [0., 40.], curve="latitude circle") + # Rejected even as one leg of a section whose other legs are constant-latitude. + with pytest.raises(ValueError, match="constant latitude"): + grid_section(grid, [0., 120., 120.], [0., 0., 40.], curve="latitude circle") - i, j, lons, lats = grid_section(grid, [60., 60.], [-80., 80.], curve="latitude circle") + +def test_latitude_and_great_circle_chooses_per_segment(): + """curve="latitude and great circle" decides segment by segment: constant-latitude + segments follow the parallel, everything else follows the geodesic. Segments that + "latitude circle" rejects are therefore traced, not raised on.""" + from sectionate.section import grid_section, _check_segment_span + + combined = "latitude and great circle" + assert _check_segment_span(0., 40., 120., 40., combined) == "latitude circle" + assert _check_segment_span(60., -80., 60., 80., combined) == "great circle" + assert _check_segment_span(0., 0., 120., 40., combined) == "great circle" + + # A meridional segment walks straight down the meridian. + i, j, lons, lats = grid_section(grid, [60., 60.], [-80., 80.], curve=combined) assert np.all([ modequal(i, np.array([1, 1, 1, 1, 1])), modequal(j, np.array([0, 1, 2, 3, 4])), @@ -157,19 +201,100 @@ def test_latitude_circle_meridional_segment_walks_the_meridian(): modequal(lats, np.array([-80., -40., 0., 40., 80.])), ]) + # An oblique segment traces exactly the great-circle path. + oblique = grid_section(grid, [0., 120.], [0., 40.], curve=combined) + assert np.all([np.array_equal(a, b) + for a, b in zip(oblique, grid_section(grid, [0., 120.], [0., 40.]))]) -def test_latitude_circle_oblique_segment_rejected(): - """No circle of constant latitude passes through two points that differ in both - latitude and longitude, so such a segment is rejected rather than traced as an - L-shaped, direction-dependent path. Great circle handles it instead.""" - from sectionate.section import grid_section, _check_segment_span + # A zonal segment traces exactly the latitude-circle path. + zonal = grid_section(grid, [0., 120.], [0., 0.], curve=combined) + assert np.all([ + np.array_equal(a, b) for a, b in + zip(zonal, grid_section(grid, [0., 120.], [0., 0.], curve="latitude circle")) + ]) - with pytest.raises(ValueError, match="oblique"): - _check_segment_span(0., 0., 120., 40., "latitude circle") - with pytest.raises(ValueError, match="oblique"): - grid_section(grid, [0., 120.], [0., 40.], curve="latitude circle") + # And a section mixing the two kinds of leg is traced end to end. + i, j, lons, lats = grid_section(grid, [0., 120., 120.], [0., 0., 40.], curve=combined) + assert np.all([ + modequal(lons, np.array([0., 60., 120., 120.])), + modequal(lats, np.array([0., 0., 0., 40.])), + ]) + + +def test_unknown_curve_is_rejected(): + """An unrecognized `curve` names all three accepted values.""" + from sectionate.section import grid_section + + with pytest.raises(ValueError, match="curve must be one of"): + grid_section(grid, [0., 120.], [0., 0.], curve="rhumb line") + + +def _fine_global_grid(dlon=2., dlat=2.): + """A 2-degree global 'outer' grid, fine enough that a parallel and a geodesic between + the same two waypoints are visibly different paths. The 7x5 fixture above cannot show + this: its corner rows are 40 degrees apart, so both curves snap to the same one.""" + xq = np.arange(0., 360. + dlon, dlon) + yq = np.arange(-90., 90. + dlat, dlat) + lon_f, lat_f = np.meshgrid(xq, yq) + ds_f = xr.Dataset({}, coords={ + "xq": xr.DataArray(xq, dims=("xq",)), + "yq": xr.DataArray(yq, dims=("yq",)), + "lon_c": xr.DataArray(lon_f, dims=("yq", "xq")), + "lat_c": xr.DataArray(lat_f, dims=("yq", "xq")), + }) + return xgcm.Grid( + ds_f, + coords={"X": {"outer": "xq"}, "Y": {"outer": "yq"}}, + padding={"X": "periodic", "Y": "extend"}, + autoparse_metadata=False, + ) + + +def test_parallel_is_held_where_the_geodesic_bows(): + """The test that actually separates the two metrics. Between (0, 40N) and (120, 40N) + the geodesic bows a long way poleward -- it is the shorter path -- while the parallel + does not. On a 2-degree grid the constant-latitude walk holds lat 40 for all 61 of its + points; the great-circle walk climbs to 60N and takes 83.""" + from sectionate.section import grid_section + + fine = _fine_global_grid() + + i, j, lons, lats = grid_section(fine, [0., 120.], [40., 40.], curve="latitude circle") + assert np.all(lats == 40.) + assert len(i) == 61 + assert np.array_equal(lons, np.arange(0., 122., 2.)) + + i_gc, j_gc, lons_gc, lats_gc = grid_section(fine, [0., 120.], [40., 40.]) + assert lats_gc.max() == 60. # bows ~20 degrees poleward + assert len(i_gc) == 83 + + # The combined option classifies this segment as constant-latitude, so it must + # reproduce the parallel exactly, not the geodesic. + combined = grid_section(fine, [0., 120.], [40., 40.], + curve="latitude and great circle") + assert np.all([np.array_equal(a, b) for a, b in zip(combined, (i, j, lons, lats))]) + + +def test_float32_waypoints_still_count_as_constant_latitude(): + """Model corner coordinates are routinely stored as float32 (ECCO's are), so + latitudes that a user reads off a grid come back a few 1e-7 degrees -- a few cm -- + from the value they were written as. That must not turn a zonal segment into an + oblique one, which is what a nanodegree classification tolerance would do.""" + from sectionate.section import ( + grid_section, _is_constant_latitude, _check_segment_span, + CONSTANT_LATITUDE_ATOL_DEG, + ) + + lat32 = float(np.float32(26.4)) + offset = abs(lat32 - 26.4) + assert 1.e-9 < offset < CONSTANT_LATITUDE_ATOL_DEG # ~4.e-7 degrees, ~4 cm + assert _is_constant_latitude(lat32, 26.4) + assert _check_segment_span(0., lat32, 120., 26.4, "latitude circle") == "latitude circle" + assert _check_segment_span(0., lat32, 120., 26.4, + "latitude and great circle") == "latitude circle" - # The same waypoints are fine as a great circle, and as a zonal leg plus a - # meridional leg under latitude circle. - grid_section(grid, [0., 120.], [0., 40.]) - grid_section(grid, [0., 120., 120.], [0., 0., 40.], curve="latitude circle") + # End to end: both waypoints snap to the same corner row and the path holds it. + fine = _fine_global_grid() + i, j, lons, lats = grid_section(fine, [0., 120.], [lat32, 26.4], curve="latitude circle") + assert np.all(lats == lats[0]) + assert len(i) == 61 diff --git a/sectionate/tests/test_section_reversibility.py b/sectionate/tests/test_section_reversibility.py index 49e841d..8be57c7 100644 --- a/sectionate/tests/test_section_reversibility.py +++ b/sectionate/tests/test_section_reversibility.py @@ -82,10 +82,16 @@ def test_steep_and_shallow_oblique_sections_reversible(): np.array([-7.0, -5.0, -3.0, -1.0])) # shallow -def test_oblique_latitude_circle_segments_reversible(): +def test_zonal_and_meridional_segments_reversible(): grid = _fine_grid(seed=5) - # Latitude-circle segments (constant-latitude hops at two different latitudes, - # connected by a meridional step) must also reverse exactly. + # Constant-latitude hops at two different latitudes, joined by a meridional step. + # Under "latitude and great circle" the two zonal legs follow their parallels and the + # meridional leg follows the geodesic; both metrics must reverse exactly, and so must + # the per-segment choice between them. lon = np.array([2.0, 10.0, 10.0, 20.0]) lat = np.array([-6.0, -6.0, 4.0, 4.0]) - _assert_reverse_equal(grid, lon, lat, curve="latitude circle") + _assert_reverse_equal(grid, lon, lat, curve="latitude and great circle") + + # A purely zonal section reverses exactly under "latitude circle" too. + _assert_reverse_equal(grid, np.array([2.0, 10.0, 20.0]), np.array([-6.0, -6.0, -6.0]), + curve="latitude circle") From 124fa06f09d3121ec95b614b657e532b51ca361c Mon Sep 17 00:00:00 2001 From: Henri Drake Date: Tue, 11 Aug 2026 13:13:28 -0700 Subject: [PATCH 3/4] Update section.py --- sectionate/section.py | 17 ++++------------- 1 file changed, 4 insertions(+), 13 deletions(-) diff --git a/sectionate/section.py b/sectionate/section.py index 238d5b4..2dcc707 100644 --- a/sectionate/section.py +++ b/sectionate/section.py @@ -289,27 +289,18 @@ def grid_section(grid, lons, lats, curve="great circle"): - "great circle" (default): every segment follows the geodesic. - "latitude circle": every segment follows a circle of constant latitude, marching in longitude. A segment whose endpoints do not share a latitude (to within - 1.e-6 degrees) lies on no such circle, so it raises a ValueError -- meridional - segments included. + 1.e-6 degrees) lies on no such circle, so it raises a ValueError. - "latitude and great circle": decided per segment. A segment whose endpoints share a latitude follows the parallel; every other segment follows the geodesic. - This is the option for a section that is zonal in places and joined up by - meridional or slanted legs elsewhere. + This is the option for a section that is zonal in places and joined up by arbitrarily + oriented legs elsewhere. Under every option each segment takes the **shortest** path between its two vertices. Raw longitudes are never read as a request to go the long way round, so - a segment written 0 -> 270 along a parallel runs 90 degrees *west*. Encircle the + a segment written 0 -> 270 along the equator runs 90 degrees *west*. Encircle the globe by giving intermediate vertices (e.g. 0 -> 120 -> 240 -> 360), which is also what says which way round it goes. - A ValueError is also raised for a segment whose endpoints are written exactly half - a circle apart -- antipodal for a geodesic segment, or exactly +/-180 degrees of - longitude apart for one along a parallel -- because then neither way round is the - shorter. That check only catches endpoints typed as exactly half a circle apart: - a hair short of it, which way the segment runs is decided by the sign of the - round-off, so subdivide near-half-circle segments with an intermediate vertex - rather than expecting to be warned about them. - Returns ------- i_c, j_c[, f_c], lons_c, lats_c: `np.ndarray` From 007a5340692fb0d5c1fc862611c094cb82e6ca9d Mon Sep 17 00:00:00 2001 From: Henri Drake Date: Tue, 11 Aug 2026 13:16:25 -0700 Subject: [PATCH 4/4] Update section.py --- sectionate/section.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/sectionate/section.py b/sectionate/section.py index 2dcc707..dcbffab 100644 --- a/sectionate/section.py +++ b/sectionate/section.py @@ -707,9 +707,6 @@ def neighbor(direction, f, j, i): lon1, lat1 = coord(gridlon, f1, j1, i1), coord(gridlat, f1, j1, i1) lon2, lat2 = coord(gridlon, f2, j2, i2), coord(gridlat, f2, j2, i2) - # Which curve this segment follows. `infer_grid_path_from_geo` has normally resolved - # it already, from the requested waypoints; resolving again here is a no-op for it and - # makes all three `CURVES` usable when this function is called directly with indices. segment_curve = _segment_curve(lat1, lat2, curve) # Per-curve metrics used by the deterministic neighbor selection below.