diff --git a/CLAUDE.md b/CLAUDE.md index 3d90752..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"`); 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. 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 b730346..dcbffab 100644 --- a/sectionate/section.py +++ b/sectionate/section.py @@ -17,6 +17,28 @@ # index, so the path is independent of platform and of travel direction. WALK_DEVIATION_ATOL = 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""" def __init__(self, name, coords, children = {}, parent = None): @@ -262,9 +284,22 @@ 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 must span - less than 180 degrees, otherwise the direction is ambiguous and a ValueError is raised. + 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. + - "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 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 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. Returns ------- @@ -351,6 +386,11 @@ 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), + "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: ------- @@ -452,43 +492,79 @@ 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 _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 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. + """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": - 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.: - 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." - ) + 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.)) - 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." - ) + measure = "degrees of arc" + + 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 ({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"): @@ -525,7 +601,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: @@ -544,7 +625,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, ) @@ -578,6 +659,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: ------- @@ -622,6 +707,8 @@ 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) + 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. @@ -631,23 +718,27 @@ 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": + 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): - # 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( - 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 8299ba9..e3e658c 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,213 @@ 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.""" +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 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 = [ + # (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 + (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 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"): + 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_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 - 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") + + 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") + + +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])), + modequal(lons, np.array([60., 60., 60., 60., 60.])), + 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.]))]) + + # 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")) + ]) + + # 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" + + # 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")