diff --git a/gplately/__main__.py b/gplately/__main__.py index 11f86a6c..bd6813e2 100644 --- a/gplately/__main__.py +++ b/gplately/__main__.py @@ -2,17 +2,14 @@ import os import sys import warnings -from typing import ( - List, - Optional, - Sequence, - Union, -) +from typing import List, Optional, Sequence, Union import pygplates from gplately import __version__, feature_filter +from .ptt import fix_crossovers, remove_plate_rotations + def combine_feature_collections(input_files: List[str], output_file: str): """Combine multiple feature collections into one.""" @@ -34,53 +31,6 @@ def _run_combine_feature_collections(args): ) -def filter_feature_collection(args): - """Filter the input feature collection according to command line arguments.""" - input_feature_collection = pygplates.FeatureCollection(args.filter_input_file) - - filters = [] - if args.names: - filters.append( - feature_filter.FeatureNameFilter( - args.names, - exact_match=args.exact_match, - case_sensitive=args.case_sensitive, - ) - ) - elif args.exclude_names: - filters.append( - feature_filter.FeatureNameFilter( - args.exclude_names, - exclude=True, - exact_match=args.exact_match, - case_sensitive=args.case_sensitive, - ) - ) - - if args.pids: - filters.append(feature_filter.PlateIDFilter(args.pids)) - elif args.exclude_pids: - filters.append(feature_filter.PlateIDFilter(args.exclude_pids, exclude=True)) - - # print(args.max_birth_age) - if args.max_birth_age is not None: - filters.append( - feature_filter.BirthAgeFilter(args.max_birth_age, keep_older=False) - ) - elif args.min_birth_age is not None: - filters.append(feature_filter.BirthAgeFilter(args.min_birth_age)) - - new_fc = feature_filter.filter_feature_collection( - input_feature_collection, - filters, - ) - - new_fc.write(args.filter_output_file) - print( - f"Done! The filtered feature collection has been saved to {args.filter_output_file}." - ) - - def create_agegrids( input_filenames: Union[str, Sequence[str]], continents_filenames: Union[str, Sequence[str]], @@ -97,15 +47,9 @@ def create_agegrids( unmasked: bool = False, ) -> None: """Create age grids for a plate model.""" - from gplately import ( - PlateReconstruction, - PlotTopologies, - SeafloorGrid, - ) + from gplately import PlateReconstruction, PlotTopologies, SeafloorGrid - features = pygplates.FeaturesFunctionArgument( - input_filenames - ).get_features() + features = pygplates.FeaturesFunctionArgument(input_filenames).get_features() rotations = [] topologies = [] for i in features: @@ -123,9 +67,7 @@ def create_agegrids( continents = pygplates.FeatureCollection() else: continents = pygplates.FeatureCollection( - pygplates.FeaturesFunctionArgument( - continents_filenames - ).get_features() + pygplates.FeaturesFunctionArgument(continents_filenames).get_features() ) with warnings.catch_warnings(): @@ -192,16 +134,17 @@ def main(): title="subcommands", description="valid subcommands", ) + + # add combine feature sub-command combine_cmd = subparser.add_parser( "combine", help=combine_feature_collections.__doc__, description=combine_feature_collections.__doc__, ) - filter_cmd = subparser.add_parser( - "filter", - help=filter_feature_collection.__doc__, - description=filter_feature_collection.__doc__, - ) + + # add feature filter sub-command + feature_filter.add_parser(subparser) + agegrid_cmd = subparser.add_parser( "agegrid", aliases=("ag",), @@ -210,38 +153,28 @@ def main(): description=create_agegrids.__doc__, ) + # add fix crossovers sub-command + fix_crossovers_cmd = subparser.add_parser( + "fix_crossovers", + help="fix crossovers", + add_help=True, + ) + fix_crossovers.add_arguments(fix_crossovers_cmd) + + # add remove plate rotations sub-command + remove_plate_rotations_cmd = subparser.add_parser( + "remove_rotations", + help="remove plate rotations", + add_help=True, + ) + remove_plate_rotations.add_arguments(remove_plate_rotations_cmd) + # combine command arguments combine_cmd.set_defaults(func=_run_combine_feature_collections) combine_cmd.add_argument("combine_first_input_file", type=str) combine_cmd.add_argument("combine_other_input_files", nargs="+", type=str) combine_cmd.add_argument("combine_output_file", type=str) - # feature filter command arguments - filter_cmd.set_defaults(func=filter_feature_collection) - filter_cmd.add_argument("filter_input_file", type=str) - filter_cmd.add_argument("filter_output_file", type=str) - - name_group = filter_cmd.add_mutually_exclusive_group() - name_group.add_argument("-n", "--names", type=str, dest="names", nargs="+") - name_group.add_argument( - "--exclude-names", type=str, dest="exclude_names", nargs="+" - ) - - pid_group = filter_cmd.add_mutually_exclusive_group() - pid_group.add_argument("-p", "--pids", type=int, dest="pids", nargs="+") - pid_group.add_argument("--exclude-pids", type=int, dest="exclude_pids", nargs="+") - - birth_age_group = filter_cmd.add_mutually_exclusive_group() - birth_age_group.add_argument( - "-a", "--min-birth-age", type=float, dest="min_birth_age" - ) - birth_age_group.add_argument("--max-birth-age", type=float, dest="max_birth_age") - - filter_cmd.add_argument( - "--case-sensitive", dest="case_sensitive", action="store_true" - ) - filter_cmd.add_argument("--exact-match", dest="exact_match", action="store_true") - # agegrid command arguments agegrid_cmd.set_defaults(func=_run_create_agegrids) agegrid_cmd.add_argument( diff --git a/gplately/feature_filter.py b/gplately/feature_filter.py index 70ac30b5..300f1ead 100644 --- a/gplately/feature_filter.py +++ b/gplately/feature_filter.py @@ -131,3 +131,83 @@ def filter_feature_collection( if keep_flag: new_feature_collection.add(feature) return new_feature_collection + + +def add_parser(subparser): + """add feature filter command line argument parser""" + filter_cmd = subparser.add_parser( + "filter", + help=filter_feature_collection.__doc__, + description=filter_feature_collection.__doc__, + ) + + # feature filter command arguments + filter_cmd.set_defaults(func=run_filter_feature_collection) + filter_cmd.add_argument("filter_input_file", type=str) + filter_cmd.add_argument("filter_output_file", type=str) + + name_group = filter_cmd.add_mutually_exclusive_group() + name_group.add_argument("-n", "--names", type=str, dest="names", nargs="+") + name_group.add_argument( + "--exclude-names", type=str, dest="exclude_names", nargs="+" + ) + + pid_group = filter_cmd.add_mutually_exclusive_group() + pid_group.add_argument("-p", "--pids", type=int, dest="pids", nargs="+") + pid_group.add_argument("--exclude-pids", type=int, dest="exclude_pids", nargs="+") + + birth_age_group = filter_cmd.add_mutually_exclusive_group() + birth_age_group.add_argument( + "-a", "--min-birth-age", type=float, dest="min_birth_age" + ) + birth_age_group.add_argument("--max-birth-age", type=float, dest="max_birth_age") + + filter_cmd.add_argument( + "--case-sensitive", dest="case_sensitive", action="store_true" + ) + filter_cmd.add_argument("--exact-match", dest="exact_match", action="store_true") + + +def run_filter_feature_collection(args): + """Filter the input feature collection according to command line arguments.""" + input_feature_collection = pygplates.FeatureCollection(args.filter_input_file) + + filters = [] + if args.names: + filters.append( + FeatureNameFilter( + args.names, + exact_match=args.exact_match, + case_sensitive=args.case_sensitive, + ) + ) + elif args.exclude_names: + filters.append( + FeatureNameFilter( + args.exclude_names, + exclude=True, + exact_match=args.exact_match, + case_sensitive=args.case_sensitive, + ) + ) + + if args.pids: + filters.append(PlateIDFilter(args.pids)) + elif args.exclude_pids: + filters.append(PlateIDFilter(args.exclude_pids, exclude=True)) + + # print(args.max_birth_age) + if args.max_birth_age is not None: + filters.append(BirthAgeFilter(args.max_birth_age, keep_older=False)) + elif args.min_birth_age is not None: + filters.append(BirthAgeFilter(args.min_birth_age)) + + new_fc = filter_feature_collection( + input_feature_collection, + filters, + ) + + new_fc.write(args.filter_output_file) + print( + f"Done! The filtered feature collection has been saved to {args.filter_output_file}." + ) diff --git a/gplately/plot.py b/gplately/plot.py index 6b213fe1..dcc9612b 100644 --- a/gplately/plot.py +++ b/gplately/plot.py @@ -2367,7 +2367,7 @@ def plot_subduction_teeth( projection = ax.projection except AttributeError: print( - "The ax.projection does not exist. You must set project to plot Cartopy maps, such as ax = plt.subplot(211, projection=cartopy.crs.PlateCarree())" + "The ax.projection does not exist. You must set projection to plot Cartopy maps, such as ax = plt.subplot(211, projection=cartopy.crs.PlateCarree())" ) projection = None diff --git a/gplately/ptt/fix_crossovers.py b/gplately/ptt/fix_crossovers.py index 14322d8f..b0ff4b34 100644 --- a/gplately/ptt/fix_crossovers.py +++ b/gplately/ptt/fix_crossovers.py @@ -1,4 +1,3 @@ - """ Copyright (C) 2014 The University of Sydney, Australia @@ -17,104 +16,181 @@ """ from __future__ import print_function + import argparse -import math -import sys import os.path +import sys + import pygplates -DEFAULT_OUTPUT_FILENAME_SUFFIX = '_fixed_crossovers' +DEFAULT_OUTPUT_FILENAME_SUFFIX = "_fixed_crossovers" def _fix_crossovers( - rotation_feature_collections, - crossover_threshold_degrees, - crossover_type_function, - ignore_moving_plates, - debug): - + rotation_feature_collections, + crossover_threshold_degrees, + crossover_type_function, + ignore_moving_plates, + debug, +): crossover_filter = None # Ignore a subset of moving plates if requested. if ignore_moving_plates: - crossover_filter = lambda crossover: crossover.moving_plate_id not in ignore_moving_plates - + crossover_filter = ( + lambda crossover: crossover.moving_plate_id not in ignore_moving_plates + ) + # Synchronise crossovers. crossover_results = [] synchronise_crossovers_success = pygplates.synchronise_crossovers( - rotation_feature_collections, - crossover_filter, - crossover_threshold_degrees, - crossover_type_function, - crossover_results) - + rotation_feature_collections, + crossover_filter, + crossover_threshold_degrees, + crossover_type_function, + crossover_results, + ) + # Print debug output if requested. if debug: num_crossovers_error = 0 num_crossovers_ignored = 0 num_crossovers_passed = 0 num_crossovers_corrected = 0 - + for crossover_result in crossover_results: crossover = crossover_result[0] result = crossover_result[1] - + if crossover.type == pygplates.CrossoverType.synch_old_crossover_and_stages: - type_str = 'synch old crossover and stages' + type_str = "synch old crossover and stages" elif crossover.type == pygplates.CrossoverType.synch_old_crossover_only: - type_str = 'synch old crossover only' - elif crossover.type == pygplates.CrossoverType.synch_young_crossover_and_stages: - type_str = 'synch young crossover and stages' + type_str = "synch old crossover only" + elif ( + crossover.type + == pygplates.CrossoverType.synch_young_crossover_and_stages + ): + type_str = "synch young crossover and stages" elif crossover.type == pygplates.CrossoverType.synch_young_crossover_only: - type_str = 'synch young crossover only' + type_str = "synch young crossover only" elif crossover.type == pygplates.CrossoverType.ignore: - type_str = 'ignore' + type_str = "ignore" else: - type_str = 'unknown' - + type_str = "unknown" + if result == pygplates.CrossoverResult.not_synchronised: num_crossovers_passed += 1 - result_str = 'passed' + result_str = "passed" elif result == pygplates.CrossoverResult.synchronised: num_crossovers_corrected += 1 - result_str = 'corrected' + result_str = "corrected" elif result == pygplates.CrossoverResult.ignored: num_crossovers_ignored += 1 - result_str = 'ignored' + result_str = "ignored" else: num_crossovers_error += 1 - result_str = 'error' - - print('Time({0}), moving_pid({1}), young_fixed_pid({2}), old_fixed_pid({3}), type({4}): {5}'.format( + result_str = "error" + + print( + "Time({0}), moving_pid({1}), young_fixed_pid({2}), old_fixed_pid({3}), type({4}): {5}".format( crossover.time, crossover.moving_plate_id, crossover.young_crossover_fixed_plate_id, crossover.old_crossover_fixed_plate_id, type_str, - result_str)) - - print('Results:') - print(' Total number of crossovers = {0}'.format(len(crossover_results))) - print(' Total errors = {0}'.format(num_crossovers_error)) - print(' Total ignored = {0}'.format(num_crossovers_ignored)) - print(' Total passed = {0}'.format(num_crossovers_passed)) - print(' Total corrected = {0}'.format(num_crossovers_corrected)) - + result_str, + ) + ) + + print("Results:") + print(" Total number of crossovers = {0}".format(len(crossover_results))) + print(" Total errors = {0}".format(num_crossovers_error)) + print(" Total ignored = {0}".format(num_crossovers_ignored)) + print(" Total passed = {0}".format(num_crossovers_passed)) + print(" Total corrected = {0}".format(num_crossovers_corrected)) + return synchronise_crossovers_success - -if __name__ == "__main__": - # Check the imported pygplates version. - required_version = pygplates.Version(12) - if not hasattr(pygplates, 'Version') or pygplates.Version.get_imported_version() < required_version: - print('{0}: Error - imported pygplates version {1} but version {2} or greater is required'.format( - os.path.basename(__file__), pygplates.Version.get_imported_version(), required_version), - file=sys.stderr) - sys.exit(1) - - - __description__ = \ - """Loads one or more input rotation files, fixes any crossovers and saves the rotations to \ +def parse_positive_number(value_string): + """parse and return a positive number""" + try: + value = float(value_string) + except ValueError: + raise argparse.ArgumentTypeError("%s is not a number" % value_string) + + if value < 0: + raise argparse.ArgumentTypeError("%g is not a positive number" % value) + + return value + + +def add_arguments(parser: argparse.ArgumentParser): + """add command line argument parser""" + + parser.formatter_class = argparse.RawDescriptionHelpFormatter + parser.description = __description__ + + parser.set_defaults(func=main) + + parser.add_argument( + "-d", "--debug", action="store_true", help="Print debug output." + ) + parser.add_argument( + "-c", + "--crossover_threshold_degrees", + type=parse_positive_number, + help="If specified then crossovers are fixed only if post-crossover rotation latitude, " + "longitude or angle differ from those in pre-crossover rotation by the specified amount " + "(in degrees). This is useful for some PLATES rotation files that are typically accurate " + "to 2 decimal places (or threshold of 0.01).", + ) + + # Can specify only one of '-x' or '-g'. + crossover_type_group = parser.add_mutually_exclusive_group() + crossover_type_group.add_argument( + "-x", + "--default_xo_ys", + action="store_true", + dest="crossover_type_default_xo_ys", + help="If specified, then if a crossover's type is unknown it will default to " + '"synch old crossover and stages", which is equivalent to the "@xo_ys" comment tag.', + ) + crossover_type_group.add_argument( + "-g", + "--default_xo_ig", + action="store_true", + dest="crossover_type_default_xo_ig", + help="If specified, then if a crossover's type is unknown it will default to " + 'ignoring the crossover, which is equivalent to the "@xo_ig" comment tag.', + ) + + parser.add_argument( + "-i", + "--ignore_moving_plates", + type=parse_positive_number, + nargs="+", + metavar="MOVING_PLATE_ID", + help="If specified then is a list of moving plate ids to ignore when fixing crossovers.", + ) + parser.add_argument( + "-s", + "--output_filename_suffix", + type=str, + default="{0}".format(DEFAULT_OUTPUT_FILENAME_SUFFIX), + help="The suffix to append to each input rotation filename to get each output rotation " + "filename - the default suffix is '{0}'".format(DEFAULT_OUTPUT_FILENAME_SUFFIX), + ) + + parser.add_argument( + "input_rotation_filenames", + type=str, + nargs="+", + metavar="input_rotation_filename", + help="One or more input rotation filenames (original files).", + ) + + +__description__ = """Loads one or more input rotation files, fixes any crossovers and saves the rotations to \ output rotation files. The name of each output file is the associated input filename with a suffix appended. For example: @@ -147,89 +223,96 @@ def _fix_crossovers( NOTE: Separate the positional and optional arguments with '--' (workaround for bug in argparse module). For example... - python %(prog)s -d -c 0.01 -i 201 701 -- input_rotations1.rot input_rotations2.rot - """.format(DEFAULT_OUTPUT_FILENAME_SUFFIX) + %(prog)s -d -c 0.01 -i 201 701 -- input_rotations1.rot input_rotations2.rot + """.format( + DEFAULT_OUTPUT_FILENAME_SUFFIX +) + + +def main(args): + # Check the imported pygplates version. + required_version = pygplates.Version(12) + if ( + not hasattr(pygplates, "Version") + or pygplates.Version.get_imported_version() < required_version + ): + print( + "{0}: Error - imported pygplates version {1} but version {2} or greater is required".format( + os.path.basename(__file__), + pygplates.Version.get_imported_version(), + required_version, + ), + file=sys.stderr, + ) + sys.exit(1) - # The command-line parser. - parser = argparse.ArgumentParser(description = __description__, formatter_class=argparse.RawDescriptionHelpFormatter) - - def parse_positive_number(value_string): - try: - value = float(value_string) - except ValueError: - raise argparse.ArgumentTypeError("%s is not a number" % value_string) - - if value < 0: - raise argparse.ArgumentTypeError("%g is not a positive number" % value) - - return value - - parser.add_argument('-d', '--debug', action="store_true", help='Print debug output.') - parser.add_argument('-c', '--crossover_threshold_degrees', type=parse_positive_number, - help='If specified then crossovers are fixed only if post-crossover rotation latitude, ' - 'longitude or angle differ from those in pre-crossover rotation by the specified amount ' - '(in degrees). This is useful for some PLATES rotation files that are typically accurate ' - 'to 2 decimal places (or threshold of 0.01).') - - # Can specify only one of '-x' or '-g'. - crossover_type_group = parser.add_mutually_exclusive_group() - crossover_type_group.add_argument('-x', '--default_xo_ys', action="store_true", - dest='crossover_type_default_xo_ys', - help='If specified, then if a crossover\'s type is unknown it will default to ' - '"synch old crossover and stages", which is equivalent to the "@xo_ys" comment tag.') - crossover_type_group.add_argument('-g', '--default_xo_ig', action="store_true", - dest='crossover_type_default_xo_ig', - help='If specified, then if a crossover\'s type is unknown it will default to ' - 'ignoring the crossover, which is equivalent to the "@xo_ig" comment tag.') - - parser.add_argument('-i', '--ignore_moving_plates', type=parse_positive_number, nargs='+', - metavar='MOVING_PLATE_ID', - help='If specified then is a list of moving plate ids to ignore when fixing crossovers.') - parser.add_argument('-s', '--output_filename_suffix', type=str, - default='{0}'.format(DEFAULT_OUTPUT_FILENAME_SUFFIX), - help="The suffix to append to each input rotation filename to get each output rotation " - "filename - the default suffix is '{0}'".format(DEFAULT_OUTPUT_FILENAME_SUFFIX)) - - parser.add_argument('input_rotation_filenames', type=str, nargs='+', - metavar='input_rotation_filename', - help='One or more input rotation filenames (original files).') - - # Parse command-line options. - args = parser.parse_args() - file_registry = pygplates.FeatureCollectionFileFormatRegistry() - + # Read/parse the input rotation feature collections. - rotation_feature_collections = [file_registry.read(input_rotation_filename) - for input_rotation_filename in args.input_rotation_filenames] - + rotation_feature_collections = [ + file_registry.read(input_rotation_filename) + for input_rotation_filename in args.input_rotation_filenames + ] + # Whether to crossover types should default to '@xo_ys' or '@xo_ig' if type not found... if args.crossover_type_default_xo_ys: - crossover_type_function = pygplates.CrossoverTypeFunction.type_from_xo_tags_in_comment_default_xo_ys + crossover_type_function = ( + pygplates.CrossoverTypeFunction.type_from_xo_tags_in_comment_default_xo_ys + ) elif args.crossover_type_default_xo_ig: - crossover_type_function = pygplates.CrossoverTypeFunction.type_from_xo_tags_in_comment_default_xo_ig + crossover_type_function = ( + pygplates.CrossoverTypeFunction.type_from_xo_tags_in_comment_default_xo_ig + ) else: - crossover_type_function = pygplates.CrossoverTypeFunction.type_from_xo_tags_in_comment - + crossover_type_function = ( + pygplates.CrossoverTypeFunction.type_from_xo_tags_in_comment + ) + # Fix crossovers. # If any errors occurred we will still write to the output files. if not _fix_crossovers( - rotation_feature_collections, - args.crossover_threshold_degrees, - crossover_type_function, - args.ignore_moving_plates, - args.debug): - print('Warning: One or more crossovers were not processed since unable to determine crossover ' - 'correction method or infinite cycle detected.', - file=sys.stderr) - + rotation_feature_collections, + args.crossover_threshold_degrees, + crossover_type_function, + args.ignore_moving_plates, + args.debug, + ): + print( + "Warning: One or more crossovers were not processed since unable to determine crossover " + "correction method or infinite cycle detected.", + file=sys.stderr, + ) + # Write the modified rotation feature collections to disk. for rotation_feature_collection_index in range(len(rotation_feature_collections)): - rotation_feature_collection = rotation_feature_collections[rotation_feature_collection_index] + rotation_feature_collection = rotation_feature_collections[ + rotation_feature_collection_index + ] # Each output filename is the input filename with a suffix appended. - input_rotation_filename = args.input_rotation_filenames[rotation_feature_collection_index] + input_rotation_filename = args.input_rotation_filenames[ + rotation_feature_collection_index + ] filename_root, filename_ext = os.path.splitext(input_rotation_filename) - output_rotation_filename = ''.join((filename_root, args.output_filename_suffix, filename_ext)) - + output_rotation_filename = "".join( + (filename_root, args.output_filename_suffix, filename_ext) + ) + file_registry.write(rotation_feature_collection, output_rotation_filename) + + +if __name__ == "__main__": + # The command-line parser. + parser = argparse.ArgumentParser( + description=__description__, + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + + # add arguments + add_arguments(parser) + + # Parse command-line options. + args = parser.parse_args() + + # call main function + main(args) diff --git a/gplately/ptt/remove_plate_rotations.py b/gplately/ptt/remove_plate_rotations.py index 0ea86c79..9367fef6 100644 --- a/gplately/ptt/remove_plate_rotations.py +++ b/gplately/ptt/remove_plate_rotations.py @@ -1,4 +1,3 @@ - """ Copyright (C) 2019 The University of Sydney, Australia @@ -23,10 +22,14 @@ from __future__ import print_function -import sys + +import argparse import math -import pygplates +import os.path +import sys +import traceback +import pygplates # Required pygplates version. # Need pygplates.RotationModel to clone rotation features by default (revision 12). @@ -34,23 +37,20 @@ PYGPLATES_VERSION_REQUIRED = pygplates.Version(12) -def remove_plates( - rotation_feature_collections, - plate_ids, - accuracy_parameters=None): +def remove_plates(rotation_feature_collections, plate_ids, accuracy_parameters=None): # Docstring in numpydoc format... """Remove one or more plate IDs from a rotation model (consisting of one or more rotation feature collections). - + Any rotations with a fixed plate referencing one of the removed plates will be adjusted such that the rotation model effectively remains unchanged. - + Optional accuracy threshold parameters can be specified to ensure the rotation model after removing plate rotations is very similar to the rotation model before removal. - + The results are returned as a list of pygplates.FeatureCollection (one per input rotation feature collection). - + Ensure you specify all input rotation feature collections that contain the plate IDs to be removed (either as a moving or fixed plate ID). - + Parameters ---------- rotation_feature_collections : sequence of (str, or sequence of pygplates.Feature, or pygplates.FeatureCollection, or pygplates.Feature) @@ -68,23 +68,24 @@ def remove_plates( This mid-way adaptive bisection is repeated (when there is inaccuracy) until the interval between samples becomes smaller than the second parameter (threshold time interval). Rotation threshold is in degrees and threshold time interval is in millions of years. - + Returns ------- list of pygplates.FeatureCollection The (potentially) modified feature collections. Returned list is same length as ``rotation_feature_collections``. """ - + # Convert each feature collection into a list of features so we can more easily remove features # and insert features at arbitrary locations within each feature collection (for example removing # a plate sequence and replacing it with a sequence with the same plate ID). - rotation_feature_collections = [list(pygplates.FeatureCollection(rotation_feature_collection)) - for rotation_feature_collection in rotation_feature_collections] - + rotation_feature_collections = [ + list(pygplates.FeatureCollection(rotation_feature_collection)) + for rotation_feature_collection in rotation_feature_collections + ] + # Iterate over the plates to be removed and remove each one separately. for remove_plate_id in plate_ids: - # Rotation model before any modifications to rotation features. # # Previously we only created one instance of RotationModel in this function (to serve all removed plates). @@ -101,7 +102,7 @@ def remove_plates( # However the RotationModel in the next loop iteration will be affected of course. # UPDATE: Since pygplates revision 25 cloning is no longer necessary (and has been deprecated). rotation_model = pygplates.RotationModel(rotation_feature_collections) - + # Rotation sequences with the current remove plate ID as the *moving* plate ID. # Each sequence will have a different *fixed* plate ID. remove_plate_sequences = [] @@ -109,97 +110,152 @@ def remove_plates( rotation_feature_index = 0 while rotation_feature_index < len(rotation_feature_collection): rotation_feature = rotation_feature_collection[rotation_feature_index] - total_reconstruction_pole = rotation_feature.get_total_reconstruction_pole() + total_reconstruction_pole = ( + rotation_feature.get_total_reconstruction_pole() + ) if total_reconstruction_pole: - fixed_plate_id, moving_plate_id, rotation_sequence = total_reconstruction_pole + ( + fixed_plate_id, + moving_plate_id, + rotation_sequence, + ) = total_reconstruction_pole if moving_plate_id == remove_plate_id: - sample_times = [pygplates.GeoTimeInstant(sample.get_time()) - for sample in rotation_sequence.get_enabled_time_samples()] + sample_times = [ + pygplates.GeoTimeInstant(sample.get_time()) + for sample in rotation_sequence.get_enabled_time_samples() + ] if sample_times: remove_plate_sequences.append( - (fixed_plate_id, sample_times)) + (fixed_plate_id, sample_times) + ) # Remove plate sequences whose moving plate is the current remove plate. # Note that this won't affect 'rotation_model' (since it used a cloned version of all features). del rotation_feature_collection[rotation_feature_index] rotation_feature_index -= 1 - + rotation_feature_index += 1 - + # Sort the remove plate sequences in time order. # This helps out below, to find the max sample time over the removed sequences (all having same *moving* plate ID). # # Easiest way to do this is to sort based on the first time sample in each sequence # (since each sequence should already be sorted internally). - remove_plate_sequences.sort(key = lambda sequence: sequence[1][0]) - + remove_plate_sequences.sort(key=lambda sequence: sequence[1][0]) + # Find those sequences that need adjustment due to the plate removal. # These are sequences whose *fixed* plate is the plate currently being removed. for rotation_feature_collection in rotation_feature_collections: rotation_feature_index = 0 while rotation_feature_index < len(rotation_feature_collection): rotation_feature = rotation_feature_collection[rotation_feature_index] - total_reconstruction_pole = rotation_feature.get_total_reconstruction_pole() + total_reconstruction_pole = ( + rotation_feature.get_total_reconstruction_pole() + ) if total_reconstruction_pole: - fixed_plate_id, moving_plate_id, rotation_sequence = total_reconstruction_pole + ( + fixed_plate_id, + moving_plate_id, + rotation_sequence, + ) = total_reconstruction_pole if fixed_plate_id == remove_plate_id: child_remove_plate_id = moving_plate_id child_remove_plate_rotation_feature = rotation_feature - child_remove_plate_samples = rotation_sequence.get_time_samples() - - child_remove_plate_sample_times = [pygplates.GeoTimeInstant(sample.get_time()) - for sample in child_remove_plate_samples] - child_remove_plate_min_sample_time = child_remove_plate_sample_times[0] - child_remove_plate_max_sample_time = child_remove_plate_sample_times[-1] - + child_remove_plate_samples = ( + rotation_sequence.get_time_samples() + ) + + child_remove_plate_sample_times = [ + pygplates.GeoTimeInstant(sample.get_time()) + for sample in child_remove_plate_samples + ] + child_remove_plate_min_sample_time = ( + child_remove_plate_sample_times[0] + ) + child_remove_plate_max_sample_time = ( + child_remove_plate_sample_times[-1] + ) + # Iterate over the removed sequences whose moving plate matched the current plate being removed. - for remove_plate_sequence_index, (parent_remove_plate_id, remove_plate_sample_times) in enumerate(remove_plate_sequences): + for remove_plate_sequence_index, ( + parent_remove_plate_id, + remove_plate_sample_times, + ) in enumerate(remove_plate_sequences): remove_plate_min_sample_time = remove_plate_sample_times[0] remove_plate_max_sample_time = remove_plate_sample_times[-1] - + # Find the time overlap of the removed sequence and the (child) sequence requiring modification. - min_sample_time = max(remove_plate_min_sample_time, child_remove_plate_min_sample_time) - if remove_plate_sequence_index == len(remove_plate_sequences) - 1: + min_sample_time = max( + remove_plate_min_sample_time, + child_remove_plate_min_sample_time, + ) + if ( + remove_plate_sequence_index + == len(remove_plate_sequences) - 1 + ): # We want the last remove plate sequence to go back to the child max sample time. # If it doesn't go that far back then we will artificially extend the remove plate sequence that far back. max_sample_time = child_remove_plate_max_sample_time else: - max_sample_time = min(remove_plate_max_sample_time, child_remove_plate_max_sample_time) - + max_sample_time = min( + remove_plate_max_sample_time, + child_remove_plate_max_sample_time, + ) + # Note that the remove sequences are ordered by time (ie, first sequence should start at 0Ma, etc). # The two sequences must overlap. # Note that this excludes the case where the min of one sequence equals the max of the other (or max and min). if min_sample_time < max_sample_time: sample_times = [] # Find those sample times of the child sequence within the overlap range. - for child_remove_plate_sample_time in child_remove_plate_sample_times: - if (child_remove_plate_sample_time >= min_sample_time and - child_remove_plate_sample_time <= max_sample_time): - sample_times.append(child_remove_plate_sample_time) + for ( + child_remove_plate_sample_time + ) in child_remove_plate_sample_times: + if ( + child_remove_plate_sample_time + >= min_sample_time + and child_remove_plate_sample_time + <= max_sample_time + ): + sample_times.append( + child_remove_plate_sample_time + ) # Find those sample times of the remove sequence within the overlap range. # Also avoiding duplicating sample times (times already in the child sequence). - for remove_plate_sample_time in remove_plate_sample_times: + for ( + remove_plate_sample_time + ) in remove_plate_sample_times: # Only add the sample time if it's not already in the list. - if (remove_plate_sample_time not in child_remove_plate_sample_times and - remove_plate_sample_time >= min_sample_time and - remove_plate_sample_time <= max_sample_time): + if ( + remove_plate_sample_time + not in child_remove_plate_sample_times + and remove_plate_sample_time >= min_sample_time + and remove_plate_sample_time <= max_sample_time + ): sample_times.append(remove_plate_sample_time) # Need to sort the sample times (since they're likely interleaved between remove and child sequences). sample_times.sort() - + # Gather the rotation samples from the child's moving plate to the removed plate's fixed plate. - parent_to_child_rotation_samples = _merge_rotation_samples( - rotation_model, - child_remove_plate_id, - remove_plate_id, - parent_remove_plate_id, - child_remove_plate_samples, - child_remove_plate_sample_times, - sample_times, - remove_plate_max_sample_time) - + parent_to_child_rotation_samples = ( + _merge_rotation_samples( + rotation_model, + child_remove_plate_id, + remove_plate_id, + parent_remove_plate_id, + child_remove_plate_samples, + child_remove_plate_sample_times, + sample_times, + remove_plate_max_sample_time, + ) + ) + # Insert new samples at times where the difference between original and new rotation models exceeds a threshold. if accuracy_parameters is not None: - threshold_rotation_accuracy_degrees, threshold_time_interval, use_uniform_accuracy_times = accuracy_parameters + ( + threshold_rotation_accuracy_degrees, + threshold_time_interval, + use_uniform_accuracy_times, + ) = accuracy_parameters _ensure_sequence_accuracy( rotation_model, parent_to_child_rotation_samples, @@ -209,32 +265,43 @@ def remove_plates( remove_plate_max_sample_time, threshold_rotation_accuracy_degrees, threshold_time_interval, - use_uniform_accuracy_times) - + use_uniform_accuracy_times, + ) + # Create a new rotation sequence. parent_to_child_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( parent_remove_plate_id, child_remove_plate_id, - pygplates.GpmlIrregularSampling(parent_to_child_rotation_samples), - child_remove_plate_rotation_feature.get_name(None), # Note: specifying None avoids a pygplates crash in revs < 20 - child_remove_plate_rotation_feature.get_description(None)) # Note: specifying None avoids a pygplates crash in revs < 20 - + pygplates.GpmlIrregularSampling( + parent_to_child_rotation_samples + ), + child_remove_plate_rotation_feature.get_name( + None + ), # Note: specifying None avoids a pygplates crash in revs < 20 + child_remove_plate_rotation_feature.get_description( + None + ), + ) # Note: specifying None avoids a pygplates crash in revs < 20 + # Insert the new rotation feature to the current location in the feature collection. # This is better than adding to the end of the collection and thus reordering the order of rotation sequences # in the output collection/file (making it harder to visually find it in a text editor). # Also note that this won't affect 'rotation_model' (since it used a cloned version of all features). - rotation_feature_collection.insert(rotation_feature_index, parent_to_child_rotation_feature) + rotation_feature_collection.insert( + rotation_feature_index, + parent_to_child_rotation_feature, + ) rotation_feature_index += 1 - + # The original rotation feature will no longer be needed because we remove plate sequences # whose fixed plate is the current remove plate. # We would have added one or more sequences above to replace it though. # Also note that this won't affect 'rotation_model' (since it used a cloned version of all features). del rotation_feature_collection[rotation_feature_index] rotation_feature_index -= 1 - + rotation_feature_index += 1 - + # Note that we don't join rotation sequences having the same moving/fixed plates. # However they will show up as a 'duplicate geo-time' warning when loading into GPlates. # TODO: Remove duplicate geo-times and join the offending rotation sequences. @@ -244,23 +311,26 @@ def remove_plates( # by the fixed plate of the removed sequence. In this situation the original crossover sequence # (really two sequences) could now have the same fixed plate ID, and since it also has the # same moving plate ID it should really be one sequence. - + # Return our (potentially) modified feature collections as a list of pygplates.FeatureCollection. - return [pygplates.FeatureCollection(rotation_feature_collection) - for rotation_feature_collection in rotation_feature_collections] + return [ + pygplates.FeatureCollection(rotation_feature_collection) + for rotation_feature_collection in rotation_feature_collections + ] def _merge_rotation_samples( - rotation_model, - child_remove_plate_id, - remove_plate_id, - parent_remove_plate_id, - child_remove_plate_samples, - child_remove_plate_sample_times, - sample_times, - remove_plate_max_sample_time): + rotation_model, + child_remove_plate_id, + remove_plate_id, + parent_remove_plate_id, + child_remove_plate_samples, + child_remove_plate_sample_times, + sample_times, + remove_plate_max_sample_time, +): """Gather the rotation samples from the child's moving plate to the removed plate's fixed plate.""" - + # Gather the rotation samples from the child's moving plate to the removed plate's fixed plate. parent_to_child_rotation_samples = [] for sample_time in sample_times: @@ -273,7 +343,7 @@ def _merge_rotation_samples( # not the 'fixed_plate_id' argument, in case there is no plate circuit path to anchor plate 000. # This also means the user doesn't have to load all rotations in the model, # only those that have the remove plate IDs as a moving or fixed plate. - + if sample_time > remove_plate_max_sample_time: # The time span of the (oldest) removed plate sequence is too short, so extend its oldest sample further into # the past (ie, assume a constant rotation). We do this by calculating R(0->t,parent_plate->remove_plate) @@ -281,75 +351,97 @@ def _merge_rotation_samples( # # Note that the remove sequences are ordered by time (ie, first sequence should start at 0Ma, etc). # So 'max_sample_time' should be the oldest time of the oldest removed plate sequence. - + # R(0->t,parent_plate->remove_plate) parent_to_remove_rotation = rotation_model.get_rotation( # Note the time is 'remove_plate_max_sample_time' and not 'sample_time'... - remove_plate_max_sample_time, remove_plate_id, anchor_plate_id=parent_remove_plate_id) + remove_plate_max_sample_time, + remove_plate_id, + anchor_plate_id=parent_remove_plate_id, + ) # R(0->t,remove_plate->child_plate) remove_to_child_rotation = rotation_model.get_rotation( - sample_time, child_remove_plate_id, anchor_plate_id=remove_plate_id) - - parent_to_child_rotation = parent_to_remove_rotation * remove_to_child_rotation + sample_time, child_remove_plate_id, anchor_plate_id=remove_plate_id + ) + + parent_to_child_rotation = ( + parent_to_remove_rotation * remove_to_child_rotation + ) else: # Note that here we don't actually need to compose rotations as in the above equation because both rotations # are at the same (sample) time so we can just get pygplates.RotationModel.get_rotation() to compose them for us. # # R(0->t,parent_plate->child_plate) parent_to_child_rotation = rotation_model.get_rotation( - sample_time, child_remove_plate_id, anchor_plate_id=parent_remove_plate_id) - + sample_time, + child_remove_plate_id, + anchor_plate_id=parent_remove_plate_id, + ) + # If the sample time corresponds to an existing sample then use its description, # otherwise create a new sample description noting that the new sample is due to # the removal of a specific fixed plate. if sample_time in child_remove_plate_sample_times: - child_remove_plate_sample = child_remove_plate_samples[child_remove_plate_sample_times.index(sample_time)] + child_remove_plate_sample = child_remove_plate_samples[ + child_remove_plate_sample_times.index(sample_time) + ] sample_description = child_remove_plate_sample.get_description() else: - sample_description = 'Removed fixed plate {0}'.format(remove_plate_id) - + sample_description = "Removed fixed plate {0}".format(remove_plate_id) + parent_to_child_rotation_samples.append( pygplates.GpmlTimeSample( pygplates.GpmlFiniteRotation(parent_to_child_rotation), sample_time, - sample_description)) - + sample_description, + ) + ) + return parent_to_child_rotation_samples def _ensure_sequence_accuracy( - rotation_model, - parent_to_child_rotation_samples, - child_remove_plate_id, - remove_plate_id, - parent_remove_plate_id, - remove_plate_max_sample_time, - threshold_rotation_accuracy_degrees, - threshold_time_interval, - insert_poles_at_integer_multiples_of_time_interval): + rotation_model, + parent_to_child_rotation_samples, + child_remove_plate_id, + remove_plate_id, + parent_remove_plate_id, + remove_plate_max_sample_time, + threshold_rotation_accuracy_degrees, + threshold_time_interval, + insert_poles_at_integer_multiples_of_time_interval, +): """Insert new samples at times where the difference between original and new rotation models exceeds a threshold.""" - + num_original_samples = len(parent_to_child_rotation_samples) sample_pair_stack = [] - + # Add the stage rotation intervals to the stack for later processing. - for sample_index in range(num_original_samples-1): - sample1, sample2 = parent_to_child_rotation_samples[sample_index], parent_to_child_rotation_samples[sample_index + 1] + for sample_index in range(num_original_samples - 1): + sample1, sample2 = ( + parent_to_child_rotation_samples[sample_index], + parent_to_child_rotation_samples[sample_index + 1], + ) sample_time1, sample_time2 = sample1.get_time(), sample2.get_time() - if pygplates.GeoTimeInstant(sample_time2 - sample_time1) > threshold_time_interval: + if ( + pygplates.GeoTimeInstant(sample_time2 - sample_time1) + > threshold_time_interval + ): sample_pair_stack.append((sample1, sample2)) - + # Process the stage rotation intervals on the stack until it is empty. while sample_pair_stack: sample1, sample2 = sample_pair_stack.pop() sample_time1, sample_time2 = sample1.get_time(), sample2.get_time() - + mid_sample_time = 0.5 * (sample_time1 + sample_time2) - + if insert_poles_at_integer_multiples_of_time_interval: # Round to the nearest uniformly spaced interval. - interpolated_sample_time = threshold_time_interval * math.floor((mid_sample_time / threshold_time_interval) + 0.5) + interpolated_sample_time = threshold_time_interval * math.floor( + (mid_sample_time / threshold_time_interval) + 0.5 + ) if interpolated_sample_time > mid_sample_time: if interpolated_sample_time >= pygplates.GeoTimeInstant(sample_time2): # We rounded up and the time was greater-or-equal to the end sample time so subtract one time interval. @@ -360,11 +452,11 @@ def _ensure_sequence_accuracy( # We rounded down and the time was less-or-equal to the start sample time so add one time interval. # This is guaranteed to remain within the start/end range since that range should exceed the time interval. interpolated_sample_time += threshold_time_interval - + else: # Just use the sample midway between 'sample1' and 'sample2'. interpolated_sample_time = mid_sample_time - + interpolated_sample = _create_accurate_sample( rotation_model, interpolated_sample_time, @@ -374,47 +466,61 @@ def _ensure_sequence_accuracy( remove_plate_id, parent_remove_plate_id, remove_plate_max_sample_time, - threshold_rotation_accuracy_degrees) - + threshold_rotation_accuracy_degrees, + ) + if interpolated_sample: parent_to_child_rotation_samples.append(interpolated_sample) - + # Recurse if the time interval between start sample and the interpolated sample exceeds threshold interval. - if pygplates.GeoTimeInstant(interpolated_sample_time - sample_time1) > threshold_time_interval: + if ( + pygplates.GeoTimeInstant(interpolated_sample_time - sample_time1) + > threshold_time_interval + ): sample_pair_stack.append((sample1, interpolated_sample)) - + # Recurse if the time interval between the interpolated sample and end sample exceeds threshold interval. - if pygplates.GeoTimeInstant(sample_time2 - interpolated_sample_time) > threshold_time_interval: + if ( + pygplates.GeoTimeInstant(sample_time2 - interpolated_sample_time) + > threshold_time_interval + ): sample_pair_stack.append((interpolated_sample, sample2)) - + # Sort the sample by time (if we added any new samples). if len(parent_to_child_rotation_samples) > num_original_samples: - parent_to_child_rotation_samples.sort(key = lambda sample: sample.get_time()) + parent_to_child_rotation_samples.sort(key=lambda sample: sample.get_time()) def _create_accurate_sample( - rotation_model, - interpolated_sample_time, - sample1, - sample2, - child_remove_plate_id, - remove_plate_id, - parent_remove_plate_id, - remove_plate_max_sample_time, - threshold_rotation_accuracy_degrees): + rotation_model, + interpolated_sample_time, + sample1, + sample2, + child_remove_plate_id, + remove_plate_id, + parent_remove_plate_id, + remove_plate_max_sample_time, + threshold_rotation_accuracy_degrees, +): """Create a new accurate interpolated sample if the difference between original and new rotation models exceeds a threshold, otherwise None.""" - + # Find the *original* rotation from parent plate to child plate (through removed plate). # # R(0->t,parent_plate->remove_plate) parent_to_remove_rotation = rotation_model.get_rotation( # If the time span of the (oldest) removed plate sequence is too short then extend its oldest rotation to the interpolated-sample time... - min(interpolated_sample_time, remove_plate_max_sample_time), remove_plate_id, anchor_plate_id=parent_remove_plate_id) + min(interpolated_sample_time, remove_plate_max_sample_time), + remove_plate_id, + anchor_plate_id=parent_remove_plate_id, + ) # R(0->t,remove_plate->child_plate) remove_to_child_rotation = rotation_model.get_rotation( - interpolated_sample_time, child_remove_plate_id, anchor_plate_id=remove_plate_id) - original_parent_to_child_rotation = parent_to_remove_rotation * remove_to_child_rotation - + interpolated_sample_time, child_remove_plate_id, anchor_plate_id=remove_plate_id + ) + original_parent_to_child_rotation = ( + parent_to_remove_rotation * remove_to_child_rotation + ) + # Find the *new* rotation from parent plate to child plate (through removed plate). # # This interpolates the newly calculated samples (that go directly from parent to child, ie, not via removed plate). @@ -423,64 +529,61 @@ def _create_accurate_sample( sample2.get_value().get_finite_rotation(), sample1.get_time(), sample2.get_time(), - interpolated_sample_time) - + interpolated_sample_time, + ) + # If original and new parent-to-child rotations differ too much then add a new (accurate) sample at the interpolated-sample time. interpolated_sample = None if not pygplates.FiniteRotation.are_equal( original_parent_to_child_rotation, new_parent_to_child_rotation, - threshold_rotation_accuracy_degrees): - + threshold_rotation_accuracy_degrees, + ): interpolated_sample = pygplates.GpmlTimeSample( pygplates.GpmlFiniteRotation(original_parent_to_child_rotation), interpolated_sample_time, - 'Inserted pole to improve accuracy after removing fixed plate {0}'.format(remove_plate_id)) - + "Inserted pole to improve accuracy after removing fixed plate {0}".format( + remove_plate_id + ), + ) + return interpolated_sample -if __name__ == '__main__': - - import os.path - - - # Check the imported pygplates version. - if not hasattr(pygplates, 'Version') or pygplates.Version.get_imported_version() < PYGPLATES_VERSION_REQUIRED: - print('{0}: Error - imported pygplates version {1} but version {2} or greater is required'.format( - os.path.basename(__file__), pygplates.Version.get_imported_version(), PYGPLATES_VERSION_REQUIRED), - file=sys.stderr) - sys.exit(1) - - - import argparse - - # Action to parse a tuple of accuracy parameters. - class ArgParseAccuracyAction(argparse.Action): - def __call__(self, parser, namespace, values, option_string=None): - # Need two numbers (rotation threshold and threshold time interval). - if len(values) != 2: - parser.error('accuracy must be specified as two numbers (rotation threshold and threshold time interval)') - - try: - # Convert strings to float. - threshold_rotation_accuracy_degrees = float(values[0]) - threshold_time_interval = float(values[1]) - except ValueError: - raise argparse.ArgumentTypeError("encountered a rotation threshold and threshold time interval that is not a number") - - if threshold_rotation_accuracy_degrees < 0 or threshold_rotation_accuracy_degrees > 90: - parser.error('rotation threshold must be in the range [0, 90]') - if threshold_time_interval <= 0: - parser.error('threshold time interval must be positive') - - setattr(namespace, self.dest, (threshold_rotation_accuracy_degrees, threshold_time_interval)) - - - def main(): - - __description__ = \ - """Remove one or more plate IDs from a rotation model (consisting of one or more rotation files). +# Action to parse a tuple of accuracy parameters. +class ArgParseAccuracyAction(argparse.Action): + def __call__(self, parser, namespace, values, option_string=None): + # Need two numbers (rotation threshold and threshold time interval). + if len(values) != 2: + parser.error( + "accuracy must be specified as two numbers (rotation threshold and threshold time interval)" + ) + + try: + # Convert strings to float. + threshold_rotation_accuracy_degrees = float(values[0]) + threshold_time_interval = float(values[1]) + except ValueError: + raise argparse.ArgumentTypeError( + "encountered a rotation threshold and threshold time interval that is not a number" + ) + + if ( + threshold_rotation_accuracy_degrees < 0 + or threshold_rotation_accuracy_degrees > 90 + ): + parser.error("rotation threshold must be in the range [0, 90]") + if threshold_time_interval <= 0: + parser.error("threshold time interval must be positive") + + setattr( + namespace, + self.dest, + (threshold_rotation_accuracy_degrees, threshold_time_interval), + ) + + +__description__ = """Remove one or more plate IDs from a rotation model (consisting of one or more rotation files). Any rotations with a fixed plate referencing one of the removed plates will be adjusted such that the rotation model effectively remains unchanged. @@ -495,85 +598,154 @@ def main(): NOTE: Separate the positional and optional arguments with '--' (workaround for bug in argparse module). For example... - python %(prog)s -p 70 4 3 1 -o removed_ref_frames_ -- rotations.rot - """ - - # The command-line parser. - parser = argparse.ArgumentParser(description = __description__, formatter_class=argparse.RawDescriptionHelpFormatter) - - parser.add_argument('-p', '--plates', type=int, nargs='+', required=True, - metavar='remove_plate_ID', - dest='plate_ids', - help='Plate IDs of one or more plates to remove.') - - parser.add_argument( - '-a', '--accuracy', nargs=2, action=ArgParseAccuracyAction, - metavar=('threshold_rotation_accuracy_degrees', 'threshold_time_interval_My'), - help='Optional accuracy parameters. ' - 'If specified then the first parameter is the threshold rotation accuracy (in degrees) and ' - 'the second parameter is the threshold time interval. ' - 'The first parameter is used to compare the latitude, longitude and angle of two rotations before and ' - 'after removing a plate rotation. If any of those three parameters differ by more than the rotation accuracy (in degrees) then ' - 'samples at times mid-way between samples are inserted to ensure before/after accuracy of rotations. ' - 'This mid-way adaptive bisection is repeated (when there is inaccuracy) until the interval between samples ' - 'becomes smaller than the second parameter (threshold time interval). ' - 'Rotation threshold is in degrees and threshold time interval is in millions of years.') - - parser.add_argument('-u', '--use_uniform_accuracy_times', action="store_true", - help='If specified then rotation poles inserted for accuracy (according to "-a" option) will be restricted to times ' - 'that are integer multiples of the threshold time interval (specified in the "-a" option).') - - parser.add_argument('-o', '--output_filename_prefix', type=str, - metavar='output_filename_prefix', - help='Optional output filename prefix. If one is provided then an output rotation file ' - 'is created for each input rotation file by prefixing the input filenames. ' - 'If no filename prefix is provided then the input files are overwritten.') - - parser.add_argument('input_rotation_filenames', type=str, nargs='+', - metavar='input_rotation_filename', - help='One or more rotation files of a rotation model.') - - # Parse command-line options. - args = parser.parse_args() - - # Initialise accuracy parameters (if used). - accuracy_parameters = None - if args.accuracy: - accuracy_parameters = args.accuracy[0], args.accuracy[1], args.use_uniform_accuracy_times - - # Read the input rotation feature collections. - input_rotation_feature_collections = [pygplates.FeatureCollection(input_rotation_filename) - for input_rotation_filename in args.input_rotation_filenames] - - # Remove plate rotations. - output_rotation_feature_collections = remove_plates( - input_rotation_feature_collections, - args.plate_ids, - accuracy_parameters) - - # Write the modified rotation feature collections to disk. - for rotation_feature_collection_index in range(len(output_rotation_feature_collections)): - output_rotation_feature_collection = output_rotation_feature_collections[rotation_feature_collection_index] - - # Each output filename is the input filename with an optional prefix prepended. - input_rotation_filename = args.input_rotation_filenames[rotation_feature_collection_index] - if args.output_filename_prefix: - dir, file_basename = os.path.split(input_rotation_filename) - output_rotation_filename = os.path.join(dir, '{0}{1}'.format(args.output_filename_prefix, file_basename)) - else: - output_rotation_filename = input_rotation_filename - - output_rotation_feature_collection.write(output_rotation_filename) - - sys.exit(0) - - import traceback - + %(prog)s -p 70 4 3 1 -o removed_ref_frames_ -- rotations.rot + """ + + +def main(args): + # Check the imported pygplates version. + if ( + not hasattr(pygplates, "Version") + or pygplates.Version.get_imported_version() < PYGPLATES_VERSION_REQUIRED + ): + print( + "{0}: Error - imported pygplates version {1} but version {2} or greater is required".format( + os.path.basename(__file__), + pygplates.Version.get_imported_version(), + PYGPLATES_VERSION_REQUIRED, + ), + file=sys.stderr, + ) + sys.exit(1) + + # Initialise accuracy parameters (if used). + accuracy_parameters = None + if args.accuracy: + accuracy_parameters = ( + args.accuracy[0], + args.accuracy[1], + args.use_uniform_accuracy_times, + ) + + # Read the input rotation feature collections. + input_rotation_feature_collections = [ + pygplates.FeatureCollection(input_rotation_filename) + for input_rotation_filename in args.input_rotation_filenames + ] + + # Remove plate rotations. + output_rotation_feature_collections = remove_plates( + input_rotation_feature_collections, args.plate_ids, accuracy_parameters + ) + + # Write the modified rotation feature collections to disk. + for rotation_feature_collection_index in range( + len(output_rotation_feature_collections) + ): + output_rotation_feature_collection = output_rotation_feature_collections[ + rotation_feature_collection_index + ] + + # Each output filename is the input filename with an optional prefix prepended. + input_rotation_filename = args.input_rotation_filenames[ + rotation_feature_collection_index + ] + if args.output_filename_prefix: + dir, file_basename = os.path.split(input_rotation_filename) + output_rotation_filename = os.path.join( + dir, "{0}{1}".format(args.output_filename_prefix, file_basename) + ) + else: + output_rotation_filename = input_rotation_filename + + output_rotation_feature_collection.write(output_rotation_filename) + + sys.exit(0) + + +def add_arguments(parser): + parser.formatter_class = argparse.RawDescriptionHelpFormatter + parser.description = __description__ + + parser.set_defaults(func=main) + + parser.add_argument( + "-p", + "--plates", + type=int, + nargs="+", + required=True, + metavar="remove_plate_ID", + dest="plate_ids", + help="Plate IDs of one or more plates to remove.", + ) + + parser.add_argument( + "-a", + "--accuracy", + nargs=2, + action=ArgParseAccuracyAction, + metavar=( + "threshold_rotation_accuracy_degrees", + "threshold_time_interval_My", + ), + help="Optional accuracy parameters. " + "If specified then the first parameter is the threshold rotation accuracy (in degrees) and " + "the second parameter is the threshold time interval. " + "The first parameter is used to compare the latitude, longitude and angle of two rotations before and " + "after removing a plate rotation. If any of those three parameters differ by more than the rotation accuracy (in degrees) then " + "samples at times mid-way between samples are inserted to ensure before/after accuracy of rotations. " + "This mid-way adaptive bisection is repeated (when there is inaccuracy) until the interval between samples " + "becomes smaller than the second parameter (threshold time interval). " + "Rotation threshold is in degrees and threshold time interval is in millions of years.", + ) + + parser.add_argument( + "-u", + "--use_uniform_accuracy_times", + action="store_true", + help='If specified then rotation poles inserted for accuracy (according to "-a" option) will be restricted to times ' + 'that are integer multiples of the threshold time interval (specified in the "-a" option).', + ) + + parser.add_argument( + "-o", + "--output_filename_prefix", + type=str, + metavar="output_filename_prefix", + help="Optional output filename prefix. If one is provided then an output rotation file " + "is created for each input rotation file by prefixing the input filenames. " + "If no filename prefix is provided then the input files are overwritten.", + ) + + parser.add_argument( + "input_rotation_filenames", + type=str, + nargs="+", + metavar="input_rotation_filename", + help="One or more rotation files of a rotation model.", + ) + + +if __name__ == "__main__": + # The command-line parser. + parser = argparse.ArgumentParser( + description=__description__, + formatter_class=argparse.RawDescriptionHelpFormatter, + ) + + # add arguments + add_arguments(parser) + + # Parse command-line options. + args = parser.parse_args() + + # call main function try: - main() + main(args) sys.exit(0) except Exception as exc: - print('ERROR: {0}'.format(exc), file=sys.stderr) + print("ERROR: {0}".format(exc), file=sys.stderr) # Uncomment this to print traceback to location of raised exception. # traceback.print_exc() sys.exit(1)