diff --git a/gplately/__init__.py b/gplately/__init__.py index f2b442e9..c450cb53 100644 --- a/gplately/__init__.py +++ b/gplately/__init__.py @@ -52,6 +52,25 @@ from .auxiliary import get_plate_reconstruction, get_gplot from .data_server import DataServer from .raster import Raster +from .grids.paleobathymetry import ( + AGE_DEPTH_MODELS, + DUTKIEWICZ_2017_SEDIMENT_THICKNESS, + age_to_basement_depth, + dutkiewicz_2017_sediment_thickness, + sediment_isostatic_correction, + paleobathymetry, + simple_paleobathymetry, +) +from .grids.sediment_thickness import ( + generate_input_points_grid, + generate_distance_grids, + generate_sediment_thickness_grids, +) +from .grids.pybacktrack_paleobathymetry import merge_pybacktrack_paleobathymetry +from .grids.continent_contouring import ( + passive_margin_polylines, + generate_passive_margins, +) from .grids import ( read_netcdf_grid, write_netcdf_grid, @@ -151,6 +170,20 @@ "load_feature_collection", "get_plate_reconstruction", "get_gplot", + # paleobathymetry workflow (steps 1-5) + "age_to_basement_depth", + "generate_input_points_grid", + "generate_distance_grids", + "generate_passive_margins", + "passive_margin_polylines", + "generate_sediment_thickness_grids", + "dutkiewicz_2017_sediment_thickness", + "sediment_isostatic_correction", + "paleobathymetry", + "simple_paleobathymetry", + "merge_pybacktrack_paleobathymetry", # constants "EARTH_RADIUS", + "AGE_DEPTH_MODELS", + "DUTKIEWICZ_2017_SEDIMENT_THICKNESS", ] diff --git a/gplately/__main__.py b/gplately/__main__.py index 4784afed..25e528d7 100644 --- a/gplately/__main__.py +++ b/gplately/__main__.py @@ -26,11 +26,14 @@ from .commands import ( seafloor_grids, + continent_contouring, feature_filter_cmd, list_models, + paleobathymetry, regrid, reset_feature_type, rotate_grid, + sediment_thickness, ) from .ptt import ( cleanup_topologies, @@ -133,6 +136,15 @@ def main(): # add "rotate_grid" sub-command rotate_grid.add_parser(subparser) + # add "generate-distance-grids"/"generate-sediment-grids" sub-commands + sediment_thickness.add_parser(subparser) + + # add "paleobathymetry" sub-command + paleobathymetry.add_parser(subparser) + + # add "generate-passive-margins" sub-command + continent_contouring.add_parser(subparser) + # add "fix crossovers" sub-command fix_crossovers_cmd = subparser.add_parser( "fix-crossovers", diff --git a/gplately/commands/continent_contouring.py b/gplately/commands/continent_contouring.py new file mode 100644 index 00000000..1a7bf37e --- /dev/null +++ b/gplately/commands/continent_contouring.py @@ -0,0 +1,245 @@ +# +# Copyright (C) 2026 The University of Sydney, Australia +# +# This program is free software; you can redistribute it and/or modify it under +# the terms of the GNU General Public License, version 2, as published by +# the Free Software Foundation. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +# + +import argparse +import logging + +from plate_model_manager import PlateModelManager + +from ..grids.continent_contouring import generate_passive_margins +from .sediment_thickness import _time_range + +_logger = logging.getLogger("gplately") + + +def _run_generate_passive_margins(args): + times = _time_range(args.min_time, args.max_time, args.time_step) + + plate_model = None + if args.model_name: + plate_model = PlateModelManager().get_model( + args.model_name, data_dir=args.plate_model_repo + ) + if not plate_model: + raise Exception( + f"Unable to create PlateModel object for model {args.model_name}." + ) + + rotation_files = args.rotation_filenames or ( + plate_model.get_rotation_model() if plate_model else None + ) + topology_files = args.topology_filenames or ( + plate_model.get_layer("Topologies", return_none_if_not_exist=True) + if plate_model + else None + ) + continent_files = args.continent_filenames or ( + plate_model.get_layer("ContinentalPolygons", return_none_if_not_exist=True) + if plate_model + else None + ) + if continent_files and plate_model and "Cratons" in plate_model.get_avail_layers(): + # Cratons are stored as a separate layer for some models; merge them in, matching + # gplately agegrid's own continent-file resolution (commands/seafloor_grids.py). + continent_files = continent_files + plate_model.get_layer("Cratons") + if not rotation_files or not topology_files: + raise Exception( + "No rotation/topology files found: use -m/--model, or --rotations/--topologies." + ) + if not continent_files: + raise Exception( + "No continental polygon files found: use -m/--model, or --continents." + ) + + _logger.info(f"Using rotation files: {rotation_files}") + _logger.info(f"Using topology files: {topology_files}") + _logger.info(f"Using continent files: {continent_files}") + + generate_passive_margins( + rotation_model=rotation_files, + continent_features=continent_files, + topological_features=topology_files, + times=times, + point_spacing_degrees=args.point_spacing, + area_threshold_square_kms=args.area_threshold_km2, + buffer_and_gap_distance_kms=args.buffer_and_gap_km, + exclusion_area_threshold_square_kms=args.exclusion_area_threshold_km2, + max_distance_of_subduction_from_active_margin_kms=args.max_distance_active_margin_km, + anchor_plate_id=args.anchor_plate_id or 0, + time_step=args.time_step, + output_directory=args.output_dir, + decimal_places_in_time=args.decimal_places_in_time, + ) + _logger.info(f"Passive margins written to {args.output_dir}") + + +def add_parser(parser): + """add command line argument parser for 'generate-passive-margins'""" + + cmd = parser.add_parser( + "generate-passive-margins", + aliases=("gpm",), + help="Dynamically contour continents through time and split each contour into " + "passive margins.", + add_help=True, + description=( + "For each time, contour continental polygons into continents (gplately's " + "ContinentContouring engine), then split each contour into passive-margin " + "segments by removing the parts close to a subduction zone. A port of " + "EarthByte's continent-contouring create_passive_margins.py.\n\n" + "Example usage:\n" + " gplately gpm output_dir -m muller2025 -e 0 -s 100\n" + ), + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + cmd.set_defaults(func=_run_generate_passive_margins) + cmd.add_argument( + metavar="output_dir", + help="(required) output directory", + dest="output_dir", + ) + cmd.add_argument( + "-m", + "--model", + metavar="model_name", + dest="model_name", + default=None, + help="reconstruction model name (fetched via the Plate Model Manager); " + "supplies rotations/topologies/continental polygons unless overridden below", + ) + cmd.add_argument( + "-f", + "--plate-model-repo", + metavar="plate_model_repo", + dest="plate_model_repo", + default="plate-model-repo", + help="local cache directory for -m/--model", + ) + cmd.add_argument( + "--rotations", + metavar="rotation_filenames", + nargs="+", + dest="rotation_filenames", + default=[], + help="alternative to -m/--model", + ) + cmd.add_argument( + "--topologies", + metavar="topology_filenames", + nargs="+", + dest="topology_filenames", + default=[], + help="alternative to -m/--model", + ) + cmd.add_argument( + "--continents", + metavar="continent_filenames", + nargs="+", + dest="continent_filenames", + default=[], + help="continental polygon (or craton) file(s) to contour; alternative to -m/--model", + ) + cmd.add_argument( + "-e", + "--min-time", + metavar="min_time", + type=float, + default=0, + dest="min_time", + help="minimum time (Ma); default: 0", + ) + cmd.add_argument( + "-s", + "--max-time", + metavar="max_time", + type=float, + default=100, + dest="max_time", + help="maximum time (Ma); default: 100", + ) + cmd.add_argument( + "--time-step", + metavar="time_step", + type=float, + default=1, + dest="time_step", + help="time increment (Myr); default: 1", + ) + cmd.add_argument( + "--decimal-places-in-time", + metavar="decimal_places_in_time", + type=int, + default=None, + dest="decimal_places_in_time", + help="decimal places of the reconstruction time in output filenames; default: 0, " + "which reproduces the original workflow's names. Raise it when using a fractional " + "time step, otherwise consecutive times share a filename and only the last is kept", + ) + cmd.add_argument( + "-r", + "--point-spacing", + metavar="point_spacing_degrees", + type=float, + default=0.25, + dest="point_spacing", + help="grid spacing (degrees) used to contour/aggregate continental polygons; " + "default: 0.25", + ) + cmd.add_argument( + "--area-threshold-km2", + metavar="area_threshold_km2", + type=float, + default=0.0, + dest="area_threshold_km2", + help="exclude contoured continents smaller than this (km^2); default: 0", + ) + cmd.add_argument( + "--buffer-and-gap-km", + metavar="buffer_and_gap_km", + type=float, + default=0.0, + dest="buffer_and_gap_km", + help="expand continents ocean-ward by this distance (km) before contouring; " + "default: 0", + ) + cmd.add_argument( + "--exclusion-area-threshold-km2", + metavar="exclusion_area_threshold_km2", + type=float, + default=800000.0, + dest="exclusion_area_threshold_km2", + help="drop enclosed interior gaps (e.g. lakes) smaller than this (km^2); " + "default: 800000", + ) + cmd.add_argument( + "--max-distance-active-margin-km", + metavar="max_distance_active_margin_km", + type=float, + default=500.0, + dest="max_distance_active_margin_km", + help="a contour segment within this distance (km) of a subduction zone is an " + "active margin; default: 500", + ) + cmd.add_argument( + "-a", + "--anchor-plate-id", + metavar="anchor_plate_id", + type=int, + default=None, + dest="anchor_plate_id", + help="anchor plate ID; default: 0", + ) diff --git a/gplately/commands/paleobathymetry.py b/gplately/commands/paleobathymetry.py new file mode 100644 index 00000000..f02a4daf --- /dev/null +++ b/gplately/commands/paleobathymetry.py @@ -0,0 +1,189 @@ +# +# Copyright (C) 2026 The University of Sydney, Australia +# +# This program is free software; you can redistribute it and/or modify it under +# the terms of the GNU General Public License, version 2, as published by +# the Free Software Foundation. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +# + +import argparse +import logging +import os +import tempfile + +import pygplates + +from ..grids.paleobathymetry import AGE_DEPTH_MODELS, simple_paleobathymetry +from .sediment_thickness import ( + _add_common_arguments, + _add_distance_arguments, + _resolve_age_grid_filenames_and_times, + _resolve_distance_grid_kwargs, + _resolve_rotation_topology_proximity_files, +) + +_logger = logging.getLogger("gplately") + + +def _run_paleobathymetry(args): + # The merged static-polygons file below, if one is needed, is written inside this and + # goes away with it. It used to be a NamedTemporaryFile(delete=False), which nothing + # ever deleted -- one abandoned .gpmlz in the system temp directory per run. + # ignore_cleanup_errors: the grids are already written by the time this unwinds, so a + # file still held open here (Windows, in particular) must not turn a finished run into a + # traceback. + with tempfile.TemporaryDirectory( + prefix="gplately-paleobathymetry-", ignore_cleanup_errors=True + ) as scratch_dir: + _run_paleobathymetry_in(args, scratch_dir) + + +def _run_paleobathymetry_in(args, scratch_dir): + age_grid_filenames_and_times, plate_model = _resolve_age_grid_filenames_and_times( + args + ) + rotation_files, topology_files, proximity_files = ( + _resolve_rotation_topology_proximity_files(args, plate_model) + ) + + kwargs = _resolve_distance_grid_kwargs(args, plate_model) + if args.pybacktrack: + static_polygon_filename = args.static_polygons or ( + plate_model.get_layer("StaticPolygons") if plate_model else None + ) + present_day_age_grid_filename = args.present_day_age_grid or ( + plate_model.get_age_grid(0) if plate_model else None + ) + if not static_polygon_filename: + raise Exception( + "--pybacktrack requires --static-polygons (unless -m/--model's plate " + "model provides a StaticPolygons layer)." + ) + if not present_day_age_grid_filename: + raise Exception( + "--pybacktrack requires --present-day-age-grid (unless -m/--model is used, " + "which fetches the age grid at 0 Ma automatically)." + ) + if isinstance(static_polygon_filename, (list, tuple)): + if len(static_polygon_filename) > 1: + # pyBacktrack's static_polygon_filename takes exactly one file; merge + # rather than silently dropping every file but the first. + merged = pygplates.FeatureCollection() + for filename in static_polygon_filename: + merged.add(pygplates.FeatureCollection(filename)) + merged_filename = os.path.join( + scratch_dir, "merged_static_polygons.gpmlz" + ) + merged.write(merged_filename) + _logger.info( + f"Merged {len(static_polygon_filename)} StaticPolygons files into " + f"{merged_filename} for --pybacktrack" + ) + static_polygon_filename = merged_filename + else: + static_polygon_filename = static_polygon_filename[0] + kwargs.update( + pybacktrack=True, + static_polygon_filename=static_polygon_filename, + present_day_age_grid_filename=present_day_age_grid_filename, + ) + + simple_paleobathymetry( + rotation_model=rotation_files, + proximity_features=proximity_files, + topological_features=topology_files, + age_grid_filenames_and_times=age_grid_filenames_and_times, + age_depth_model=args.age_depth_model, + grid_spacing=args.grid_spacing, + time_increment=args.time_increment, + max_reconstruction_time=args.max_reconstruction_time, + anchor_plate_id=args.anchor_plate_id or 0, + clamp_distance_km=args.clamp_distance_km, + richards_table_filename=args.richards_table, + output_directory=args.output_dir, + decimal_places_in_time=args.decimal_places_in_time, + **kwargs, + ) + _logger.info(f"Paleobathymetry grids written to {args.output_dir}") + + +def add_parser(parser): + """add command line argument parser for 'paleobathymetry'""" + + cmd = parser.add_parser( + "paleobathymetry", + aliases=("pb",), + help="Run the simple_paleobathymetry workflow (Steps 1-4, optionally 5) end to end.", + add_help=True, + description=( + "Reconstruct paleobathymetry of ocean crust from a seafloor-age grid: age -> " + "basement depth (Step 1), distance to the nearest passive continental margin " + "(Step 2), predicted sediment thickness (Step 3), and isostatically-compensated " + "paleobathymetry (Step 4). Optionally (--pybacktrack) also merge in pyBacktrack's " + "present-day paleobathymetry (Step 5) to also cover submerged continental crust " + "and crust that has since subducted. A port of EarthByte's simple_paleobathymetry " + "workflow; see gplately.grids.paleobathymetry for the Python API. Step 2 " + "routes distances around continents by default, as the original workflow " + "does; pass --no-route-around-continents for straight-line great-circle " + "distances.\n\n" + "Example usage:\n" + " gplately pb output_dir -m muller2025 --proximity-features cobs.gpml -e 0 -s 10\n" + ), + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + _add_common_arguments(cmd) + _add_distance_arguments(cmd) + cmd.add_argument( + "--age-depth-model", + metavar="age_depth_model", + choices=AGE_DEPTH_MODELS, + default="gdh1", + dest="age_depth_model", + help="thermal-subsidence (age -> basement depth) model; default: gdh1. Also used " + "by --pybacktrack, so Steps 1-4 and Step 5 agree", + ) + cmd.add_argument( + "--richards-table", + metavar="richards_table", + default=None, + dest="richards_table", + help="age-depth lookup table for --age-depth-model rhcw18, used by --pybacktrack " + "as well; " + "default: the table shipped with gplately", + ) + cmd.add_argument( + "--pybacktrack", + action="store_true", + dest="pybacktrack", + help="also run Step 5: merge in pyBacktrack's present-day paleobathymetry, to also " + "cover submerged continental crust and crust that has since subducted. Requires the " + "optional 'pybacktrack' package (pip install pybacktrack, or " + "gplately[paleobathymetry]).", + ) + cmd.add_argument( + "--static-polygons", + metavar="static_polygon_filename", + default=None, + dest="static_polygons", + help="static polygons, for --pybacktrack (pyBacktrack uses these to assign plate IDs); " + "required unless -m/--model's plate model provides a StaticPolygons layer", + ) + cmd.add_argument( + "--present-day-age-grid", + metavar="present_day_age_grid_filename", + default=None, + dest="present_day_age_grid", + help="the seafloor-age grid at 0 Ma, for --pybacktrack (regardless of -e/-s, since " + "pyBacktrack backtracks from the present day); default: fetched automatically " + "when -m/--model is used", + ) + cmd.set_defaults(func=_run_paleobathymetry) diff --git a/gplately/commands/sediment_thickness.py b/gplately/commands/sediment_thickness.py new file mode 100644 index 00000000..cd2955ed --- /dev/null +++ b/gplately/commands/sediment_thickness.py @@ -0,0 +1,484 @@ +# +# Copyright (C) 2026 The University of Sydney, Australia +# +# This program is free software; you can redistribute it and/or modify it under +# the terms of the GNU General Public License, version 2, as published by +# the Free Software Foundation. +# +# This program is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +# for more details. +# +# You should have received a copy of the GNU General Public License along +# with this program; if not, write to Free Software Foundation, Inc., +# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. +# + +import argparse +import logging +import os + +from plate_model_manager import PlateModelManager + +from ..grids import read_netcdf_grid +from ..grids._utils import ( + DEFAULT_DECIMAL_PLACES_IN_TIME, + DEFAULT_DISTANCE_GRID_DECIMAL_PLACES_IN_TIME, + check_times_are_distinct_in_filenames, + distance_grid_filename, + resolve_decimal_places_in_time, +) +from ..grids.sediment_thickness import ( + generate_distance_grids, + generate_sediment_thickness_grids, +) + +_logger = logging.getLogger("gplately") + + +def _time_range(min_time, max_time, time_step): + """min_time, max_time, time_step -> a list of times, min_time to max_time inclusive. + + Unlike ``range(int(min_time), int(max_time) + 1, int(time_step))``, this works for + fractional (sub-Myr) values -- --min-time/--max-time/--time-step are all declared as + floats on the CLI, so truncating them to int here would silently corrupt (or, for a + time_step below 1, crash outright: ``range()`` rejects a step of 0) a perfectly valid + request like ``--time-step 0.5``. + """ + if time_step <= 0: + raise ValueError("time_step must be positive.") + num_steps = int(round((max_time - min_time) / time_step)) + 1 + return [min_time + i * time_step for i in range(num_steps)] + + +def _resolve_age_grid_filenames_and_times(args): + times = _time_range(args.min_time, args.max_time, args.time_step) + if args.model_name: + plate_model = PlateModelManager().get_model( + args.model_name, data_dir=args.plate_model_repo + ) + if not plate_model: + raise Exception( + f"Unable to create PlateModel object for model {args.model_name}." + ) + age_grid_filenames_and_times = [ + (plate_model.get_age_grid(t), float(t)) for t in times + ] + elif args.age_grid_template: + age_grid_filenames_and_times = [ + (args.age_grid_template.format(time=t), float(t)) for t in times + ] + else: + raise Exception( + "No age grid source given: use -m/--model or --age-grid-template." + ) + return age_grid_filenames_and_times, plate_model if args.model_name else None + + +def _resolve_rotation_topology_proximity_files(args, plate_model): + rotation_files = args.rotation_filenames or ( + plate_model.get_rotation_model() if plate_model else None + ) + topology_files = args.topology_filenames or ( + plate_model.get_layer("Topologies", return_none_if_not_exist=True) + if plate_model + else None + ) + # Deliberately not falling back to plate_model.get_layer("COBs"), as the original + # workflow refuses to in three separate places. A COBs layer traces the whole + # continent-ocean boundary, active margins included, where this workflow wants passive + # margins only -- so ocean points beside a subduction zone come out close to a "margin". + # What the layer contains also varies by model: muller2019's is line segments, while + # merdith2021's and the COB Terranes sets are largely polygons, which are worse again. + proximity_files = args.proximity_filenames or None + if not rotation_files or not topology_files: + raise Exception( + "No rotation/topology files found: use -m/--model, or --rotations/--topologies." + ) + if not proximity_files: + raise Exception( + "No proximity feature files found: use --proximity-features to supply " + "passive-margin continent-ocean-boundary line segments. These are deliberately " + "not taken from -m/--model's plate model: a COBs layer traces the entire " + "continent-ocean boundary, active margins included, so using it reports ocean " + "points beside a subduction zone as close to a passive margin. Use " + "'gplately generate-passive-margins' to produce passive margins, or pass the " + "model's COBs explicitly if that approximation is what you want." + ) + + _logger.info(f"Using rotation files: {rotation_files}") + _logger.info(f"Using topology files: {topology_files}") + _logger.info(f"Using proximity feature files: {proximity_files}") + return rotation_files, topology_files, proximity_files + + +def _resolve_distance_grid_kwargs(args, plate_model): + """Resolve the continent-obstacle-routing kwargs shared by 'generate-distance-grids' and + 'paleobathymetry'. Returns a dict suitable for **-splatting into generate_distance_grids() + (or gplately.grids.paleobathymetry.simple_paleobathymetry()) -- empty if obstacle routing + was not requested. + """ + if args.continent_obstacle_filenames and args.route_around_continents is False: + raise Exception( + "--continent-obstacles and --no-route-around-continents contradict each other: " + "the files would be ignored. Drop one." + ) + if args.route_around_continents is False: + return {} + + continent_obstacle_files = args.continent_obstacle_filenames + if not continent_obstacle_files and plate_model: + try: + continent_obstacle_files = plate_model.get_layer( + "Coastlines", return_none_if_not_exist=True + ) + except Exception as exc: + # This is reached on every -m/--model run now that routing is on by default, and + # obtaining the layer can go to the network. A failure here should not take down + # a run that would otherwise have worked -- unless routing was asked for. + if args.route_around_continents is True: + raise + _logger.warning( + "Could not obtain the plate model's Coastlines layer (%s). Falling back to " + "straight-line great-circle distances; pass --continent-obstacles to route " + "around continents, or --no-route-around-continents to silence this.", + exc, + ) + return {} + + if not continent_obstacle_files: + # Only reachable without --continent-obstacles: files given are always used. + if args.route_around_continents is True: + raise Exception( + "--route-around-continents requires continent/coastline files: use " + "--continent-obstacles, or -m/--model's plate model must provide a " + "Coastlines layer." + ) + # Routing is on by default, but a plate model without coastlines cannot do it. Fall + # back rather than refusing to run at all, and say so -- the distances will differ. + _logger.warning( + "Routing around continents is on by default, but no continent/coastline files " + "are available (no --continent-obstacles, and -m/--model has no Coastlines " + "layer). Falling back to straight-line great-circle distances, which will differ " + "from the original workflow's output. Pass --no-route-around-continents to " + "silence this." + ) + return {} + _logger.info(f"Using continent obstacle files: {continent_obstacle_files}") + + return dict( + continent_obstacle_features=continent_obstacle_files, + shortest_path_grid_subdivision_depth=args.shortest_path_grid_depth, + ) + + +def _run_generate_distance_grids(args): + age_grid_filenames_and_times, plate_model = _resolve_age_grid_filenames_and_times( + args + ) + rotation_files, topology_files, proximity_files = ( + _resolve_rotation_topology_proximity_files(args, plate_model) + ) + distance_grid_kwargs = _resolve_distance_grid_kwargs(args, plate_model) + + generate_distance_grids( + rotation_model=rotation_files, + proximity_features=proximity_files, + topological_features=topology_files, + age_grid_filenames_and_times=age_grid_filenames_and_times, + grid_spacing=args.grid_spacing, + time_increment=args.time_increment, + max_reconstruction_time=args.max_reconstruction_time, + anchor_plate_id=args.anchor_plate_id or 0, + clamp_distance_km=args.clamp_distance_km, + output_directory=args.output_dir, + decimal_places_in_time=args.decimal_places_in_time, + **distance_grid_kwargs, + ) + _logger.info(f"Distance grids written to {args.output_dir}") + + +def _run_generate_sediment_grids(args): + age_grid_filenames_and_times, _ = _resolve_age_grid_filenames_and_times(args) + + # Check the output names before reading anything: the library raises on a collision, + # but only after this function has already read every distance grid off disk. + check_times_are_distinct_in_filenames( + [time for _, time in age_grid_filenames_and_times], + resolve_decimal_places_in_time( + args.decimal_places_in_time, DEFAULT_DECIMAL_PLACES_IN_TIME + ), + "sediment_thickness_{}Ma.nc", + ) + + # Must match the names generate_distance_grids() wrote, which default to one decimal + # place of time rather than the zero used by the paleobathymetry grids. + distance_grid_decimal_places = resolve_decimal_places_in_time( + args.decimal_places_in_time, DEFAULT_DISTANCE_GRID_DECIMAL_PLACES_IN_TIME + ) + distance_grids = {} + missing = [] + for _, t in age_grid_filenames_and_times: + distance_path = os.path.join( + args.distance_grids_dir, + distance_grid_filename(args.grid_spacing, t, distance_grid_decimal_places), + ) + if not os.path.isfile(distance_path): + # 'generate-distance-grids' writes nothing for an age grid it found no usable + # data in, so one absent file is a skip rather than a mistake. All of them + # absent is a mistake. + missing.append(distance_path) + continue + grid, lon, lat = read_netcdf_grid(distance_path, return_grids=True) + distance_grids[t] = (lon, lat, grid) + + if not distance_grids: + raise Exception( + f"No distance grids found in {args.distance_grids_dir!r} for times " + f"{[t for _, t in age_grid_filenames_and_times]}. Check --distance-grids-dir, " + "--grid-spacing and --decimal-places-in-time match the " + "'generate-distance-grids' run that produced them." + ) + if missing: + _logger.warning( + "Skipping %d time(s) with no distance grid: %s", len(missing), missing + ) + age_grid_filenames_and_times = [ + (filename, t) + for filename, t in age_grid_filenames_and_times + if t in distance_grids + ] + + generate_sediment_thickness_grids( + age_grid_filenames_and_times, + distance_grids, + output_directory=args.output_dir, + decimal_places_in_time=args.decimal_places_in_time, + ) + _logger.info(f"Sediment thickness grids written to {args.output_dir}") + + +def _add_common_arguments(cmd): + cmd.add_argument( + metavar="output_dir", + help="(required) output directory", + dest="output_dir", + ) + cmd.add_argument( + "-m", + "--model", + metavar="model_name", + dest="model_name", + default=None, + help="reconstruction model name (fetched via the Plate Model Manager); " + "supplies rotations/topologies/age-grids unless overridden below", + ) + cmd.add_argument( + "-f", + "--plate-model-repo", + metavar="plate_model_repo", + dest="plate_model_repo", + default="plate-model-repo", + help="local cache directory for -m/--model", + ) + cmd.add_argument( + "--age-grid-template", + metavar="age_grid_template", + dest="age_grid_template", + default=None, + help="alternative to -m/--model: a filename template for local age grids, " + "using '{time}' for the reconstruction age in Ma, e.g. 'agegrids/age_{time:.0f}Ma.nc'", + ) + cmd.add_argument( + "-e", + "--min-time", + metavar="min_time", + type=float, + default=0, + dest="min_time", + help="minimum time (Ma); default: 0", + ) + cmd.add_argument( + "-s", + "--max-time", + metavar="max_time", + type=float, + default=100, + dest="max_time", + help="maximum time (Ma); default: 100", + ) + cmd.add_argument( + "--time-step", + metavar="time_step", + type=float, + default=1, + dest="time_step", + help="spacing (Myr) of the times to generate output for, between --min-time and " + "--max-time; default: 1. This selects which times get output; how finely each " + "point's lifetime is sampled is --time-increment, which every output time must be " + "a multiple of (so a fractional step below 1 Myr needs --time-increment set to " + "match)", + ) + cmd.add_argument( + "--decimal-places-in-time", + metavar="decimal_places_in_time", + type=int, + default=None, + dest="decimal_places_in_time", + help="decimal places of the reconstruction time in output filenames; by default 1 " + "for the mean-distance grids and 0 for all the others, which reproduces the " + "original workflows' names. Raise it when using a fractional time step, otherwise " + "consecutive times share a filename and only the last is kept", + ) + cmd.add_argument( + "-r", + "--grid-spacing", + metavar="grid_spacing", + type=float, + default=0.5, + dest="grid_spacing", + help="grid spacing (degrees); default: 0.5", + ) + cmd.add_argument( + "-a", + "--anchor-plate-id", + metavar="anchor_plate_id", + type=int, + default=None, + dest="anchor_plate_id", + help="anchor plate ID; default: 0", + ) + + +def _add_distance_arguments(cmd): + cmd.add_argument( + "--time-increment", + metavar="time_increment", + type=float, + default=1, + dest="time_increment", + help="increment (Myr) used to step the backward reconstruction and sample distance " + "along each ocean point's lifetime; default: 1, as in the original workflow. This is " + "independent of --time-step: coarser output does not mean coarser sampling", + ) + cmd.add_argument( + "--proximity-features", + metavar="proximity_filenames", + nargs="+", + dest="proximity_filenames", + default=[], + help="passive-margin continent-ocean-boundary line-segment file(s). Always " + "required: these are deliberately not taken from -m/--model's plate model, whose " + "COBs layer traces the whole continent-ocean boundary including active margins. " + "Use 'gplately generate-passive-margins' to produce passive margins", + ) + cmd.add_argument( + "--rotations", + metavar="rotation_filenames", + nargs="+", + dest="rotation_filenames", + default=[], + help="alternative to -m/--model", + ) + cmd.add_argument( + "--topologies", + metavar="topology_filenames", + nargs="+", + dest="topology_filenames", + default=[], + help="alternative to -m/--model", + ) + cmd.add_argument( + "--max-reconstruction-time", + metavar="max_reconstruction_time", + type=float, + default=None, + dest="max_reconstruction_time", + help="do not reconstruct ocean points older than this (Ma); default: unlimited", + ) + cmd.add_argument( + "--clamp-distance-km", + metavar="clamp_distance_km", + type=float, + default=3000.0, + dest="clamp_distance_km", + help="clamp mean distances above this (km); default: 3000", + ) + cmd.add_argument( + "--route-around-continents", + action=argparse.BooleanOptionalAction, + # None means "not specified": on, but willing to fall back if the plate model has no + # coastlines. An explicit --route-around-continents is not willing to fall back. + default=None, + dest="route_around_continents", + help="route distances around continents instead of a great-circle straight line, as " + "the original workflow does; on by default. Continent files come from " + "--continent-obstacles, or from -m/--model's Coastlines layer. Pass " + "--no-route-around-continents for straight-line distances", + ) + cmd.add_argument( + "--continent-obstacles", + metavar="continent_obstacle_filenames", + nargs="+", + dest="continent_obstacle_filenames", + default=[], + help="continent/coastline file(s) to route around. Routing is on by default, so " + "this chooses where the files come from rather than switching anything on, and it " + "cannot be combined with --no-route-around-continents", + ) + cmd.add_argument( + "--shortest-path-grid-depth", + metavar="shortest_path_grid_subdivision_depth", + type=int, + default=6, + dest="shortest_path_grid_depth", + help="subdivision depth of the grid used for continent-obstacle routing " + "(spacing = 90/2^depth degrees); default: 6", + ) + + +def add_parser(parser): + """add command line argument parsers for 'generate-distance-grids' and 'generate-sediment-grids'""" + + distance_cmd = parser.add_parser( + "generate-distance-grids", + aliases=("gdg",), + help="Generate distance-to-passive-margin grids (for predicting sediment thickness).", + add_help=True, + description=( + "For each ocean point in a seafloor-age grid, reconstruct it backward through " + "time and compute its lifetime-mean distance to the nearest passive-margin " + "continent-ocean-boundary line segment.\n\n" + "Example usage:\n" + " gplately gdg output_dir -m muller2025 --proximity-features cobs.gpml -e 0 -s 10\n" + ), + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + _add_common_arguments(distance_cmd) + _add_distance_arguments(distance_cmd) + distance_cmd.set_defaults(func=_run_generate_distance_grids) + + sediment_cmd = parser.add_parser( + "generate-sediment-grids", + aliases=("gsg",), + help="Predict sediment-thickness grids from seafloor age and distance-to-margin grids.", + add_help=True, + description=( + "Combine a seafloor-age grid with the distance-to-passive-margin grids from " + "'generate-distance-grids' into predicted compacted sediment-thickness grids " + "(Dutkiewicz et al., 2017).\n\n" + "Example usage:\n" + " gplately gsg output_dir -m muller2025 --distance-grids-dir distances/ -e 0 -s 10\n" + ), + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + _add_common_arguments(sediment_cmd) + sediment_cmd.add_argument( + "--distance-grids-dir", + metavar="distance_grids_dir", + required=True, + dest="distance_grids_dir", + help="directory of mean_distance_d_