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.

[michaelchin]

Wow, this is quite some work - nice!

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

Dietmar wanted it be included in the next GPlately release. In fact, I had completely forgotten the paleobathymetry task until he reminded me yesterday. It is the reason why I suddenly asked Claude Code to do this after I've tagged v2.1.0rc0. Take your time. Either in or not in the next release is fine with me. I am also hesitate to rush it in.

[michaelchin]

Wow, this is quite some work - nice!

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

By the way, Claude wrote almost all the long comments here. You may want your Claude give you a brief summary of those long comments. I found sometimes it took some time to read the reports from Claude.

[michaelchin]

@copilot resolve the merge conflicts in this pull request

[Copilot]

@copilot resolve the merge conflicts in this pull request

Resolved by merging master into this branch and fixing the conflict in sphinx-doc/source/command_line_interface.rst in commit 2a52a58.

[jcannon-gplates]

Thanks @michaelchin — this is a lot of ground covered, and the parts that are straight ports are genuinely well done.

I reviewed it by diffing against the original source rather than reading the diff, since the question I care about is fidelity to the workflows rather than whether the new code is internally sensible. Worth flagging up front, because it's why most of what follows wasn't visible in your own review pass — several of these need the ancestor to see at all.

How I'd suggest we split this up. I'm happy to take nearly all of it. I wrote most of the original code, so the fixes and the decisions will be quicker for me than for you — and it saves you the tokens, which I know get tight. I'd rather not hand you a long list of things to interpret either:

• The blocking items below — I'll fix them on this branch, unless you'd rather do them yourself. They're all small; the largest is a one-line change repeated in three places.

• Everything under "Decisions rather than fixes" — treat as decided, they're calls about my own original code. I've written each one up so it's on the record, and I'll check with you at the time if I'm unsure about anything.

• The follow-up issues — mine, except the test plumbing. I've drafted them all under tracking issue #457 so you don't have to write any of this up:

  	Issue	Owner
  #458	Blocking fixes for this PR (checklist; closed by #449)	me
  #459	Restore the shared time loop in generate_distance_grids() — the 5–8× win	me
  #460	Restore the internal/output grid-spacing split and upscaling filter	me
  #461	Restore parallelism across Steps 1–3 and Step 5	me
  #462	Contouring fidelity: time-dependent thresholds, buffer ramp, subduction-zone output	me
  #463	Test plumbing: test-cli.sh into CI, pybacktrack in test-env.yml, GPLATELY_TEST_LEVEL gating	you, if you're happy to take it
  #464	Tests for the changed surface (the porting delta)	me
  #465	--max_distance_threshold: help text vs. behaviour	me
  #466	Module placement and naming (physics functions in grids/, the continent_contouring / continent_contours clash)	me
  #467	Example notebook for the workflow	me

  Issue #463 is CI and packaging rather than workflow science, which is much more your area than mine — that's the only reason I've put your name on it. Issue #466 touches GPlately's own layout rather than the ported code, and on issue #467 I can get the notebook started but you may well want to tidy it up to match the style of the others — so say the word if you'd rather own either of those.

What's solid

• lib/shortest_path.py is a genuinely faithful port. A token-level diff against the original shows no behavioural deviations at all — constants, the 16-neighbour construction with dateline wrap and pole clamping, the quad-tree hemisphere bounding, the lazy-deletion Dijkstra heap, the inverse-distance blend in shortest_distance, and both memo caches are intact. Dropping the __main__ demo block was the right call.
• Steps 1/3/4: every constant and branch boundary matches, including the Dutkiewicz polynomial with its standardisation constants, and the Sykes correction.
• RHCW18_age_depth.dat is byte-identical to both my copy and pyBacktrack's bundled one, so Step 1 and Step 5 are genuinely consistent for rhcw18. Nice property to have. (gdh1 too — I checked the port's analytic form against pyBacktrack's and they agree to 0.00 m.)
• Step 2's core maths — seeding, backward stepping, retirement, the unweighted mean, the math.pi fallback, mean-only clamping — matches statement for statement.
• pyBacktrack Step 5 passes every argument with the same name and value as the original, and the lazy import is right.
• Good reuse of ptt/continent_contours.py rather than re-implementing it.

Blocking

1. Age grids with descending latitude are sampled upside-down. Three sites: grids/sediment_thickness.py:409-420 and :515-526, grids/paleobathymetry.py:511-526. Each builds extent=(min(lon), max(lon), min(lat), max(lat)), which throws away the sign sample_grid needs. read_netcdf_grid only flips latitude when max(lon) > 180 (_grids.py:371-376), so a −180..180 north-up grid keeps descending rows, and sample_grid then reads row 0 as −90° when it is +90° (_grids.py:858-860).

Reproduced with the PR's exact sampling call — one synthetic field where age increases with latitude, written twice with the two equally valid on-disk row orders:

query latitude      : [-75. -45. -15.  15.  45.  75.]
expected age (Myr)  : [ 15.  45.  75. 105. 135. 165.]
PR code, lat ASC    : [ 15.  45.  75. 105. 135. 165.]
PR code, lat DESC   : [165. 135. 105.  75.  45.  15.]   <-- exactly mirrored
signed extent, DESC : [ 15.  45.  75. 105. 135. 165.]   <-- _grids.py:390-395 form

max error on the descending grid: 150.0 Myr (field spans 15-165 Myr)

It's silent — no warning, and the output keeps entirely plausible ranges. gplately's own internal callers already use the correct signed form at _grids.py:390-395, so the fix is to use (grid_lon[0], grid_lon[-1], grid_lat[0], grid_lat[-1]), or to pass realign=True. The original workflow was immune because grdtrack and xarray.interp both work off labelled coordinates.

On how likely this is to bite in practice, see the first answer under "Questions I'd flagged" below — short version, GMT- and gplately-written grids are always safe, GDAL-written ones aren't.

2. {:.0f} output filenames collide and silently overwrite. grids/paleobathymetry.py:538, grids/sediment_thickness.py:535, grids/continent_contouring.py:292 — while the distance grids correctly use {:.1f} (sediment_thickness.py:459). The format strings are inherited; what changed is the precondition. The original forced integer times, so {:.0f} was injective, whereas _time_range exists specifically to allow fractional steps. Feeding it the times from --time-step 0.5 over 0→2 Ma:

   0.0 Ma  ->  paleobathymetry_0Ma.nc
   0.5 Ma  ->  paleobathymetry_0Ma.nc     <-- overwrites 0.0
   1.0 Ma  ->  paleobathymetry_1Ma.nc
   1.5 Ma  ->  paleobathymetry_2Ma.nc     <-- banker's rounding
   2.0 Ma  ->  paleobathymetry_2Ma.nc     <-- overwrites 1.5

grids computed : 5     files written : 3     silently lost : 2

So paleobathymetry_0Ma.nc ends up holding the 0.5 Ma grid. The returned dict is correct, so this is disk-only and easy to miss. The distance-grid writer survives the same times unscathed because of its {:.1f} — worth making the other three match it.

3. --time-step is bound to time_increment. commands/sediment_thickness.py:146 and commands/paleobathymetry.py:92. These are two different quantities in the original: time.step chooses the output times, while time_increment=1 is hardwired for the lifetime sampling loop. --time-step 10 samples each point's distance-to-margin every 10 Myr instead of 1 Myr, and degrades the reconstruction trajectory as well. Plausible output, wrong numbers. Needs a separate --time-increment defaulting to 1.

4. The same fault again in Step 5. grids/paleobathymetry.py:584 forwards the reconstruction increment as pyBacktrack's output increment. With times [0, 10, 20] and the default time_increment=1, Step 5 asks to merge paleobathymetry_1Ma.nc through _19Ma.nc, none of which exist — after the full backtrack has already run.

5. COBs are auto-sourced from the plate model, which the original explicitly refuses to do. commands/sediment_thickness.py:82-85 falls back to plate_model.get_layer("COBs"). The original refuses this in three places — the load_plate_model docstring, a capitalised warning in config.yml ("COB POLYGONS MUST NOT BE USED — a polygon outline wraps the whole continent, including active margins, creating false passive-margin proximities"), and the _check_cob_line_segments() advisory. That polygon guard was dropped at the same time, so there is now nothing to catch it. The error message at :91-95 half-knows this already.

6. An all-NaN age grid still writes an output file. sediment_thickness.py:423-427 falls through to the write at :453-461. The original warns and continues. A bad age grid silently yields a valid-looking all-NaN .nc that Step 3 then consumes.

7. No validation of time_increment in the library (the original raises). 0 gives ZeroDivisionError at :148; negative gives an infinite loop. Only the CLI path is guarded.

8. Valid-time epsilon dropped. create_passive_margins.py writes time + 0.5*interval - 1e-4 in three places, with the comment "epsilon to avoid overlap at interval boundaries". continent_contouring.py:275-287 drops it, so adjacent time slices now overlap at every boundary when loaded into GPlates.

9. crosby09 shouldn't be accepted for Step 5. pybacktrack_paleobathymetry.py:48 maps it to pyBacktrack's AGE_TO_DEPTH_MODEL_CROSBY_2007, but those are two different models — the numbers are under "Questions I'd flagged" below. gdh1 and rhcw18 genuinely are consistent across the merge; this one isn't, and it should raise rather than silently produce a grid with a model discontinuity along the vanished-crust boundary.

10. The public API should be locked down before release, not after. Neither generate_distance_grids nor simple_paleobathymetry has a * marker, so all 11 / 16 optional parameters can be passed positionally — including shortest_path_grid_subdivision_depth, distance_threshold_radians and plate_boundary_obstacle_feature_types. That depth was a deliberate speed/quality compromise on my side and is fine to expose, but as an advanced keyword. A bare *, after age_grid_filenames_and_times is one character per signature, and it's the kind of thing that can't be added once people are calling it. Same argument for the 14 new public names being missing from __all__ (so from gplately import * omits the whole new API — the hasattr-based export tests pass and mask this), and for only 4 of the 14 reaching Sphinx.

11. Minor. commands/paleobathymetry.py:69 leaks a NamedTemporaryFile(delete=False). commands/paleobathymetry.py:119 says obstacle routing isn't included, but it is. RHCW18_age_depth.dat is vendored without a licence statement.

Decisions rather than fixes

These are all calls about code I wrote, so I've taken them rather than asking — flagging them here so they're on the record, and happy to be argued out of any of them.

• Restore continent-obstacle routing as the default. config.yml ships use_continent_obstacles: true; the port makes it opt-in. This also means the batching regression under "Performance" is on the default path rather than an edge case.

• Provenance of grids/continent_contouring.py. The docstring says it ports create_passive_margins.py, but it actually descends from simple_paleobathymetry's generate_continent_contours.py after my fix in EarthByte/simple_paleobathymetry#2 — it reproduces wording that exists only in that fix. That's good news: it inherited both corrections (the swapped area/exclusion thresholds, and .inf separation → engine default), so it never had the merge-every-continent-into-one-blob bug. Two consequences worth recording. The docstring should cite the actual ancestor and note that the defaults follow the paleobathymetry-workflow parameter set rather than continent-contouring at main. And it missed the buffer_and_gap_mode: ramp from that same fix — and can't express it, since buffer_and_gap_distance_kms is a float where the engine accepts callables. Deep-time (pre-Pangea) contouring isn't reachable as a result; that part is issue #462 rather than this PR.

• separation_distance_threshold_radians=None is safe, but only by accident of style. It works because :209-213 omits the kwarg so the constructor default applies. The engine treats an explicit None as zero. Worth a comment there so nobody "simplifies" it into a direct pass-through later.

• Keep the closed-ring seam rotation (:97-111), but document it as a deliberate deviation. It emits one polyline where every ancestor emits two. I think it's an improvement, but it changes feature counts and per-feature lengths against the reference passive_margin_features.gpmlz, so anyone comparing will otherwise think something is broken. It also evaluates _near_subduction twice per arc; the memo is in issue #462.

Testing

I agree with the principle that verbatim ports carry the original's track record — shortest_path.py and the exact-constant maths don't need new tests for port fidelity. But that makes the case for concentrating effort on what actually changed, and at the moment the suite doesn't.

Classifying the ~2,200 new code lines: excluding shortest_path.py, only about a quarter is verbatim, grids/*.py is roughly 57% changed, and the numerically risky part is about 90 lines across six grid-I/O call sites — which is exactly where bug 1 lives, three times over. All the confirmed bugs are in changed or new code; none are in verbatim code.

Of the 33 tests, 3 make an exact numerical claim, 19 are sign/shape/NaN properties, and 11 are hasattr/raises/isfile. None use golden values from the original workflow, though reference outputs are committed in the source repos. Concretely: the coefficients, the spacing, and the valid-time epsilon could all change and the suite would stay green. Every distance-grid test passes a single age grid, so the batching behaviour is untestable by construction; and the fixture is written with write_netcdf_grid, which always emits ascending latitude, so bug 1 is invisible to it. test-cli.sh isn't referenced by any workflow and pybacktrack isn't in test-env.yml, so the CLI has no CI coverage and both Step 5 tests always skip.

One nuance on "verbatim needs no tests": it does where the port widened the envelope. generate_input_points_grid's floor(360/spacing) is verbatim, but only ever ran at 0.5 and 0.1 — --grid-spacing 0.3 now yields an actual spacing of 0.30025° while the filename still says 0.3d. Similarly Grid(6) was hardwired and is now a user-settable int, and proximity_geometries can now be empty or very large.

Tests covering the blocking fixes will land with those fixes here. The rest of the changed surface is written up as issue #464, and the CI/packaging side as issue #463.

Performance

Confirmed — and it's a regression introduced by the port rather than something pre-existing. ocean_basin_proximity.py is deliberately structured to batch, and says so in the API docstring, in the task-chunking comment that motivates min_num_age_grids_per_task = 10, and in the CLI help. One outer time loop services all age grids via ocean_basin_reconstructions, with the sharing at :1145-1152. grids/sediment_thickness.py:405 inverts that into a per-age-grid time loop.

Two corrections to the framing. "O(N²)" overstates it: total work is O(M·T) either way; what goes from O(T) to O(M·T) is the shared term, which my own profiling puts at 86–88% of per-time-step cost — so realistically 5× to 7.7×, not unbounded. And with obstacle routing restored as the default, this is the default path rather than an edge case.

Separately, TopologicalModel is built inside the per-time-step loop at :156-158 despite having no time dependence. That one's free.

None of the performance work blocks the merge, on the basis that the port reproduces the original workflow and just runs slower — that's the bar I'm applying. That includes the dropped internal/output grid-spacing split: that split and its upscaling filter were a speed optimisation on my side, not part of the method, so computing directly at output resolution is correct. Two things worth documenting rather than fixing: output computed directly at 0.1° won't be bit-identical to the published product, since the smoothing filter is absent; and the default grid_spacing=0.5 is coarser than the workflow's 0.1° output.

Questions I'd flagged — and the answers

Three things I couldn't settle from the code alone. I've since checked all three, so nothing here needs an answer from you.

Do any age grids in this chain actually store latitude descending? Effectively no, but not never. I swept every grid on my machine — 486 directories' worth of age grids, distance grids, published paleobathymetry and GPlates sample data — and 483 are ascending. Two of the three exceptions are cartopy's bundled sample data. The third is a real EarthByte age grid, and its history attribute says GDAL CreateCopy: GDAL writes north-up, i.e. descending, and its longitudes were outside −180..180 as well, so read_netcdf_grid's realign gate wouldn't have caught it either.

So anything written by GMT or by gplately is always safe, and anything that has been through GDAL, a GeoTIFF or a GIS is not. That downgrades bug 1 from "corrupts the shipped workflow" to "silently corrupts a plausible user-supplied input" — still worth fixing in a public API, and it means the test has to build its descending grid with netCDF4 directly, since write_netcdf_grid can't produce one.

Should crosby09 be rejected for Step 5? Yes. Steps 1–4's crosby09 is the Crosby & McKenzie (2009) empirical piecewise fit; pyBacktrack's CROSBY_2007 is the plate-cooling model from Crosby's 2007 thesis, with explicit mantle density, diffusivity and plate thickness. Different models, and the numbers say so:

 age   crosby09 (Steps 1-4)   CROSBY_2007 (Step 5)      diff
   0            2652.0 m              2601.4 m       +50.6 m
  20            4101.0                4150.7         -49.7
  75            5457.9                5517.8         -59.8
 100            5369.0                5461.8         -92.9
 140            5557.5                5656.8         -99.3
 160            5793.7                5703.5         +90.3

max |difference| over 0-200 Ma : 99.3 m
same check for gdh1            :  0.00 m   (port vs pyBacktrack, exact)

So the merged grid would carry a model discontinuity of up to ~100 m wherever present-day crust meets vanished crust. gdh1 is exact and rhcw18 shares the identical table, so it really is just this one model. That's blocking item 9 — better to raise now than to add a rejection after release.

--max_distance_threshold: help text vs. behaviour. Confirmed, and it's my bug carried over faithfully. The help says over-threshold distances are ignored and "will not contribute to mean / standard deviation"; both code paths set min_distance = math.pi and include it. The reason is that the lookup returns None for "unreachable through obstacles" and for "beyond threshold", and the unreachable case genuinely needs a value, so the two got conflated on the same line. The behaviour I actually want is π·R for unreachable and a real exclusion for over-threshold — but that changes published numbers, so for this PR I'd just make the docstring describe what the code does, and settle the semantics in issue #465 (upstream as well as here).

[michaelchin]

convert back to "Draft" so that the tests don't run for every push.

[jcannon-gplates]

Hi @michaelchin, the blocking items in #458 are complete, so you can mark it ready-for-review if you like.

I did one Ubuntu/3.12 dispatch after the last commit. Marking it will then unlock macOS, Windows and Python 3.10/3.11/3.13.

After this is merged I will work on the remaining PRs listed in #457.