Port simple_paleobathymetry (Steps 1-4) and its distance/sediment engine into gplately

[michaelchin]

Summary

Implements #444 and #445 together (they turned out to be tightly coupled - see the discussion on #444).

gplately.paleobathymetry (#444) - Steps 1, 3, 4 of EarthByte's simple_paleobathymetry workflow, plus an orchestrator:

• age_to_basement_depth() - Step 1: seafloor age -> basement depth (gdh1/rhcw18/parsons_sclater/crosby09). Ships the RHCW18 lookup table.
• dutkiewicz_2017_sediment_thickness() - Step 3: Dutkiewicz et al. (2017) sediment-thickness polynomial fit.
• sediment_isostatic_correction() / paleobathymetry() - Step 4: Sykes (1996) isostatic correction, combined into final paleobathymetry.
• simple_paleobathymetry() - runs Steps 1-4 end to end, calling into gplately.sediment_thickness for Step 2.

gplately.sediment_thickness (#445) - a port of predicting-sediment-thickness's ocean_basin_proximity.py engine:

• generate_distance_grids() - Step 2: reconstructs each ocean point backward through time and computes its lifetime-mean distance to the nearest proximity feature (e.g. passive-margin COB line segments), using pygplates.TopologicalModel and gplately's own ptt.utils.proximity_query.
• generate_sediment_thickness_grids() - Step 3's grid-level driver, composing the above with dutkiewicz_2017_sediment_thickness().
• Exposed on the CLI: gplately generate-distance-grids / generate-sediment-grids (aliases gdg/gsg).

Validation

• Steps 1, 3, 4 numerically cross-checked against the original workflow's own reference code (diffs at floating-point noise, ~1e-12).
• Step 2 + the full simple_paleobathymetry() pipeline run end-to-end against a real plate model (Müller et al. 2025, via the Plate Model Manager) with a downloaded age grid: produces physically sensible distance-to-margin (0-3000 km) and paleobathymetry (-2 to -5.5 km) values.
• CLI (generate-distance-grids -> generate-sediment-grids) smoke-tested end-to-end against the same real data.
• Unit tests added for all of the above (test_9_paleobathymetry.py, test_10_sediment_thickness.py), using the existing Muller2019 test fixtures.

Deliberately out of scope

• Continent-obstacle routing in Step 2 (shortest path around continents, via the original's shortest_path.py) - distances are currently always great-circle. This is a real difference from simple_paleobathymetry's default config (use_continent_obstacles: true).
• Topological proximity features - only static/non-topological proximity features (e.g. COB line segments) are supported.
• Step 5 (pyBacktrack merge) - tracked separately in Integrate pyBacktrack: merge present-day and subducted paleobathymetry grids into complete grids #447.

All flagged in the relevant module docstrings.

Test plan

• python -m pytest -vv tests-dir/pytestcases/test_9_paleobathymetry.py tests-dir/pytestcases/test_10_sediment_thickness.py tests-dir/pytestcases/test_0_imports.py - all passing.
• Manual cross-check of Steps 1/3/4 against reference implementation.
• Manual end-to-end run (Python API and CLI) against a real plate model, sanity-checked output ranges.
• black --check clean on all changed/added files.
• Full pytest -vv tests-dir/pytestcases (not run here - some fixtures need network access / take a while).

🤖 Generated with Claude Code

[michaelchin]

Added a `gplately paleobathymetry` CLI command (alias `pb`), wrapping `simple_paleobathymetry()` to run Steps 1-4 end to end in one command - previously there was CLI coverage for Steps 2/3 only (`generate-distance-grids`/`generate-sediment-grids`), with no CLI path to the full pipeline even though that's the actual one-command UX the source `simple_paleobathymetry` workflow offers.

Refactored the shared argparse/file-resolution logic out of `commands/sediment_thickness.py` so it's reused by the new command rather than duplicated. Verified end-to-end against a real plate model via the CLI, and added the same invocation to `tests-dir/test-cli.sh`.

[michaelchin]

Added Step 5 (pyBacktrack merge, #447) too, so this PR now covers #444, #445 and #447.

• `gplately.grids.pybacktrack_paleobathymetry.merge_pybacktrack_paleobathymetry()` - thin wrapper around `pybacktrack.reconstruct_paleo_bathymetry_grids()`'s merge support. `pybacktrack` is an optional dependency (new `gplately[paleobathymetry]` extra), imported lazily so `import gplately` never requires it.
• Wired into `simple_paleobathymetry()` as an opt-in step (`pybacktrack=True`) and into the CLI as `gplately paleobathymetry --pybacktrack` (`--static-polygons`/`--present-day-age-grid` auto-resolve from `-m/--model` when omitted).
• Actually validated with the real `pybacktrack` package (I had it install cleanly in the test env): ran the full pipeline end-to-end (Python API and CLI) against `muller2025`, and the merged output goes from 330 to 441 finite cells (out of 703, at a coarse 10° test grid) - i.e. pyBacktrack is correctly filling in the submerged-continental-crust / already-subducted-crust gaps that Steps 1-4 leave as NaN, which is the whole point of Step 5.
• Unit tests added, skipped automatically via `pytest.importorskip` when `pybacktrack` isn't installed (kept it out of the core test env since it's a large optional dependency).

Only thing still open for the whole paleobathymetry effort now is #446 (continent-contouring passive-margins use case) and the continent-obstacle routing gap noted in #445's module docstring.

[michaelchin]

Added #446 too (continent-contouring passive margins), so this PR now covers #444, #445, #446 and #447 - the whole paleobathymetry tracking effort (#448) except the continent-obstacle routing gap noted in #445's docstring.

• `gplately.grids.continent_contouring` - port of `continent-contouring`'s `create_passive_margins.py`, built on gplately's own already-vendored `ptt.continent_contours.ContinentContouring` engine. `passive_margin_polylines()` (the core segment-splitting algorithm) + `generate_passive_margins()` (the full time-stepping driver).
• CLI: `gplately generate-passive-margins` (alias `gpm`).
• Validated the actual point of this: `generate_passive_margins()`'s `passive_margin_features` output composes directly with Move Predicting Sediment Thickness workflow into GPlately, expose via API/CLI #445's `generate_distance_grids()` as a dynamically-contoured alternative to a static COB file - zero glue code needed (it's already the right `pygplates.FeaturesFunctionArgument`-compatible type). Ran this composition end-to-end against `muller2025` and confirmed it produces sensible distance grids.
• Unit tests added, including hand-verified geometric cases for the segment-splitting logic (e.g. a contour with an active-margin gap in the middle correctly splits into two separate passive-margin polylines).

[michaelchin]

Added continent-obstacle routing too, closing the one remaining documented gap from #445.

• `gplately.lib.shortest_path` - a faithful port of `predicting-sediment-thickness`'s `shortest_path.py`: a spherical grid + Dijkstra's algorithm that computes shortest-path distances around obstacle geometries instead of a great-circle line that can cut through land.
• Wired into `generate_distance_grids()` as new optional args (`continent_obstacle_features`, `plate_boundary_obstacle_feature_types`, `shortest_path_grid_subdivision_depth`), and propagated through to `simple_paleobathymetry()` and the CLI (`--route-around-continents`/`--continent-obstacles`/`--shortest-path-grid-depth`, on both `generate-distance-grids` and `paleobathymetry`).
• Validated: a clean synthetic-obstacle unit test confirms routed distance > great-circle when blocked, ≈ great-circle when clear. Against real data (`muller2025` + real coastlines), obstacle-routed distances average ~70-80km longer than great-circle across a test grid, with a small number of individual-point exceptions that shrink as grid resolution increases - consistent with the algorithm's own known grid-interpolation smoothing behaviour, not a port defect.

With this, the only thing left from the original scope (per #448) is #446... which is also already in this PR. So as far as I can tell, everything from the original email thread is now implemented, at least at a first-pass level of completeness.

[michaelchin]

Ran a thorough multi-angle code review over the whole diff before asking for review. Fixed everything that was a genuine bug:

• `simple_paleobathymetry()` crashed with `TypeError: got multiple values for keyword argument 'max_distance_km'` if `sediment_thickness_kwargs` used its own documented override of that name.
• `age_to_basement_depth(model="rhcw18")` mapped tiny negative ages (common floating-point noise right at a ridge) to 0 m (sea level) instead of ridge-crest depth (~-2500 m) - the other three models already special-cased this, rhcw18 didn't.
• `simple_paleobathymetry(pybacktrack=True)` now fails with a clear error instead of an obscure one if `rotation_model` is passed as an already-built `pygplates.RotationModel` (pyBacktrack needs raw filenames).
• `passive_margin_polylines()` could incorrectly split a single passive-margin stretch into two output polylines when it straddled a closed contour ring's arbitrary start/end point.
• `--time-step`/`--min-time`/`--max-time` (declared `type=float`) were silently truncated via `int()` before building the time range - `--time-step 0.5` crashed outright, `--time-step 2.5` silently became 2. Fixed with a shared float-aware helper.
• A plate model missing a `Topologies` layer raised an opaque plate-model-manager exception instead of the intended friendly error message.
• `generate-passive-margins` didn't merge in a model's `Cratons` layer the way `seafloor_grids.py`'s equivalent resolution does.
• `--pybacktrack` silently dropped all but the first file when a model's `StaticPolygons` layer spans multiple files - now merges them instead.
• `tests-dir/test-cli.sh` referenced `test_seafloor_gridding.sh` (underscores); the real file is `test-seafloor-gridding.sh` (hyphens) - would have aborted the CLI smoke test.

Added regression tests for the algorithm/API-level fixes; full suite (33 tests) passes, and I re-ran the CLI commands manually including the previously-crashing fractional `--time-step` case.

Deliberately not fixed, flagged for awareness rather than action in this PR:

• A real performance issue in `generate_distance_grids()`: when continent-obstacle routing is on and multiple age grids are processed together, each age grid's backward-reconstruction trajectory independently rebuilds the obstacle grid and reruns Dijkstra at every overlapping time step, so total work scales roughly O(N²) in the number of age grids instead of O(N) (this is exactly what the original `ocean_basin_proximity.py` was structured to avoid, by design, and is a legitimate restructuring, not a quick fix - didn't want to rush it in alongside everything else here).
• A batch of "duplicated code" / "could reuse a shared helper" / "could use scipy instead of hand-rolled Dijkstra" style findings across several files - real observations, but refactors rather than bugs, and `gplately/lib/shortest_path.py` in particular is a deliberately faithful port where I'd rather not deviate from the validated original.
• Two findings (the `math.pi` fallback being counted into the distance mean, and `DistanceGrid.shortest_distance()`'s point-vs-other-geometry inconsistency) turned out to be faithful copies of the original algorithm's own behaviour, not defects introduced here - confirmed against the source I ported from.

I think this is ready for human review now.

[michaelchin]

@jcannon-gplates

Most of the work in the PR was done by Claude Code.

[jcannon-gplates]

Wow, this is quite some work - nice!

I'll take a look - it might be a while though.