From 88a789f3df5a8477867f64bd82f8a59e50dc6a90 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Wed, 4 Jun 2025 14:14:49 +1000 Subject: [PATCH 1/7] First cut at adding reference frames (including opt) to un-optimised. This adds an extra file containing the reference frames, including optimised mantle (6), no-net rotation (7) and an approximation to true polar wander (8) which is paleomag minus optimised mantle reference frames. Plate ID 0 gives the original un-optimised reference frame. --- supplementary/Readme.txt | 11 ++ .../add_ref_frames_to_unoptimised.py | 117 ++++++++++++++++++ 2 files changed, 128 insertions(+) create mode 100644 supplementary/Readme.txt create mode 100644 supplementary/add_ref_frames_to_unoptimised.py diff --git a/supplementary/Readme.txt b/supplementary/Readme.txt new file mode 100644 index 0000000..d3e4f5f --- /dev/null +++ b/supplementary/Readme.txt @@ -0,0 +1,11 @@ +This folder contains supplementary scripts that are not part of the main optimisation workflow. + +add_ref_frames_to_unoptimised.py +-------------------------------- + +This script reads the original un-optimised rotation files, and the optimised and +the no-net rotation reference frames, and generates a new set of rotation files +with various reference frames added. +Currently the added reference frames are the original paleomag (PMAG), the optimised mantle, +the no-net rotation, and an approximation of true polar wander (TPW) calculated by replacing +the PMAG reference frame with the optimised mantle reference frame. \ No newline at end of file diff --git a/supplementary/add_ref_frames_to_unoptimised.py b/supplementary/add_ref_frames_to_unoptimised.py new file mode 100644 index 0000000..cb70f8f --- /dev/null +++ b/supplementary/add_ref_frames_to_unoptimised.py @@ -0,0 +1,117 @@ +import glob +import os.path +import pygplates + + +#### Input parameters #### + +# Original (un-optimised) rotation files. +original_rotation_filenames = [ + '1000_0_rotfile.rot', + '1800_1000_rotfile.rot', +] +#original_rotation_filenames = glob.glob('*.rot') + +# This file contains ALL rotations, including the 005-000 optimised reference frame. +optimised_rotation_file = 'optimisation/optimised_rotation_model_20240725.rot' + +# This file contains ALL rotations, including the 005-000 no-net-rotation reference frame. +no_net_rotation_file = 'optimisation/no_net_rotation_model_20240725.rot' + +# Output rotation file to contain the reference frames being added (optimised, NNR, optimised-minus-pmag). +output_ref_frames_file = 'optimisation/reference_frames_20240725.rot' + +# Plate IDs of reference frames to store in output rotation file. +optimised_ref_frame_plate_id = 6 +no_net_rotation_ref_frame_plate_id = 7 +optimised_minus_pmag_ref_frame_plate_id = 8 + +########################## + + +optimised_rotation_feature_collection = pygplates.FeatureCollection(optimised_rotation_file) +no_net_rotation_feature_collection = pygplates.FeatureCollection(no_net_rotation_file) + +# Original (un-optimised) rotation model. +original_rotation_model = pygplates.RotationModel(original_rotation_filenames) + +def find_005_000_rotation_sequence(rotation_feature_collection): + rotation_feature_005_000 = None + + for rotation_feature in rotation_feature_collection: + total_reconstruction_pole = rotation_feature.get_total_reconstruction_pole() + if total_reconstruction_pole: + fixed_plate_id, moving_plate_id, _ = total_reconstruction_pole + + if fixed_plate_id == 0 and moving_plate_id == 5: + if rotation_feature_005_000 is not None: + raise ValueError("Found multiple rotation sequences with moving plate 5 and fixed plate 0") + rotation_feature_005_000 = rotation_feature + + if rotation_feature_005_000 is None: + raise ValueError("Found zero rotation sequences with moving plate 5 and fixed plate 0") + + return rotation_feature_005_000 + +# Find the optimised (005-000) rotation sequence. +optimised_rotation_feature = find_005_000_rotation_sequence(optimised_rotation_feature_collection) +optimised_rotation_feature = optimised_rotation_feature.clone() +# Change its moving plate ID (from 005 to ). +optimised_rotation_feature.set(pygplates.PropertyName.gpml_moving_reference_frame, pygplates.GpmlPlateId(optimised_ref_frame_plate_id)) +# Reverse its rotations since it'll get activated by setting the anchor plate to instead of 000. +for rotation_sample in optimised_rotation_feature.get_total_reconstruction_pole()[2].get_enabled_time_samples(): + finite_rotation = rotation_sample.get_value().get_finite_rotation() + rotation_sample.get_value().set_finite_rotation(finite_rotation.get_inverse()) + rotation_sample.set_description('optAPM @absage') + +# Find the no-net rotation (005-000) rotation sequence. +no_net_rotation_feature = find_005_000_rotation_sequence(no_net_rotation_feature_collection) +no_net_rotation_feature = no_net_rotation_feature.clone() +# Change its moving plate ID (from 005 to ). +no_net_rotation_feature.set(pygplates.PropertyName.gpml_moving_reference_frame, pygplates.GpmlPlateId(no_net_rotation_ref_frame_plate_id)) +# Reverse its rotations since it'll get activated by setting the anchor plate to instead of 000. +for rotation_sample in no_net_rotation_feature.get_total_reconstruction_pole()[2].get_enabled_time_samples(): + finite_rotation = rotation_sample.get_value().get_finite_rotation() + rotation_sample.get_value().set_finite_rotation(finite_rotation.get_inverse()) + rotation_sample.set_description('NNR') + +# Calculate a new optimised-minus-pmag sequence that combines: +# - a reversal of 701-000 (to remove PMAG frame), and +# - a reversal of existing optimised rotation (005-000) +# - since it'll get activated by setting the anchor plate to instead of 000. +rotation_time_samples_optimised_minus_pmag = [] +for optimised_rotation_sample in optimised_rotation_feature.get_total_reconstruction_pole()[2].get_enabled_time_samples(): + time = optimised_rotation_sample.get_time() + optimised_finite_rotation = optimised_rotation_sample.get_value().get_finite_rotation() + + # Get the 701-000 rotation. + finite_rotation_701_000 = original_rotation_model.get_rotation(time, 701) + + # Remove pmag from optimised reference frame: + # + # R(opt_minus_pmag->000) = R(opt->000) * inverse[R(000->701)] + # = inverse[R(000->opt)] * inverse[R(000->701)] + # R(000->opt_minus_pmag)) = inverse[R(opt_minus_pmag->000)] + # = inverse[inverse[R(000->opt)] * inverse[R(000->701)]] + finite_rotation_optimised_minus_pmag = (optimised_finite_rotation.get_inverse() * finite_rotation_701_000.get_inverse()).get_inverse() + + rotation_time_samples_optimised_minus_pmag.append( + pygplates.GpmlTimeSample( + pygplates.GpmlFiniteRotation(finite_rotation_optimised_minus_pmag), + time, + 'optAPM without pmag')) + +# Create a new optimised-minus-pmag rotation sequence. +optimised_minus_pmag_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( + 0, + optimised_minus_pmag_ref_frame_plate_id, + pygplates.GpmlIrregularSampling(rotation_time_samples_optimised_minus_pmag)) + +# Add reference frames to feature collection. +output_ref_frames_feature_collection = pygplates.FeatureCollection() +output_ref_frames_feature_collection.add(optimised_rotation_feature) +output_ref_frames_feature_collection.add(no_net_rotation_feature) +output_ref_frames_feature_collection.add(optimised_minus_pmag_rotation_feature) + +# Write reference frames to output rotation file. +output_ref_frames_feature_collection.write(output_ref_frames_file) From 3de29ae51c42b6a3807bc4f43d4ad6dbb8d60736 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Wed, 4 Jun 2025 14:37:35 +1000 Subject: [PATCH 2/7] Can choose which reference frame is represented by plate ID 000. - Changed plate IDs of OPT, NNR and TPW from 6, 7, 8 to 15, 16, 17 to reduce potential conflict with existing IDs (eg, hotspots). - Added plate ID 014 for PMAG. - Output rotation files match input files except 000 is changed to 014. - We still keep reference frames (14, 15, 16, 17) in extra file. - Can specify default 000 to be equivalent to any reference frame. --- .../add_ref_frames_to_unoptimised.py | 400 ++++++++++++++---- 1 file changed, 319 insertions(+), 81 deletions(-) diff --git a/supplementary/add_ref_frames_to_unoptimised.py b/supplementary/add_ref_frames_to_unoptimised.py index cb70f8f..3aca145 100644 --- a/supplementary/add_ref_frames_to_unoptimised.py +++ b/supplementary/add_ref_frames_to_unoptimised.py @@ -3,115 +3,353 @@ import pygplates +# +# This script reads the original un-optimised rotation files and the optimised and +# the no-net rotation reference frames, and generates a new set of rotation files +# with various reference frames added. +# +# Currently the added reference frames are the original paleomag (PMAG), the optimised mantle, +# the no-net rotation, and an approximation of true polar wander (TPW) calculated by replacing +# the PMAG reference frame with the optimised mantle reference frame. +# + + #### Input parameters #### # Original (un-optimised) rotation files. -original_rotation_filenames = [ +input_rotation_filenames = [ '1000_0_rotfile.rot', '1800_1000_rotfile.rot', ] -#original_rotation_filenames = glob.glob('*.rot') +#input_rotation_filenames = glob.glob('*.rot') + +# This file contains the 005-000 optimised reference frame. +input_optimised_rotation_file = 'optimisation/optimised_rotation_model_20240725.rot' -# This file contains ALL rotations, including the 005-000 optimised reference frame. -optimised_rotation_file = 'optimisation/optimised_rotation_model_20240725.rot' +# This file contains the 005-000 no-net-rotation reference frame. +input_no_net_rotation_file = 'optimisation/no_net_rotation_model_20240725.rot' -# This file contains ALL rotations, including the 005-000 no-net-rotation reference frame. -no_net_rotation_file = 'optimisation/no_net_rotation_model_20240725.rot' +# The suffix to append to the basenames of the input rotation filenames to get the output rotation filenames. +# +# Note: The output rotation files also end up in the 'optimisation/' sub-directory. +output_rotation_filenames_suffix = '_20240725_2' -# Output rotation file to contain the reference frames being added (optimised, NNR, optimised-minus-pmag). -output_ref_frames_file = 'optimisation/reference_frames_20240725.rot' +# Plate IDs of reference frames. +# +# These are PMAG, optimised, no-net rotation and true polar wander (which is optimised with PMAG removed). +pmag_reference_frame_plate_id = 14 +optimised_ref_frame_plate_id = 15 +no_net_rotation_ref_frame_plate_id = 16 +true_polar_wander_ref_frame_plate_id = 17 -# Plate IDs of reference frames to store in output rotation file. -optimised_ref_frame_plate_id = 6 -no_net_rotation_ref_frame_plate_id = 7 -optimised_minus_pmag_ref_frame_plate_id = 8 +# The default reference frame that plate 000 should correspond to. +# +# It should be one of the above reference frame plate IDs. +default_reference_frame_plate_id = optimised_ref_frame_plate_id + +# +# Note that the final rotation hierarchy will look something like... +# +# 015 016 017 +# | | | +# ----------- +# | +# 000 +# | +# 014 +# | +# ----- +# | | +# 101 701 +# | | +# +# ...where if 000 represents: +# 1. PMAG – then 014-000 is zero, and 000-015 (and 014-015) represents what 005-000 previously represented. +# 2. Optimised – then 000-015 is zero, and 014-000 (and 014-015) represents what 005-000 previously represented. +# +# That can be one version of the rotation model. +# And another can be a flattened version that removes the reference frame plate IDs (eg, 014, 015, 016, 017). +# The flattened version can be what's included with GPlately. +# ########################## -optimised_rotation_feature_collection = pygplates.FeatureCollection(optimised_rotation_file) -no_net_rotation_feature_collection = pygplates.FeatureCollection(no_net_rotation_file) +# Description for the comments of the reference frame rotations. +rotation_description_pmag = ' {:03d} is paleomag reference frame'.format(pmag_reference_frame_plate_id) +rotation_description_optimised = ' {:03d} is optimised mantle reference frame @absage'.format(optimised_ref_frame_plate_id) +rotation_description_no_net_rotation = ' {:03d} is no-net rotation reference frame'.format(no_net_rotation_ref_frame_plate_id) +rotation_description_true_polar_wander = ' {:03d} is true polar wander reference frame (approx)'.format(true_polar_wander_ref_frame_plate_id) -# Original (un-optimised) rotation model. -original_rotation_model = pygplates.RotationModel(original_rotation_filenames) +# Each input rotation file will get a corresponding output rotation file. +output_rotation_feature_collections = [pygplates.FeatureCollection(filename) + for filename in input_rotation_filenames] + +# The output rotation filenames associated with the input rotation files. +output_rotation_filenames = [] +for filename in input_rotation_filenames: + # Insert suffix into output rotation filenames. + basename, ext = os.path.splitext(filename) + output_filename = basename + output_rotation_filenames_suffix + ext + # Place in 'optimisation/' sub-directory. + output_filename = os.path.join('optimisation', output_filename) + output_rotation_filenames.append(output_filename) -def find_005_000_rotation_sequence(rotation_feature_collection): - rotation_feature_005_000 = None +# The output rotation filename to store the optimised, NNR and true polar wander reference frames. +output_reference_frames_filename = 'optimisation/reference_frames{}.rot'.format(output_rotation_filenames_suffix) +########################## + + +# Change fixed plate 000 to +# (and remove any existing leftover 005-000 rotation features). +max_time_of_sequences_with_fixed_000 = 0.0 +for rotation_feature_collection in output_rotation_feature_collections: + + # Exclude any leftover 005-000 rotation features. + def is_feature_005_000(feature): + total_reconstruction_pole = feature.get_total_reconstruction_pole() + if total_reconstruction_pole: + fixed_plate_id, moving_plate_id, _ = total_reconstruction_pole + return fixed_plate_id == 0 and moving_plate_id == 5 + return False + rotation_feature_collection.remove(is_feature_005_000) + + # Change fixed plate 000 to . for rotation_feature in rotation_feature_collection: total_reconstruction_pole = rotation_feature.get_total_reconstruction_pole() if total_reconstruction_pole: - fixed_plate_id, moving_plate_id, _ = total_reconstruction_pole + fixed_plate_id, moving_plate_id, rotation_sequence = total_reconstruction_pole + + # Change fixed plate ID 000 to . + # We want everything that references 000 to now reference . + if fixed_plate_id == 0: + # Change the fixed plate ID from 000 to . + rotation_feature.set_total_reconstruction_pole( + pmag_reference_frame_plate_id, moving_plate_id, rotation_sequence) + # The oldest time in this sequence. + max_time = max(time_sample.get_time() for time_sample in rotation_sequence.get_time_samples()) + # The oldest time in all sequences (with fixed plate 000). + max_time_of_sequences_with_fixed_000 = max(max_time_of_sequences_with_fixed_000, max_time) + +# Original (un-optimised) rotation model. +# +# Note: We use the modified input rotation features (with fixed plate 000 changed to 014). +# +# Note: It has no 000 plate ID. Instead it's anchor plate is the PMAG reference frame plate ID. +original_rotation_model_anchored_014 = pygplates.RotationModel( + output_rotation_feature_collections, + default_anchor_plate_id=pmag_reference_frame_plate_id) + + +def find_005_000_rotation_sequence(rotation_filename): + rotation_sequence_005_000 = None + + for rotation_feature in pygplates.FeatureCollection(rotation_filename): + total_reconstruction_pole = rotation_feature.get_total_reconstruction_pole() + if total_reconstruction_pole: + fixed_plate_id, moving_plate_id, rotation_sequence = total_reconstruction_pole if fixed_plate_id == 0 and moving_plate_id == 5: - if rotation_feature_005_000 is not None: + if rotation_sequence_005_000 is not None: raise ValueError("Found multiple rotation sequences with moving plate 5 and fixed plate 0") - rotation_feature_005_000 = rotation_feature + rotation_sequence_005_000 = rotation_sequence - if rotation_feature_005_000 is None: + if rotation_sequence_005_000 is None: raise ValueError("Found zero rotation sequences with moving plate 5 and fixed plate 0") - return rotation_feature_005_000 + return rotation_sequence_005_000 # Find the optimised (005-000) rotation sequence. -optimised_rotation_feature = find_005_000_rotation_sequence(optimised_rotation_feature_collection) -optimised_rotation_feature = optimised_rotation_feature.clone() -# Change its moving plate ID (from 005 to ). -optimised_rotation_feature.set(pygplates.PropertyName.gpml_moving_reference_frame, pygplates.GpmlPlateId(optimised_ref_frame_plate_id)) -# Reverse its rotations since it'll get activated by setting the anchor plate to instead of 000. -for rotation_sample in optimised_rotation_feature.get_total_reconstruction_pole()[2].get_enabled_time_samples(): - finite_rotation = rotation_sample.get_value().get_finite_rotation() - rotation_sample.get_value().set_finite_rotation(finite_rotation.get_inverse()) - rotation_sample.set_description('optAPM @absage') +optimised_rotation_sequence = find_005_000_rotation_sequence(input_optimised_rotation_file) # Find the no-net rotation (005-000) rotation sequence. -no_net_rotation_feature = find_005_000_rotation_sequence(no_net_rotation_feature_collection) -no_net_rotation_feature = no_net_rotation_feature.clone() -# Change its moving plate ID (from 005 to ). -no_net_rotation_feature.set(pygplates.PropertyName.gpml_moving_reference_frame, pygplates.GpmlPlateId(no_net_rotation_ref_frame_plate_id)) -# Reverse its rotations since it'll get activated by setting the anchor plate to instead of 000. -for rotation_sample in no_net_rotation_feature.get_total_reconstruction_pole()[2].get_enabled_time_samples(): - finite_rotation = rotation_sample.get_value().get_finite_rotation() - rotation_sample.get_value().set_finite_rotation(finite_rotation.get_inverse()) - rotation_sample.set_description('NNR') - -# Calculate a new optimised-minus-pmag sequence that combines: -# - a reversal of 701-000 (to remove PMAG frame), and -# - a reversal of existing optimised rotation (005-000) -# - since it'll get activated by setting the anchor plate to instead of 000. -rotation_time_samples_optimised_minus_pmag = [] -for optimised_rotation_sample in optimised_rotation_feature.get_total_reconstruction_pole()[2].get_enabled_time_samples(): - time = optimised_rotation_sample.get_time() - optimised_finite_rotation = optimised_rotation_sample.get_value().get_finite_rotation() - - # Get the 701-000 rotation. - finite_rotation_701_000 = original_rotation_model.get_rotation(time, 701) - - # Remove pmag from optimised reference frame: +no_net_rotation_sequence = find_005_000_rotation_sequence(input_no_net_rotation_file) + + +def create_identity_rotation_sequence(rotation_description, begin_time, end_time=0.0): + rotation_time_samples = [] + identity_rotation = pygplates.FiniteRotation() + # Only need two time samples at start and end times. + for time in (end_time, begin_time): + rotation_time_samples.append( + pygplates.GpmlTimeSample( + pygplates.GpmlFiniteRotation(identity_rotation), + time, + rotation_description)) + return rotation_time_samples + +def clone_rotation_sequence(rotation_sequence, rotation_description): + rotation_time_samples = [] + for rotation_sample in rotation_sequence.get_enabled_time_samples(): + rotation_time_samples.append( + pygplates.GpmlTimeSample( + pygplates.GpmlFiniteRotation(rotation_sample.get_value().get_finite_rotation()), + rotation_sample.get_time(), + rotation_description)) + return rotation_time_samples + +def calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, rotation_description): + rotation_time_samples = [] + for optimised_rotation_sample in optimised_rotation_sequence.get_enabled_time_samples(): + time = optimised_rotation_sample.get_time() + # R(000->005). + optimised_finite_rotation = optimised_rotation_sample.get_value().get_finite_rotation() + # R(014->701). + finite_rotation_701_014 = original_rotation_model_anchored_014.get_rotation(time, 701) + # We want TPW->701 to be equivalent to the optimised reference frame 000->005: + # + # R(TPW->701) = R(000->005) + # R(TPW->014) * R(014->701) = R(000->005) + # R(TPW->014) = R(000->005) * inverse[R(014->701)] + finite_rotation_true_polar_wander = optimised_finite_rotation * finite_rotation_701_014.get_inverse() + rotation_time_samples.append( + pygplates.GpmlTimeSample( + pygplates.GpmlFiniteRotation(finite_rotation_true_polar_wander), + time, + rotation_description)) + return rotation_time_samples + +# Generate the rotation sequence for the PMAG reference frame. +# +# This depends on which reference frame 000 is assigned to. +if default_reference_frame_plate_id == pmag_reference_frame_plate_id: + # Both 000 and 014 are PMAG. + # + # So 000->014 is zero (identity). + rotation_time_samples_pmag = create_identity_rotation_sequence(rotation_description_pmag, max_time_of_sequences_with_fixed_000) +elif default_reference_frame_plate_id == optimised_ref_frame_plate_id: + # 000 is optimised (and 014 is PMAG). + # + # So 000->014 is the existing optimised rotation (000->005). + rotation_time_samples_pmag = clone_rotation_sequence(optimised_rotation_sequence, rotation_description_pmag) +elif default_reference_frame_plate_id == no_net_rotation_ref_frame_plate_id: + # 000 is no-net rotation (and 014 is PMAG). + # + # So 000->014 is the existing no-net rotation (000->005). + rotation_time_samples_pmag = clone_rotation_sequence(no_net_rotation_sequence, rotation_description_pmag) +elif default_reference_frame_plate_id == true_polar_wander_ref_frame_plate_id: + # 000 is TPW (and 014 is PMAG). + # + # So 000->014 is the calculated TPW rotation sequence. + rotation_time_samples_pmag = calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, rotation_description_pmag) +else: + raise ValueError('Default reference frame plate ID is not one of the reference frame plate IDs') + +# Create a new PMAG rotation sequence R(000->014). +pmag_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( + fixed_plate_id=0, + moving_plate_id=pmag_reference_frame_plate_id, + total_reconstruction_pole=pygplates.GpmlIrregularSampling(rotation_time_samples_pmag)) + +# Default rotation model (anchor plate 000). +default_rotation_model = pygplates.RotationModel(output_rotation_feature_collections + [pmag_rotation_feature]) + + +def remove_pmag_from_rotation_sequence(rotation_sequence, rotation_description): + rotation_time_samples = [] + for rotation_sample in rotation_sequence: + time = rotation_sample.get_time() + finite_rotation = rotation_sample.get_value().get_finite_rotation() + # R(000->014). + finite_rotation_014_000 = default_rotation_model.get_rotation(time, pmag_reference_frame_plate_id) + # pid->014 = pid->000 * 000->014 + # pid->000 = pid->014 * inverse[000->014] + finite_rotation = finite_rotation * finite_rotation_014_000.get_inverse() + rotation_time_samples.append( + pygplates.GpmlTimeSample( + pygplates.GpmlFiniteRotation(finite_rotation), + time, + rotation_description, + rotation_sample.is_enabled())) + return rotation_time_samples + +# Generate the rotation sequences for the remaining reference frames (opt, NNR, TPW). +# +# This depends on which reference frame 000 is assigned to. +if default_reference_frame_plate_id == pmag_reference_frame_plate_id: + # Both 000 and 014 are PMAG. + # + # 015->000 is the existing optimised rotation (000->005). + rotation_time_samples_optimised = clone_rotation_sequence(optimised_rotation_sequence, rotation_description_optimised) + # 016->000 is the existing no-net rotation (000->005). + rotation_time_samples_no_net_rotation = clone_rotation_sequence(no_net_rotation_sequence, rotation_description_no_net_rotation) + # 017->000 is the calculated TPW rotation sequence. + rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, rotation_description_true_polar_wander) +elif default_reference_frame_plate_id == optimised_ref_frame_plate_id: + # 000 is optimised (and 014 is PMAG). # - # R(opt_minus_pmag->000) = R(opt->000) * inverse[R(000->701)] - # = inverse[R(000->opt)] * inverse[R(000->701)] - # R(000->opt_minus_pmag)) = inverse[R(opt_minus_pmag->000)] - # = inverse[inverse[R(000->opt)] * inverse[R(000->701)]] - finite_rotation_optimised_minus_pmag = (optimised_finite_rotation.get_inverse() * finite_rotation_701_000.get_inverse()).get_inverse() - - rotation_time_samples_optimised_minus_pmag.append( - pygplates.GpmlTimeSample( - pygplates.GpmlFiniteRotation(finite_rotation_optimised_minus_pmag), - time, - 'optAPM without pmag')) - -# Create a new optimised-minus-pmag rotation sequence. -optimised_minus_pmag_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( - 0, - optimised_minus_pmag_ref_frame_plate_id, - pygplates.GpmlIrregularSampling(rotation_time_samples_optimised_minus_pmag)) - -# Add reference frames to feature collection. -output_ref_frames_feature_collection = pygplates.FeatureCollection() -output_ref_frames_feature_collection.add(optimised_rotation_feature) -output_ref_frames_feature_collection.add(no_net_rotation_feature) -output_ref_frames_feature_collection.add(optimised_minus_pmag_rotation_feature) - -# Write reference frames to output rotation file. -output_ref_frames_feature_collection.write(output_ref_frames_file) + # 015->000 is zero (identity). + rotation_time_samples_optimised = create_identity_rotation_sequence(rotation_description_optimised, max_time_of_sequences_with_fixed_000) + # 016->014 is the existing no-net rotation (000->005). + # 016->014 = 016->000 * 000->014 + # 016->000 = 016->014 * inverse[000->014] + rotation_time_samples_no_net_rotation = remove_pmag_from_rotation_sequence(no_net_rotation_sequence, rotation_description_no_net_rotation) + # 017->014 is the calculated TPW rotation sequence. + # 017->014 = 017->000 * 000->014 + # 017->000 = 017->014 * inverse[000->014] + rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, rotation_description_true_polar_wander) + rotation_time_samples_true_polar_wander = remove_pmag_from_rotation_sequence(rotation_time_samples_true_polar_wander, rotation_description_true_polar_wander) +elif default_reference_frame_plate_id == no_net_rotation_ref_frame_plate_id: + # 000 is no-net rotation (and 014 is PMAG). + # + # 015->014 is the existing optimised rotation (000->005). + # 015->014 = 015->000 * 000->014 + # 015->000 = 015->014 * inverse[000->014] + rotation_time_samples_optimised = remove_pmag_from_rotation_sequence(optimised_rotation_sequence, rotation_description_optimised) + # 016->000 is zero (identity). + rotation_time_samples_no_net_rotation = create_identity_rotation_sequence(rotation_description_no_net_rotation, max_time_of_sequences_with_fixed_000) + # 017->014 is the calculated TPW rotation sequence. + # 017->014 = 017->000 * 000->014 + # 017->000 = 017->014 * inverse[000->014] + rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, rotation_description_true_polar_wander) + rotation_time_samples_true_polar_wander = remove_pmag_from_rotation_sequence(rotation_time_samples_true_polar_wander, rotation_description_true_polar_wander) +elif default_reference_frame_plate_id == true_polar_wander_ref_frame_plate_id: + # 000 is TPW (and 014 is PMAG). + # + # 015->014 is the existing optimised rotation (000->005). + # 015->014 = 015->000 * 000->014 + # 015->000 = 015->014 * inverse[000->014] + rotation_time_samples_optimised = remove_pmag_from_rotation_sequence(optimised_rotation_sequence, rotation_description_optimised) + # 016->014 is the existing no-net rotation (000->005). + # 016->014 = 016->000 * 000->014 + # 016->000 = 016->014 * inverse[000->014] + rotation_time_samples_no_net_rotation = remove_pmag_from_rotation_sequence(no_net_rotation_sequence, rotation_description_no_net_rotation) + # 017->000 is zero (identity). + rotation_time_samples_true_polar_wander = create_identity_rotation_sequence(rotation_description_true_polar_wander, max_time_of_sequences_with_fixed_000) +else: + raise ValueError('Default reference frame plate ID is not one of the reference frame plate IDs') + +# Create a new optimised rotation sequence R(015->000). +optimised_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( + fixed_plate_id=optimised_ref_frame_plate_id, + moving_plate_id=0, + total_reconstruction_pole=pygplates.GpmlIrregularSampling(rotation_time_samples_optimised)) + +# Create a new no-net rotation sequence R(016->000). +no_net_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( + fixed_plate_id=no_net_rotation_ref_frame_plate_id, + moving_plate_id=0, + total_reconstruction_pole=pygplates.GpmlIrregularSampling(rotation_time_samples_no_net_rotation)) + +# Create a new true polar wander rotation sequence R(017->000). +true_polar_wander_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( + fixed_plate_id=true_polar_wander_ref_frame_plate_id, + moving_plate_id=0, + total_reconstruction_pole=pygplates.GpmlIrregularSampling(rotation_time_samples_true_polar_wander)) + +# Add reference frames to a feature collection. +output_reference_frames_feature_collection = pygplates.FeatureCollection() +output_reference_frames_feature_collection.add(pmag_rotation_feature) +output_reference_frames_feature_collection.add(optimised_rotation_feature) +output_reference_frames_feature_collection.add(no_net_rotation_feature) +output_reference_frames_feature_collection.add(true_polar_wander_rotation_feature) + +# Write the reference frames to their associated output rotation file. +output_reference_frames_feature_collection.write(output_reference_frames_filename) + +# Write the output rotation filenames associated with the input rotation files. +for file_index in range(len(output_rotation_feature_collections)): + output_rotation_feature_collection = output_rotation_feature_collections[file_index] + output_rotation_filename = output_rotation_filenames[file_index] + output_rotation_feature_collection.write(output_rotation_filename) From 2540724919fae5a497af9d51fa41bc699db5a176 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Wed, 4 Jun 2025 14:48:18 +1000 Subject: [PATCH 3/7] Reverse order of subtraction: TPW=PMAG-OPT (instead of OPT-PMAG). Because paleomag includes a mantle reference frame and true polar wander (PMAG=OPT+TPW). We were using the old 005-000 which is OPT-PMAG. --- .../add_ref_frames_to_unoptimised.py | 151 ++++++++++++------ 1 file changed, 105 insertions(+), 46 deletions(-) diff --git a/supplementary/add_ref_frames_to_unoptimised.py b/supplementary/add_ref_frames_to_unoptimised.py index 3aca145..6966c30 100644 --- a/supplementary/add_ref_frames_to_unoptimised.py +++ b/supplementary/add_ref_frames_to_unoptimised.py @@ -32,12 +32,12 @@ # The suffix to append to the basenames of the input rotation filenames to get the output rotation filenames. # # Note: The output rotation files also end up in the 'optimisation/' sub-directory. -output_rotation_filenames_suffix = '_20240725_2' +output_rotation_filenames_suffix = '_20240725' # Plate IDs of reference frames. # # These are PMAG, optimised, no-net rotation and true polar wander (which is optimised with PMAG removed). -pmag_reference_frame_plate_id = 14 +pmag_ref_frame_plate_id = 14 optimised_ref_frame_plate_id = 15 no_net_rotation_ref_frame_plate_id = 16 true_polar_wander_ref_frame_plate_id = 17 @@ -75,12 +75,6 @@ ########################## -# Description for the comments of the reference frame rotations. -rotation_description_pmag = ' {:03d} is paleomag reference frame'.format(pmag_reference_frame_plate_id) -rotation_description_optimised = ' {:03d} is optimised mantle reference frame @absage'.format(optimised_ref_frame_plate_id) -rotation_description_no_net_rotation = ' {:03d} is no-net rotation reference frame'.format(no_net_rotation_ref_frame_plate_id) -rotation_description_true_polar_wander = ' {:03d} is true polar wander reference frame (approx)'.format(true_polar_wander_ref_frame_plate_id) - # Each input rotation file will get a corresponding output rotation file. output_rotation_feature_collections = [pygplates.FeatureCollection(filename) for filename in input_rotation_filenames] @@ -101,7 +95,7 @@ ########################## -# Change fixed plate 000 to +# Change fixed plate 000 to # (and remove any existing leftover 005-000 rotation features). max_time_of_sequences_with_fixed_000 = 0.0 for rotation_feature_collection in output_rotation_feature_collections: @@ -115,18 +109,18 @@ def is_feature_005_000(feature): return False rotation_feature_collection.remove(is_feature_005_000) - # Change fixed plate 000 to . + # Change fixed plate 000 to . for rotation_feature in rotation_feature_collection: total_reconstruction_pole = rotation_feature.get_total_reconstruction_pole() if total_reconstruction_pole: fixed_plate_id, moving_plate_id, rotation_sequence = total_reconstruction_pole - # Change fixed plate ID 000 to . - # We want everything that references 000 to now reference . + # Change fixed plate ID 000 to . + # We want everything that references 000 to now reference . if fixed_plate_id == 0: - # Change the fixed plate ID from 000 to . + # Change the fixed plate ID from 000 to . rotation_feature.set_total_reconstruction_pole( - pmag_reference_frame_plate_id, moving_plate_id, rotation_sequence) + pmag_ref_frame_plate_id, moving_plate_id, rotation_sequence) # The oldest time in this sequence. max_time = max(time_sample.get_time() for time_sample in rotation_sequence.get_time_samples()) # The oldest time in all sequences (with fixed plate 000). @@ -139,7 +133,7 @@ def is_feature_005_000(feature): # Note: It has no 000 plate ID. Instead it's anchor plate is the PMAG reference frame plate ID. original_rotation_model_anchored_014 = pygplates.RotationModel( output_rotation_feature_collections, - default_anchor_plate_id=pmag_reference_frame_plate_id) + default_anchor_plate_id=pmag_ref_frame_plate_id) def find_005_000_rotation_sequence(rotation_filename): @@ -181,65 +175,106 @@ def create_identity_rotation_sequence(rotation_description, begin_time, end_time def clone_rotation_sequence(rotation_sequence, rotation_description): rotation_time_samples = [] - for rotation_sample in rotation_sequence.get_enabled_time_samples(): + for rotation_sample in rotation_sequence: rotation_time_samples.append( pygplates.GpmlTimeSample( pygplates.GpmlFiniteRotation(rotation_sample.get_value().get_finite_rotation()), rotation_sample.get_time(), - rotation_description)) + rotation_description, + rotation_sample.is_enabled())) return rotation_time_samples def calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, rotation_description): rotation_time_samples = [] - for optimised_rotation_sample in optimised_rotation_sequence.get_enabled_time_samples(): + for optimised_rotation_sample in optimised_rotation_sequence: time = optimised_rotation_sample.get_time() # R(000->005). optimised_finite_rotation = optimised_rotation_sample.get_value().get_finite_rotation() # R(014->701). finite_rotation_701_014 = original_rotation_model_anchored_014.get_rotation(time, 701) - # We want TPW->701 to be equivalent to the optimised reference frame 000->005: # - # R(TPW->701) = R(000->005) - # R(TPW->014) * R(014->701) = R(000->005) - # R(TPW->014) = R(000->005) * inverse[R(014->701)] - finite_rotation_true_polar_wander = optimised_finite_rotation * finite_rotation_701_014.get_inverse() + # PMAG contain a mantle reference frame and true polar wander (TPW). + # Whereas OPT only contains a mantle reference frame. + # + # TPW + OPT = PMAG + # TPW = PMAG - OPT + # + # So looking at Africa this becomes... + # + # R(PMAG->701) - R(OPT->701) = R(PMAG->701) * R(701->OPT) + # = R(PMAG->OPT) + # + # ...and, using the old optimised rotation R(000->005) (where 000 was OPT and 005 was PMAG)... + # + # R(PMAG->OPT) = R(005->000) + # = inverse[R(000->005)] + # + # ...therefore... + # + # R(PMAG->701) - R(OPT->701) = R(PMAG->OPT) + # = inverse[R(000->005)] + # + # So we want 017->701 to be equivalent to the *inverse* of the optimised rotation 000->005: + # + # TPW = PMAG - OPT + # + # R(017->701) = R(PMAG->701) - R(OPT->701) + # = R(PMAG->OPT) + # = inverse[R(000->005)] + # + # And expressing this relative to 014 becomes... + # + # R(017->701) = inverse[R(000->005)] + # R(017->014) * R(014->701) = inverse[R(000->005)] + # R(017->014) = inverse[R(000->005)] * inverse[R(014->701)] + # + finite_rotation_true_polar_wander = optimised_finite_rotation.get_inverse() * finite_rotation_701_014.get_inverse() rotation_time_samples.append( pygplates.GpmlTimeSample( pygplates.GpmlFiniteRotation(finite_rotation_true_polar_wander), time, - rotation_description)) + rotation_description, + optimised_rotation_sample.is_enabled())) return rotation_time_samples # Generate the rotation sequence for the PMAG reference frame. # # This depends on which reference frame 000 is assigned to. -if default_reference_frame_plate_id == pmag_reference_frame_plate_id: +if default_reference_frame_plate_id == pmag_ref_frame_plate_id: # Both 000 and 014 are PMAG. # # So 000->014 is zero (identity). - rotation_time_samples_pmag = create_identity_rotation_sequence(rotation_description_pmag, max_time_of_sequences_with_fixed_000) + rotation_time_samples_pmag = create_identity_rotation_sequence( + f' Reference frames: paleomag ({pmag_ref_frame_plate_id:03d} and 000)', + max_time_of_sequences_with_fixed_000) elif default_reference_frame_plate_id == optimised_ref_frame_plate_id: # 000 is optimised (and 014 is PMAG). # # So 000->014 is the existing optimised rotation (000->005). - rotation_time_samples_pmag = clone_rotation_sequence(optimised_rotation_sequence, rotation_description_pmag) + rotation_time_samples_pmag = clone_rotation_sequence( + optimised_rotation_sequence, + f' Reference frames: paleomag ({pmag_ref_frame_plate_id:03d}) and optimised mantle (000)') elif default_reference_frame_plate_id == no_net_rotation_ref_frame_plate_id: # 000 is no-net rotation (and 014 is PMAG). # # So 000->014 is the existing no-net rotation (000->005). - rotation_time_samples_pmag = clone_rotation_sequence(no_net_rotation_sequence, rotation_description_pmag) + rotation_time_samples_pmag = clone_rotation_sequence( + no_net_rotation_sequence, + f' Reference frames: paleomag ({pmag_ref_frame_plate_id:03d}) and no-net rotation (000)') elif default_reference_frame_plate_id == true_polar_wander_ref_frame_plate_id: # 000 is TPW (and 014 is PMAG). # # So 000->014 is the calculated TPW rotation sequence. - rotation_time_samples_pmag = calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, rotation_description_pmag) + rotation_time_samples_pmag = calculate_true_polar_wander_rotation_sequence( + optimised_rotation_sequence, + f' Reference frames: paleomag ({pmag_ref_frame_plate_id:03d}) and approx true polar wander (000)') else: raise ValueError('Default reference frame plate ID is not one of the reference frame plate IDs') # Create a new PMAG rotation sequence R(000->014). pmag_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( fixed_plate_id=0, - moving_plate_id=pmag_reference_frame_plate_id, + moving_plate_id=pmag_ref_frame_plate_id, total_reconstruction_pole=pygplates.GpmlIrregularSampling(rotation_time_samples_pmag)) # Default rotation model (anchor plate 000). @@ -252,7 +287,7 @@ def remove_pmag_from_rotation_sequence(rotation_sequence, rotation_description): time = rotation_sample.get_time() finite_rotation = rotation_sample.get_value().get_finite_rotation() # R(000->014). - finite_rotation_014_000 = default_rotation_model.get_rotation(time, pmag_reference_frame_plate_id) + finite_rotation_014_000 = default_rotation_model.get_rotation(time, pmag_ref_frame_plate_id) # pid->014 = pid->000 * 000->014 # pid->000 = pid->014 * inverse[000->014] finite_rotation = finite_rotation * finite_rotation_014_000.get_inverse() @@ -267,56 +302,80 @@ def remove_pmag_from_rotation_sequence(rotation_sequence, rotation_description): # Generate the rotation sequences for the remaining reference frames (opt, NNR, TPW). # # This depends on which reference frame 000 is assigned to. -if default_reference_frame_plate_id == pmag_reference_frame_plate_id: +if default_reference_frame_plate_id == pmag_ref_frame_plate_id: # Both 000 and 014 are PMAG. # # 015->000 is the existing optimised rotation (000->005). - rotation_time_samples_optimised = clone_rotation_sequence(optimised_rotation_sequence, rotation_description_optimised) + rotation_time_samples_optimised = clone_rotation_sequence( + optimised_rotation_sequence, + f' Reference frames: paleomag (000) and optimised mantle ({optimised_ref_frame_plate_id:03d})') # 016->000 is the existing no-net rotation (000->005). - rotation_time_samples_no_net_rotation = clone_rotation_sequence(no_net_rotation_sequence, rotation_description_no_net_rotation) + rotation_time_samples_no_net_rotation = clone_rotation_sequence( + no_net_rotation_sequence, + f' Reference frames: paleomag (000) and no-net rotation ({no_net_rotation_ref_frame_plate_id:03d})') # 017->000 is the calculated TPW rotation sequence. - rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, rotation_description_true_polar_wander) + rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence( + optimised_rotation_sequence, + f' Reference frames: paleomag (000) and approx true polar wander ({true_polar_wander_ref_frame_plate_id:03d})') elif default_reference_frame_plate_id == optimised_ref_frame_plate_id: # 000 is optimised (and 014 is PMAG). # # 015->000 is zero (identity). - rotation_time_samples_optimised = create_identity_rotation_sequence(rotation_description_optimised, max_time_of_sequences_with_fixed_000) + rotation_time_samples_optimised = create_identity_rotation_sequence( + f' Reference frames: optimised mantle (000 and {optimised_ref_frame_plate_id:03d})', + max_time_of_sequences_with_fixed_000) # 016->014 is the existing no-net rotation (000->005). # 016->014 = 016->000 * 000->014 # 016->000 = 016->014 * inverse[000->014] - rotation_time_samples_no_net_rotation = remove_pmag_from_rotation_sequence(no_net_rotation_sequence, rotation_description_no_net_rotation) + rotation_time_samples_no_net_rotation = remove_pmag_from_rotation_sequence( + no_net_rotation_sequence, + f' Reference frames: optimised mantle (000) and no-net rotation ({no_net_rotation_ref_frame_plate_id:03d})') # 017->014 is the calculated TPW rotation sequence. # 017->014 = 017->000 * 000->014 # 017->000 = 017->014 * inverse[000->014] - rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, rotation_description_true_polar_wander) - rotation_time_samples_true_polar_wander = remove_pmag_from_rotation_sequence(rotation_time_samples_true_polar_wander, rotation_description_true_polar_wander) + rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, '') + rotation_time_samples_true_polar_wander = remove_pmag_from_rotation_sequence( + rotation_time_samples_true_polar_wander, + f' Reference frames: optimised mantle (000) and approx true polar wander ({true_polar_wander_ref_frame_plate_id:03d})') elif default_reference_frame_plate_id == no_net_rotation_ref_frame_plate_id: # 000 is no-net rotation (and 014 is PMAG). # # 015->014 is the existing optimised rotation (000->005). # 015->014 = 015->000 * 000->014 # 015->000 = 015->014 * inverse[000->014] - rotation_time_samples_optimised = remove_pmag_from_rotation_sequence(optimised_rotation_sequence, rotation_description_optimised) + rotation_time_samples_optimised = remove_pmag_from_rotation_sequence( + optimised_rotation_sequence, + f' Reference frames: no-net rotation (000) and optimised mantle ({optimised_ref_frame_plate_id:03d})') # 016->000 is zero (identity). - rotation_time_samples_no_net_rotation = create_identity_rotation_sequence(rotation_description_no_net_rotation, max_time_of_sequences_with_fixed_000) + rotation_time_samples_no_net_rotation = create_identity_rotation_sequence( + f' Reference frames: no-net rotation (000 and {no_net_rotation_ref_frame_plate_id:03d})', + max_time_of_sequences_with_fixed_000) # 017->014 is the calculated TPW rotation sequence. # 017->014 = 017->000 * 000->014 # 017->000 = 017->014 * inverse[000->014] - rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, rotation_description_true_polar_wander) - rotation_time_samples_true_polar_wander = remove_pmag_from_rotation_sequence(rotation_time_samples_true_polar_wander, rotation_description_true_polar_wander) + rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, '') + rotation_time_samples_true_polar_wander = remove_pmag_from_rotation_sequence( + rotation_time_samples_true_polar_wander, + f' Reference frames: no-net rotation (000) and approx true polar wander ({true_polar_wander_ref_frame_plate_id:03d})') elif default_reference_frame_plate_id == true_polar_wander_ref_frame_plate_id: # 000 is TPW (and 014 is PMAG). # # 015->014 is the existing optimised rotation (000->005). # 015->014 = 015->000 * 000->014 # 015->000 = 015->014 * inverse[000->014] - rotation_time_samples_optimised = remove_pmag_from_rotation_sequence(optimised_rotation_sequence, rotation_description_optimised) + rotation_time_samples_optimised = remove_pmag_from_rotation_sequence( + optimised_rotation_sequence, + f' Reference frames: approx true polar wander (000) and optimised mantle ({optimised_ref_frame_plate_id:03d})') # 016->014 is the existing no-net rotation (000->005). # 016->014 = 016->000 * 000->014 # 016->000 = 016->014 * inverse[000->014] - rotation_time_samples_no_net_rotation = remove_pmag_from_rotation_sequence(no_net_rotation_sequence, rotation_description_no_net_rotation) + rotation_time_samples_no_net_rotation = remove_pmag_from_rotation_sequence( + no_net_rotation_sequence, + f' Reference frames: approx true polar wander (000) and no-net rotation ({no_net_rotation_ref_frame_plate_id:03d})') # 017->000 is zero (identity). - rotation_time_samples_true_polar_wander = create_identity_rotation_sequence(rotation_description_true_polar_wander, max_time_of_sequences_with_fixed_000) + rotation_time_samples_true_polar_wander = create_identity_rotation_sequence( + f' Reference frames: approx true polar wander (000 and {true_polar_wander_ref_frame_plate_id:03d})', + max_time_of_sequences_with_fixed_000) else: raise ValueError('Default reference frame plate ID is not one of the reference frame plate IDs') From 80de38fda1e963207eed63b8a44f42f9f5990e44 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Wed, 4 Jun 2025 15:01:02 +1000 Subject: [PATCH 4/7] Use rotation sequences 015-000, 016-000, 017-000 (instead of 000-*). MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit You can’t have 000-015, 000-016, 000-017 rotation sequences (ie, with 000 as a shared moving plate) since they’ll overlap in time (which GPlates complains about). This was causing issues when anchoring to 015 and 016 (but 017 was fine). --- .../add_ref_frames_to_unoptimised.py | 66 +++++++++++++------ 1 file changed, 45 insertions(+), 21 deletions(-) diff --git a/supplementary/add_ref_frames_to_unoptimised.py b/supplementary/add_ref_frames_to_unoptimised.py index 6966c30..80967ac 100644 --- a/supplementary/add_ref_frames_to_unoptimised.py +++ b/supplementary/add_ref_frames_to_unoptimised.py @@ -308,35 +308,35 @@ def remove_pmag_from_rotation_sequence(rotation_sequence, rotation_description): # 015->000 is the existing optimised rotation (000->005). rotation_time_samples_optimised = clone_rotation_sequence( optimised_rotation_sequence, - f' Reference frames: paleomag (000) and optimised mantle ({optimised_ref_frame_plate_id:03d})') + f' Reference frames: optimised mantle ({optimised_ref_frame_plate_id:03d}) and paleomag (000)') # 016->000 is the existing no-net rotation (000->005). rotation_time_samples_no_net_rotation = clone_rotation_sequence( no_net_rotation_sequence, - f' Reference frames: paleomag (000) and no-net rotation ({no_net_rotation_ref_frame_plate_id:03d})') + f' Reference frames: no-net rotation ({no_net_rotation_ref_frame_plate_id:03d}) and paleomag (000)') # 017->000 is the calculated TPW rotation sequence. rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence( optimised_rotation_sequence, - f' Reference frames: paleomag (000) and approx true polar wander ({true_polar_wander_ref_frame_plate_id:03d})') + f' Reference frames: approx true polar wander ({true_polar_wander_ref_frame_plate_id:03d}) and paleomag (000)') elif default_reference_frame_plate_id == optimised_ref_frame_plate_id: # 000 is optimised (and 014 is PMAG). # # 015->000 is zero (identity). rotation_time_samples_optimised = create_identity_rotation_sequence( - f' Reference frames: optimised mantle (000 and {optimised_ref_frame_plate_id:03d})', + f' Reference frames: optimised mantle ({optimised_ref_frame_plate_id:03d} and 000)', max_time_of_sequences_with_fixed_000) # 016->014 is the existing no-net rotation (000->005). # 016->014 = 016->000 * 000->014 # 016->000 = 016->014 * inverse[000->014] rotation_time_samples_no_net_rotation = remove_pmag_from_rotation_sequence( no_net_rotation_sequence, - f' Reference frames: optimised mantle (000) and no-net rotation ({no_net_rotation_ref_frame_plate_id:03d})') + f' Reference frames: no-net rotation ({no_net_rotation_ref_frame_plate_id:03d}) and optimised mantle (000)') # 017->014 is the calculated TPW rotation sequence. # 017->014 = 017->000 * 000->014 # 017->000 = 017->014 * inverse[000->014] rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, '') rotation_time_samples_true_polar_wander = remove_pmag_from_rotation_sequence( rotation_time_samples_true_polar_wander, - f' Reference frames: optimised mantle (000) and approx true polar wander ({true_polar_wander_ref_frame_plate_id:03d})') + f' Reference frames: approx true polar wander ({true_polar_wander_ref_frame_plate_id:03d}) and optimised mantle (000)') elif default_reference_frame_plate_id == no_net_rotation_ref_frame_plate_id: # 000 is no-net rotation (and 014 is PMAG). # @@ -345,10 +345,10 @@ def remove_pmag_from_rotation_sequence(rotation_sequence, rotation_description): # 015->000 = 015->014 * inverse[000->014] rotation_time_samples_optimised = remove_pmag_from_rotation_sequence( optimised_rotation_sequence, - f' Reference frames: no-net rotation (000) and optimised mantle ({optimised_ref_frame_plate_id:03d})') + f' Reference frames: optimised mantle ({optimised_ref_frame_plate_id:03d}) and no-net rotation (000)') # 016->000 is zero (identity). rotation_time_samples_no_net_rotation = create_identity_rotation_sequence( - f' Reference frames: no-net rotation (000 and {no_net_rotation_ref_frame_plate_id:03d})', + f' Reference frames: no-net rotation ({no_net_rotation_ref_frame_plate_id:03d} and 000)', max_time_of_sequences_with_fixed_000) # 017->014 is the calculated TPW rotation sequence. # 017->014 = 017->000 * 000->014 @@ -356,7 +356,7 @@ def remove_pmag_from_rotation_sequence(rotation_sequence, rotation_description): rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, '') rotation_time_samples_true_polar_wander = remove_pmag_from_rotation_sequence( rotation_time_samples_true_polar_wander, - f' Reference frames: no-net rotation (000) and approx true polar wander ({true_polar_wander_ref_frame_plate_id:03d})') + f' Reference frames: approx true polar wander ({true_polar_wander_ref_frame_plate_id:03d}) and no-net rotation (000)') elif default_reference_frame_plate_id == true_polar_wander_ref_frame_plate_id: # 000 is TPW (and 014 is PMAG). # @@ -365,37 +365,61 @@ def remove_pmag_from_rotation_sequence(rotation_sequence, rotation_description): # 015->000 = 015->014 * inverse[000->014] rotation_time_samples_optimised = remove_pmag_from_rotation_sequence( optimised_rotation_sequence, - f' Reference frames: approx true polar wander (000) and optimised mantle ({optimised_ref_frame_plate_id:03d})') + f' Reference frames: optimised mantle ({optimised_ref_frame_plate_id:03d}) and approx true polar wander (000)') # 016->014 is the existing no-net rotation (000->005). # 016->014 = 016->000 * 000->014 # 016->000 = 016->014 * inverse[000->014] rotation_time_samples_no_net_rotation = remove_pmag_from_rotation_sequence( no_net_rotation_sequence, - f' Reference frames: approx true polar wander (000) and no-net rotation ({no_net_rotation_ref_frame_plate_id:03d})') + f' Reference frames: no-net rotation ({no_net_rotation_ref_frame_plate_id:03d}) and approx true polar wander (000)') # 017->000 is zero (identity). rotation_time_samples_true_polar_wander = create_identity_rotation_sequence( - f' Reference frames: approx true polar wander (000 and {true_polar_wander_ref_frame_plate_id:03d})', + f' Reference frames: approx true polar wander ({true_polar_wander_ref_frame_plate_id:03d} and 000)', max_time_of_sequences_with_fixed_000) else: raise ValueError('Default reference frame plate ID is not one of the reference frame plate IDs') + +def invert_rotation_sequence(rotation_sequence): + rotation_time_samples = [] + for rotation_sample in rotation_sequence: + rotation_time_samples.append( + pygplates.GpmlTimeSample( + pygplates.GpmlFiniteRotation(rotation_sample.get_value().get_finite_rotation().get_inverse()), + rotation_sample.get_time(), + rotation_sample.get_description(), + rotation_sample.is_enabled())) + return rotation_time_samples + # Create a new optimised rotation sequence R(015->000). +# +# But we write it out as R(000->015), which is inverse of R(015->000), otherwise 000 is the moving plate ID and it +# will overlap in time with other sequences where 000 is also the moving plate ID (when loading into GPlates/pyGPlates). optimised_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( - fixed_plate_id=optimised_ref_frame_plate_id, - moving_plate_id=0, - total_reconstruction_pole=pygplates.GpmlIrregularSampling(rotation_time_samples_optimised)) + fixed_plate_id=0, + moving_plate_id=optimised_ref_frame_plate_id, + total_reconstruction_pole=pygplates.GpmlIrregularSampling( + invert_rotation_sequence(rotation_time_samples_optimised))) # Create a new no-net rotation sequence R(016->000). +# +# But we write it out as R(000->016), which is inverse of R(016->000), otherwise 000 is the moving plate ID and it +# will overlap in time with other sequences where 000 is also the moving plate ID (when loading into GPlates/pyGPlates). no_net_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( - fixed_plate_id=no_net_rotation_ref_frame_plate_id, - moving_plate_id=0, - total_reconstruction_pole=pygplates.GpmlIrregularSampling(rotation_time_samples_no_net_rotation)) + fixed_plate_id=0, + moving_plate_id=no_net_rotation_ref_frame_plate_id, + total_reconstruction_pole=pygplates.GpmlIrregularSampling( + invert_rotation_sequence(rotation_time_samples_no_net_rotation))) # Create a new true polar wander rotation sequence R(017->000). +# +# But we write it out as R(000->016), which is inverse of R(016->000), otherwise 000 is the moving plate ID and it +# will overlap in time with other sequences where 000 is also the moving plate ID (when loading into GPlates/pyGPlates). true_polar_wander_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( - fixed_plate_id=true_polar_wander_ref_frame_plate_id, - moving_plate_id=0, - total_reconstruction_pole=pygplates.GpmlIrregularSampling(rotation_time_samples_true_polar_wander)) + fixed_plate_id=0, + moving_plate_id=true_polar_wander_ref_frame_plate_id, + total_reconstruction_pole=pygplates.GpmlIrregularSampling( + invert_rotation_sequence(rotation_time_samples_true_polar_wander))) # Add reference frames to a feature collection. output_reference_frames_feature_collection = pygplates.FeatureCollection() From 7dfe1dd47f6df66ce24216b6483cba3f2023f4ee Mon Sep 17 00:00:00 2001 From: John Cannon Date: Wed, 4 Jun 2025 17:32:27 +1000 Subject: [PATCH 5/7] Store reference frames in one of the output feature collections/files. Instead of in a separate 'reference_frames.rot' file. Otherwise users might have thought they could exclude it if they didn't want the reference frames because only using 000 (but 000 is part of the reference frames and so always needed inclusion). --- .../add_ref_frames_to_unoptimised.py | 22 ++++++++++++++----- 1 file changed, 16 insertions(+), 6 deletions(-) diff --git a/supplementary/add_ref_frames_to_unoptimised.py b/supplementary/add_ref_frames_to_unoptimised.py index 80967ac..20e507f 100644 --- a/supplementary/add_ref_frames_to_unoptimised.py +++ b/supplementary/add_ref_frames_to_unoptimised.py @@ -89,9 +89,6 @@ output_filename = os.path.join('optimisation', output_filename) output_rotation_filenames.append(output_filename) -# The output rotation filename to store the optimised, NNR and true polar wander reference frames. -output_reference_frames_filename = 'optimisation/reference_frames{}.rot'.format(output_rotation_filenames_suffix) - ########################## @@ -421,15 +418,28 @@ def invert_rotation_sequence(rotation_sequence): total_reconstruction_pole=pygplates.GpmlIrregularSampling( invert_rotation_sequence(rotation_time_samples_true_polar_wander))) -# Add reference frames to a feature collection. +# Add the reference frame features to a feature collection. output_reference_frames_feature_collection = pygplates.FeatureCollection() output_reference_frames_feature_collection.add(pmag_rotation_feature) output_reference_frames_feature_collection.add(optimised_rotation_feature) output_reference_frames_feature_collection.add(no_net_rotation_feature) output_reference_frames_feature_collection.add(true_polar_wander_rotation_feature) -# Write the reference frames to their associated output rotation file. -output_reference_frames_feature_collection.write(output_reference_frames_filename) +# Store the reference frames in one of the output feature collections/files (we arbitrarily choose the first). +# +# Previously we added to a separate file called 'reference_frames.rot'. +# But that can be confusing for novice users who might think they can exclude that file if +# they don't want any of the reference frames (and hence just rely on the default 000 plate ID). +# But that would have been erroneous because 000 only existed in 'reference_frames.rot'. +# So now the reference frames are included within the normal output rotation files (associated with the input files). +# +# Note: We could split the rotation sequences of the reference frames across the output rotation files +# (instead of dumping them all in just one of the files). For example, 0-1000 Ma into '1000_0_rotfile.rot' +# and 1000-1800 Ma into '1800_1000_rotfile.rot'. But we only have the filename to determine the time range +# of each rotation file. And that information isn't always in the filename. +# Besides, the original rotation sequences aren't always split up exactly according to the filename anyway. +output_reference_frames_feature_collection.add(output_rotation_feature_collections[0]) +output_rotation_feature_collections[0] = output_reference_frames_feature_collection # Write the output rotation filenames associated with the input rotation files. for file_index in range(len(output_rotation_feature_collections)): From 2d8c833861a8351a6ba049c279fa36e06424a846 Mon Sep 17 00:00:00 2001 From: John Cannon Date: Fri, 26 Jun 2026 23:45:16 +1000 Subject: [PATCH 6/7] Move add_ref_frames_to_unoptimised.py into main directory. No longer really need a supplementary directory. --- ...unoptimised.py => add_ref_frames_to_unoptimised.py | 0 supplementary/Readme.txt | 11 ----------- 2 files changed, 11 deletions(-) rename supplementary/add_ref_frames_to_unoptimised.py => add_ref_frames_to_unoptimised.py (100%) delete mode 100644 supplementary/Readme.txt diff --git a/supplementary/add_ref_frames_to_unoptimised.py b/add_ref_frames_to_unoptimised.py similarity index 100% rename from supplementary/add_ref_frames_to_unoptimised.py rename to add_ref_frames_to_unoptimised.py diff --git a/supplementary/Readme.txt b/supplementary/Readme.txt deleted file mode 100644 index d3e4f5f..0000000 --- a/supplementary/Readme.txt +++ /dev/null @@ -1,11 +0,0 @@ -This folder contains supplementary scripts that are not part of the main optimisation workflow. - -add_ref_frames_to_unoptimised.py --------------------------------- - -This script reads the original un-optimised rotation files, and the optimised and -the no-net rotation reference frames, and generates a new set of rotation files -with various reference frames added. -Currently the added reference frames are the original paleomag (PMAG), the optimised mantle, -the no-net rotation, and an approximation of true polar wander (TPW) calculated by replacing -the PMAG reference frame with the optimised mantle reference frame. \ No newline at end of file From 50c8b0b931fdb3152b3d57a34cbcf8b70876bdfe Mon Sep 17 00:00:00 2001 From: John Cannon Date: Fri, 26 Jun 2026 23:57:32 +1000 Subject: [PATCH 7/7] Add Merdith rev6-2 model config, auto-derive paths in add_ref_frames_to_unoptimised.py, and add optimised_max_NNR frame (016). - Optimised_config.py: add Merdith_etal_plate_model_1Ga-present_rev6-2 (rotation file, topologies, continental polygons, diagnostics) and fix Zahirovic_etal_2022_GDJ continental polygons path and reference params (seeds from paleomagnetic pole via Africa 701701 / GAPWAP). - add_ref_frames_to_unoptimised.py: refactor to read all input/output paths from Optimised_config.py automatically (no manual editing required); add optimised_max_NNR reference frame (plate ID 016, from the nr_max run), shifting NNR to 017 and TPW to 018; add unoptimised_input_is_in_pmag_frame guard that aborts rather than silently mislabelling frames; extract helper functions for the frame calculations; add a module docstring. - README.md: rewrite to reflect current capabilities (laptop runtime, GPlately PMM, 5 reference frames, updated installation and running instructions). Co-Authored-By: Claude Sonnet 4.6 --- Optimised_config.py | 81 +++- README.md | 490 +++++++++++----------- add_ref_frames_to_unoptimised.py | 671 ++++++++++++++----------------- 3 files changed, 619 insertions(+), 623 deletions(-) diff --git a/Optimised_config.py b/Optimised_config.py index cf9f604..f06f05b 100644 --- a/Optimised_config.py +++ b/Optimised_config.py @@ -30,7 +30,14 @@ #data_model = 'Global_Model_WD_Internal_Release_2019_v2' #data_model = 'Global_Model_WD_Internal_Release-EarthBytePlateMotionModel-TRUNK' #data_model = 'SM2-Merdith_et_al_1_Ga_reconstruction_v1.1' -data_model = 'Global_2000-540' +#data_model = 'Global_2000-540' +#data_model = 'Merdith_etal_plate_model_1Ga-present_rev6-2' +# Zahirovic et al. (2022) GDJ model (0-410 Ma), downloaded via the GPlately Plate Model Manager +# (PMM model name "Zahirovic2022") into 'data/Zahirovic_etal_2022_GDJ/'. The reference plate is +# Africa (701); each interval's search is seeded from the model's paleomagnetic ("pmag") pole, which +# the model stores as plate 701701 (GAPWAP; readme.txt: "Set the anchor plate ID to 701701 to use +# PMAG reference frame"). See get_reference_params() below. +data_model = 'Zahirovic_etal_2022_GDJ' # The model name is suffixed to various output filenames. @@ -44,6 +51,8 @@ model_name = "20240725_run3" elif data_model.startswith('SM2-Merdith_et_al_1_Ga_reconstruction'): model_name = "git_20231114_run1" +elif data_model == 'Merdith_etal_plate_model_1Ga-present_rev6-2': + model_name = "rev6-2_run2" elif data_model == "Global_2000-540": model_name = "20250410_90W_run1" else: @@ -168,6 +177,23 @@ ( 'EDRG_90W_2000-540Ma.rot', )] +elif data_model == 'Merdith_etal_plate_model_1Ga-present_rev6-2': + # Explicitly use the main rotation file (the directory also contains a 'SH2017' variant + # which must not be loaded at the same time). + original_rotation_filenames = [os.path.join(data_model, rot_file) for rot_file in + ( + '1000_0_rotfile_MER21.rot', + )] +elif data_model == 'Zahirovic_etal_2022_GDJ': + # Single combined rotation file from the PMM download. It lives in a 'Rotations/' sub-directory + # (so it must be listed explicitly rather than globbed from the data model's top level). + # This one file contains both reference frames: plate 701 (rel 000) is the Muller++ 2019 + # optimised mantle frame, and plate 701701 (rel 701) carries the GAPWAP paleomagnetic + # correction (Merdith et al. 2021) - anchoring 701701 gives the "pmag" reference frame. + original_rotation_filenames = [os.path.join(data_model, rot_file) for rot_file in + ( + os.path.join('Rotations', 'CombinedRotations.rot'), + )] # 2) Automatically gather all '.rot' files (and make filenames relative to the 'data/' directory). else: original_rotation_filenames = [os.path.relpath(abs_path, datadir) for abs_path in @@ -210,6 +236,29 @@ 'EDRG_topology_90W_2000-540Ma.gpml', 'EDRG_boundary_90W_2000-540Ma.gpml', )] +elif data_model == 'Merdith_etal_plate_model_1Ga-present_rev6-2': + topology_filenames = [os.path.join(data_model, rot_file) for rot_file in + ( + '250-0_plate_boundaries_MER21.gpml', + '410-250_plate_boundaries_MER21.gpml', + '1000-410-Convergence_MER21.gpml', + '1000-410-Divergence_MER21.gpml', + '1000-410-Topologies_MER21.gpml', + '1000-410-Transforms_MER21.gpml', + 'TopologyBuildingBlocks_MER21.gpml', + )] +elif data_model == 'Zahirovic_etal_2022_GDJ': + # Topology files from the PMM download (in a 'Topologies/' sub-directory). 'Feature_Geometries' + # holds the line geometries referenced by the topological boundaries/networks, so it must be + # loaded alongside them. Only time-appropriate features resolve at each reconstruction time, so + # including the inactive deforming networks is harmless. + topology_filenames = [os.path.join(data_model, 'Topologies', gpml_file) for gpml_file in + ( + 'Plate_Boundaries.gpml', + 'Deforming_Networks_Active.gpml', + 'Deforming_Networks_Inactive.gpml', + 'Feature_Geometries.gpml', + )] # 2) Automatically gather all '.gpml' and '.gpmlz' files (and make filenames relative to the 'data/' directory). else: topology_filenames = [os.path.relpath(abs_path, datadir) for abs_path in @@ -224,11 +273,14 @@ data_model.startswith('SM2-Merdith_et_al_1_Ga_reconstruction')): plate_velocity_continental_polygons_file = os.path.join(data_model, 'shapes_continents_Merdith_et_al.gpml') elif data_model == 'Zahirovic_etal_2022_GDJ': - plate_velocity_continental_polygons_file = os.path.join(data_model, 'StaticGeometries', 'ContinentalPolygons', 'Global_EarthByte_GPlates_PresentDay_ContinentalPolygons.shp') + # Continental polygons from the PMM download (in a 'ContinentalPolygons/' sub-directory). + plate_velocity_continental_polygons_file = os.path.join(data_model, 'ContinentalPolygons', 'Global_EarthByte_GPlates_PresentDay_ContinentalPolygons.shp') elif '1.8ga' in data_model.lower(): plate_velocity_continental_polygons_file = os.path.join(data_model, 'shapes_continents.gpmlz') elif data_model == "Global_2000-540": plate_velocity_continental_polygons_file = os.path.join(data_model, 'Continental_outlines.shp') +elif data_model == 'Merdith_etal_plate_model_1Ga-present_rev6-2': + plate_velocity_continental_polygons_file = os.path.join(data_model, 'COB_MER21_land_mask.shp') else: plate_velocity_continental_polygons_file = None @@ -267,6 +319,7 @@ def plate_velocity_continental_fragmentation_area_threshold_steradians(time): def plate_velocity_continental_fragmentation_gap_threshold_radians(time): if (data_model == 'Global_1000-0_Model_2017' or data_model.startswith('SM2-Merdith_et_al_1_Ga_reconstruction') or + data_model.startswith('Merdith_etal_plate_model') or data_model.startswith('Global_Model_WD_Internal_Release') or data_model == 'Zahirovic_etal_2022_GDJ' or '1.8ga' in data_model.lower() or @@ -285,6 +338,13 @@ def plate_velocity_continental_fragmentation_gap_threshold_radians(time): ridge_file = os.path.join(data_model, 'StaticGeometries', 'AgeGridInput', 'Global_EarthByte_GeeK07_Ridges.gpml') isochron_file = os.path.join(data_model, 'StaticGeometries', 'AgeGridInput', 'Global_EarthByte_GeeK07_Isochrons.gpml') isocob_file = os.path.join(data_model, 'StaticGeometries', 'AgeGridInput', 'Global_EarthByte_GeeK07_IsoCOB.gpml') +elif data_model == 'Merdith_etal_plate_model_1Ga-present_rev6-2': + # There are no static geometries (besides coastlines) for this data model. + # These files are only used when fracture zones are enabled (they're currently disabled), + # so point at files that exist in the root 'data/' directory (so the data loader doesn't fail). + ridge_file = 'Global_EarthByte_230-0Ma_GK07_AREPS_Ridges.gpml' + isochron_file = 'Global_EarthByte_230-0Ma_GK07_AREPS_Isochrons.gpmlz' + isocob_file = 'Global_EarthByte_230-0Ma_GK07_AREPS_IsoCOB.gpml' elif (data_model == 'Global_1000-0_Model_2017' or data_model.startswith('SM2-Merdith_et_al_1_Ga_reconstruction') or data_model == 'Zahirovic_etal_2022_GDJ' or @@ -350,6 +410,7 @@ def cost_function(PTLong1, PTLat1, PTangle1, SPLong, SPLat, SPangle, SPLong_NNR, if (data_model == 'Global_1000-0_Model_2017' or data_model.startswith('SM2-Merdith_et_al_1_Ga_reconstruction') or + data_model.startswith('Merdith_etal_plate_model') or data_model.startswith('Global_Model_WD_Internal_Release') or data_model == 'Zahirovic_etal_2022_GDJ' or '1.8ga' in data_model.lower() or @@ -482,6 +543,7 @@ def cost_function(trench_vel, trench_obl, tm_vel_orth, tm_mean_vel_orth, tm_mean if (data_model == 'Global_1000-0_Model_2017' or data_model.startswith('SM2-Merdith_et_al_1_Ga_reconstruction') or + data_model.startswith('Merdith_etal_plate_model') or data_model.startswith('Global_Model_WD_Internal_Release') or data_model == 'Zahirovic_etal_2022_GDJ' or '1.8ga' in data_model.lower() or @@ -568,6 +630,7 @@ def cost_function(plate_features_are_topologies, velocity_vectors_in_contours, r if (data_model == 'Global_1000-0_Model_2017' or data_model.startswith('SM2-Merdith_et_al_1_Ga_reconstruction') or + data_model.startswith('Merdith_etal_plate_model') or data_model.startswith('Global_Model_WD_Internal_Release') or data_model == 'Zahirovic_etal_2022_GDJ' or '1.8ga' in data_model.lower() or @@ -592,8 +655,22 @@ def get_reference_params(age): If reference rotation filename is None then it means the no-net-rotation model should be used. """ + # + # Zahirovic et al. (2022): reference plate is Africa (701, as for all the other models), but + # each interval's search is *seeded* from the model's paleomagnetic ("pmag") pole rather than + # from the previously-optimised interval. The model stores its pmag path as plate 701701 + # (GAPWAP, Merdith et al. 2021; readme.txt: "Set the anchor plate ID to 701701 to use PMAG + # reference frame"). That path - the rotation of 701701 rel 000 - has been written out as a + # 701-rel-000 reference rotation file (generated from the model's own 701701 chain) so that + # the auto-calc seed (optapm.py: get_rotation(start_age, ref_plate=701, anchor=0)) returns + # the pmag pole. The optimisation itself still zeroes and re-optimises plate 701 in the main + # model; 'ref_rotation_file' only sets the search seed. + # + if data_model == 'Zahirovic_etal_2022_GDJ': + return 701, os.path.join(data_model, 'pmag', 'Zahirovic2022_Palaeomagnetic_Africa.rot') if (data_model == 'Global_1000-0_Model_2017' or data_model.startswith('SM2-Merdith_et_al_1_Ga_reconstruction') or + data_model.startswith('Merdith_etal_plate_model') or data_model.startswith('Global_Model_WD_Internal_Release') or data_model == 'Zahirovic_etal_2022_GDJ' or '1.8ga' in data_model.lower() or diff --git a/README.md b/README.md index 4dc2265..3079981 100644 --- a/README.md +++ b/README.md @@ -1,364 +1,370 @@ # optAPM - Absolute Plate Motion Optimisation -This code adapts Mike Tetley's _optAPM_ code ([Tetley et al., 2019](https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2019JB017442)) to work on a High Performance Computing (HPC) cluster such as Artemis at the University of Sydney. +optAPM optimises the **absolute plate motion** (the mantle reference frame) of an existing +relative plate model. It adapts Mike Tetley's original _optAPM_ code +([Tetley et al., 2019](https://doi.org/10.1029/2019JB017442)) and has since been substantially +reworked: it now builds on [GPlately](https://github.com/GPlates/gplately), runs ~10x faster per +timestep, produces an uncertainty envelope of reference frames rather than a single solution, and +finishes by assembling one rotation file that contains **all** of those reference frames. -## Prerequisites +> **No HPC required.** Thanks to the seed-screening speedup (see below), a full 1000-0 Ma +> optimisation of a 1 Ga global plate model now runs in roughly **4-4.5 hours on a recent laptop** +> (measured: 4 h 24 m on an Apple M3). An HPC allocation is no longer a hard requirement for single +> runs - the workflow still scales to a cluster via MPI when you want it to. -This workflow optimizes the absolute plate motion of an existing plate model. So the data of the plate model needs to be copied into a subdirectory of the `data` directory. +## How the optimisation works -The data for the deforming model can be obtained [here](https://www.earthbyte.org/webdav/ftp/Data_Collections/Muller_etal_2019_Tectonics/Muller_etal_2019_PlateMotionModel/). -The data for the 1Ga model can be obtained [here](https://www.earthbyte.org/webdav/ftp/Data_Collections/Muller_etal_2022_SE/). +The workflow perturbs the absolute rotation (pole and angle) of a reference plate (Africa, 701) and +scores each candidate with an objective function built from several geodynamic penalties: -## Installation - -### Dependencies +* **Net rotation (NR)** of the lithosphere, +* **Trench migration (TM)** - trench-normal velocities at subduction zones, +* **Plate velocity (PV)** - continental speeds, and +* **Hotspot trails (HS)** - only for 0-80 Ma. -The following Python packages are required: +Each penalty is scaled by a weight (and an internal normalisation so the components are comparable at +weight 1.0), then summed into a single cost. NLopt (COBYLA) drives the cost toward a local minimum +from each starting seed, and many seeds are used to find the global minimum at each timestep. The +model is optimised one timestep at a time from present day back to the model's start age. -* gplately>=1.3 +The optimised reference frame for each timestep is stored as a single 005-000 rotation; everything +that originally referenced plate 000 is re-pointed to plate 005, and the optimised 005-000 sequence +carries the absolute plate motion. - > __Note:__ PlateTectonicTools is no longer required - its functionality is now fully integrated into [GPlately](https://github.com/GPlates/gplately) (as the `gplately.ptt` module). +## Installation -* pygplates>=1.0 -* numpy -* scipy -* scikit-image -* pandas -* pmagpy -* xlrd==1.2.0 (apparently required for Excel support; and using >= 2.0 results in an error since 2.0 removed some formats) -* NLopt -* mpi4py -* future -* cartopy -* matplotlib -* ipyparallel -* openpyxl +### Dependencies -### Install on desktop +* `gplately>=1.3` (PlateTectonicTools is **no longer required** - its functionality now lives in + `gplately.ptt`) +* `pygplates>=1.0` +* `numpy`, `scipy`, `scikit-image`, `pandas` +* `pmagpy` +* `xlrd==1.2.0` (Excel support; `>=2.0` removes formats this workflow needs) +* `NLopt`, `mpi4py`, `future`, `cartopy`, `matplotlib`, `ipyparallel`, `openpyxl` -You can install the dependencies using `conda` and `pip`. +Verified against gplately 2.0.0 / pygplates 1.0.0. -First install dependencies that are available on `conda` (in the _conda-forge_ channel): +### Install on desktop / laptop ``` -conda create -n -c conda-forge \ +conda create -n optAPM -c conda-forge \ "pygplates>=1.0" "gplately>=1.3" numpy scipy scikit-image pandas xlrd==1.2.0 NLopt \ future cartopy matplotlib ipyparallel openpyxl -``` - -Then activate the conda environment: - -``` -conda activate -``` - -...where `` should be replaced with the name of your conda environment (eg, `optAPM`). - -On desktop systems we can also use conda to install `mpi4py` (into the conda environment): - -``` +conda activate optAPM conda install mpi4py -``` - -Then use `pip` to install `pmagpy` (into the conda environment): - -``` conda install pip pip install pmagpy ``` -### Install on a HPC system - -Installation on a High Performance Computing (HPC) system can also be done with a local installation of `conda` (and `pip`). However the exception, compared with installing on desktop, is `mpi4py`. It will likely need to be installed differently to ensure that the MPI implementation of the HPC system is used (instead of conda's MPI). +`mpi4py` from conda is fine on a desktop/laptop and lets you use all your cores with +`mpiexec -n ...`. -The example [job submission script](Optimised_APM.sh) works on [NCI's Gadi](https://nci.org.au/our-systems/hpc-systems) HPC. The script assumes that [Miniconda](https://docs.conda.io/en/main/miniconda.html) has been installed (in your $HOME directory by default). +### Install on an HPC system (optional) -> Note: Miniconda currently requires an operating system based on CentOS 7 or above. For example, Gadi is based on CentOS 8 (so it's fine), however the University of Sydney's Artemis HPC is based on CentOS 6 (so it's not). - -> Note: Miniconda might consume a fair amount of your $HOME directory allocation. For example, on Gadi it consumes 4-5GB of my max 10GB allocation. So it may be worth trying to install to a shared project directory instead and also creating the `` there (eg, using `conda create -p ...`). This would require users in the project to add that shared path to their [conda config](https://conda.io/projects/conda/en/latest/user-guide/configuration/use-condarc.html#specify-environment-directories-envs-dirs). And it relies on users not accidentally modifying that conda environment. - -Similarly to installing on desktop, start by creating a conda environment: +Installation on an HPC system is the same, except `mpi4py` should be built against the system MPI so +the cluster's interconnect is used: ``` -conda create -n -c conda-forge \ +conda create -n optAPM -c conda-forge \ "pygplates>=1.0" "gplately>=1.3" numpy scipy scikit-image pandas xlrd==1.2.0 NLopt \ future cartopy matplotlib ipyparallel openpyxl -``` - -Then activate the conda environment: - -``` -conda activate -``` - -...where `` should be replaced with the name of your conda environment (eg, `optAPM`). - -Then load the HPC system's MPI. For example: - -``` -module load openmpi -``` - -Then use `pip` to compile `mpi4py` using the system MPI with: - -``` +conda activate optAPM +module load openmpi # load the system MPI conda install pip -pip install mpi4py -``` - -Then, similarly to installing on desktop, use `pip` to install `pmagpy` (into the conda environment): - -``` +pip install mpi4py # compile against the system MPI pip install pmagpy ``` -#### Job submission script - -In our example [job submission script](Optimised_APM.sh) (that runs on [NCI's Gadi](https://nci.org.au/our-systems/hpc-systems)) we have the following commands (after the various PBS directives) that: - -- Load the system's MPI, -- configure the shell to use `conda activate`, -- activate our previously created conda environnment named `optAPM`, -- run our Python MPI program `Optimised_APM.py` across the CPUs specified in a previous PBS directive. - -``` -#!/bin/bash - -... +The example PBS job script [`Optimised_APM.sh`](Optimised_APM.sh) works on +[NCI's Gadi](https://nci.org.au/our-systems/hpc-systems). See the comments in that script (it sets +`LD_PRELOAD=libmpi.so` to avoid an OpenMPI hierarchy error). -# Use the system MPI implementation (not conda's MPI). -module load openmpi +## The data (relative plate model) -# Initialise the shell for conda environments to avoid the error: -# "CommandNotFoundError: Your shell has not been properly configured to use 'conda activate'." -source ~/.bashrc +optAPM optimises an existing relative plate model, so its data (rotation files, topology files, +continental polygons) must be placed in a sub-directory of `data/` whose name matches the `data_model` +setting in `Optimised_config.py`. -# Activate the "optAPM" conda environment in our home Miniconda installation. -conda activate optAPM +You can obtain plate models from the +[EarthByte data collections](https://www.earthbyte.org/category/resources/data-models/) or, more +conveniently, with the [GPlately Plate Model Manager (PMM)](https://github.com/GPlates/gplately), +which downloads a published model into a tidy directory layout: -# -# Run the job. -# -# Note: It seems "LD_PRELOAD=libmpi.so" is needed to prevent the error: -# "[LOG_CAT_ML] component basesmuma is not available but requested in hierarchy". -# See https://www.mail-archive.com/users@lists.open-mpi.org/msg35048.html -mpirun -x LD_PRELOAD=libmpi.so -np $PBS_NCPUS python Optimised_APM.py +```python +from plate_model_manager import PlateModelManager +PlateModelManager().get_model("", data_dir="data/") ``` -## Configuration - -All settings are now in "Optimised_config.py", such as the model name, start/end times, number of models to use, etc. So the first step is to edit that to ensure it is configured how you like. The most important parameter is `data_model` which refers to the plate model data (mentioned above). - -### Optimisation overview +When listing the model's rotation and topology files in `Optimised_config.py`, point at wherever the +files actually live (PMM, for example, stores them in `Rotations/` and `Topologies/` sub-directories). -Essentially the workflow optimises the absolute plate motion by perturbing the absolute rotation (pole and angle) and then calculating penalties derived from that rotation (eg, net rotation (NR), trench migration (TM) and plate velocity (PV); and a hotspot (HS) penalty for 0-80Ma). Then these penalties are scaled by their weights (eg, NR=1, TM=0.5, PV=0.5) and also internal scaling to make sure each penalty is roughly the same when the weights are all 1.0 (but that’s an implementation detail). Then the penalties are added to give one number. The optimization workflow then perturbs the rotation further to try to reduce that number. It does this many times until it settles on a local minimum (called a seed), and does this for a bunch of seeds to find global minimum. +## Configuration -### Seed screening ("screen then polish") +All settings live in [`Optimised_config.py`](Optimised_config.py). The most important is `data_model` +(the model sub-directory in `data/`). The config also sets the model name (suffixed onto output +filenames), the start/end ages and interval, the number of seeds, the per-component weights/bounds, +and the reference-plate / seed scheme. Most other knobs have sensible defaults. -By default the workflow now *screens* all seeds before optimising. The objective function is evaluated once at every seed (cheap - roughly the cost of fully optimising just one or two seeds, since a full NLopt optimisation typically uses ~80-200 objective evaluations), and then only the most promising seeds are fully optimised: the `seed_screen_top_n` best-screened seeds plus `seed_screen_uniform_n` seeds spread uniformly across the search space (insurance in case screening misses a basin). Both parameters are in "Optimised_config.py". +### Seed screening ("screen then polish") - the ~10x speedup -An empirical convergence study on a 1Ga global plate model found that screening 400 seeds and fully optimising the best 16 (plus 16 uniform) reproduces the full 400-seed multistart minimum to within ~0.5% cost, while reducing the per-timestep optimisation cost by roughly a factor of 10. The study also found that the optimised pole is only loosely constrained near the minimum (several near-degenerate minima can differ by 10-20 degrees in pole position while differing by less than 0.5% in cost), so the residual error introduced by screening is well below the intrinsic uncertainty of the objective function itself. +The objective function is first evaluated once at **every** seed (cheap screening - roughly the cost +of fully optimising one or two seeds, since a full COBYLA run uses ~80-200 evaluations). Only the most +promising seeds are then fully optimised: the `seed_screen_top_n` best-screened seeds plus +`seed_screen_uniform_n` seeds spread uniformly across the search space (insurance in case screening +misses a basin). -Summary of the convergence study (two representative timesteps; full results in "model_output/seed_study/"): +An empirical convergence study on a 1 Ga global plate model (~1,600 full COBYLA optimisations at two +representative timesteps) found that screening 400 seeds and fully optimising the best 16 (+16 uniform) +reproduces the full 400-seed multistart minimum to within ~0.5% cost at roughly 10x less CPU. Near the +minimum the cost surface is close to degenerate (poles 10-20° apart can differ by <0.5% in cost), so +the residual error from screening sits below the intrinsic resolving power of the objective function. -| Strategy | Full optimisations per timestep | Best cost vs 400-seed multistart | Approx. CPU cost | +| Strategy | Full optimisations per timestep | Best cost vs 400-seed multistart | Approx. CPU | |---|---|---|---| -| 400-seed multistart (original default) | 400 | - (reference) | 1.0x | +| 400-seed multistart (original) | 400 | - (reference) | 1.0x | | 100-seed uniform grid | 100 | -0.2% / -0.3% | 0.25x | | 49-seed uniform grid | 49 | +1.5% / -0.3% | 0.12x | -| screen all 400, polish best 16 + 16 uniform | ~30 | +0.5% / +0.3% | ~0.1x | - -(Negative values mean the reduced strategy found a slightly *better* minimum than the full multistart - the -multistart minimum itself is reproducible only to within ~0.5% because the optimiser converges to one of -several near-degenerate minima.) +| screen 400, polish best 16 + 16 uniform | ~30 | +0.5% / +0.3% | ~0.1x | -To reproduce or extend the study (eg, for a different plate model or different seed counts), use the standalone, resumable harness "seed_study.py": +Set `seed_screen_top_n = None` to disable screening and fully optimise every seed (original behaviour). +The reproducible study harness ships as [`seed_study.py`](seed_study.py): ``` - python seed_study.py prep --age 100 # prepare data for one timestep (re-run until 'PREP DONE') - python seed_study.py probe --age 100 # time a single objective evaluation / full optimisation - python seed_study.py run --age 100 # run the study (resumable; re-run until 'ALL TASKS DONE') - python seed_study.py report --age 100 # summarise results (writes a CSV) +python seed_study.py prep --age 100 # prepare one timestep (re-run until 'PREP DONE') +python seed_study.py probe --age 100 # time one objective evaluation / full optimisation +python seed_study.py run --age 100 # run the study (resumable; re-run until 'ALL TASKS DONE') +python seed_study.py report --age 100 # summarise results (writes a CSV) ``` -To disable screening (and fully optimise every seed, as in the original workflow), set `seed_screen_top_n = None`. - -The workflow also caps each NLopt (COBYLA) optimisation at `nlopt_max_eval_safety` objective evaluations (default 1000). COBYLA occasionally fails to converge and can otherwise cycle for thousands of evaluations, stalling an entire MPI rank (and hence the whole timestep, since all ranks must finish before the next timestep begins). +A safety cap (`nlopt_max_eval_safety`, default 1000) stops COBYLA from cycling for thousands of +evaluations on the rare timestep where it fails to converge (which would otherwise stall a whole MPI +rank, and hence the whole timestep). ### Trench migration scheme -The original trench migration component drives trench-normal migration velocities toward zero. However, since zero net rotation also implies near-zero trench migration, that constraint is largely redundant with the net rotation minimisation. It is also at odds with observations: most trenches roll back toward the ocean basin behind them rather than advancing toward the overriding plate - across reference frames, 62-78% of trench segments retreat, with mean trench-normal velocities of +1.3-1.5 cm/yr and medians of +0.9-1.3 cm/yr (Schellart et al. 2008, Earth-Science Reviews; see also Williams et al. 2015, EPSL). +The original TM component drove trench-normal velocities toward zero, which is largely redundant with +the net-rotation minimisation (zero NR implies near-zero TM) and contradicts observations: across +reference frames 62-78% of trench segments retreat, with mean trench-normal velocities of +1.3-1.5 +cm/yr (Schellart et al. 2008; Williams et al. 2015). -The default scheme is therefore now `trench_migration_scheme = 'rollback'` (in "Optimised_config.py"): per-trench orthogonal velocities are driven toward a target of +10 mm/yr retreat, and the bounds on the *mean* orthogonal velocity are tightened to (0, 20) mm/yr - the global mean must be retreating, but at less than 2 cm/yr. Individual trenches can still advance; only the mean is required to retreat. This makes trench kinematics an independent constraint that competes with the net rotation minimisation (instead of duplicating it). Set `trench_migration_scheme = 'minimise'` (or `OPTAPM_TM_SCHEME=minimise`) to restore the original behaviour. +The default is therefore now `trench_migration_scheme = 'rollback'`: per-trench orthogonal velocities +are driven toward a target of **+10 mm/yr retreat**, and the bounds on the *mean* orthogonal velocity +are tightened to (0, 20) mm/yr - the global mean must retreat, but at less than 2 cm/yr. Individual +trenches may still advance. This makes trench kinematics an independent constraint that genuinely +competes with the NR minimisation instead of duplicating it. -A single-timestep validation (5 Ma) of the two schemes, showing that the rollback scheme makes trench kinematics an independent constraint (net rotation rises instead of being co-minimised, and the trench population becomes retreat-dominated as observed): - -| Scheme | Mean trench-normal velocity | Segments retreating | Net rotation at optimum | +| Scheme | Mean trench-normal velocity | Segments retreating | NR at optimum | |---|---|---|---| | `minimise` (original) | -1.8 mm/yr (net advance) | 56% | 0.144 deg/Myr | | `rollback` (new default) | +1.8 mm/yr (net retreat) | 63% | 0.175 deg/Myr | -The parameter choices are: target retreat of +10 mm/yr = the centre of the observed median trench-normal retreat of +0.9 to +1.3 cm/yr across reference frames (Schellart et al. 2008); mean bounds of (0, 20) mm/yr = the global mean must retreat (62-78% of segments retreat in all reference frames examined by Schellart et al. 2008) but more slowly than ~2 cm/yr (just above the observed mean retreat of +1.3-1.5 cm/yr). - -### Uncertainty quantification (run variants) +Set `trench_migration_scheme = 'minimise'` (or `OPTAPM_TM_SCHEME=minimise`) to restore the original +behaviour. -The optimisation is dominated by the net rotation (NR) minimisation: subduction zone migration largely depends on the net rotation optimisation (it is not an independent parameter), and the same holds for limiting the speed of continents (see [Muller et al. 2022, Solid Earth](https://doi.org/10.5194/se-13-1127-2022)). Perturbing the component weights only slightly changes the outcome, so a defensible uncertainty envelope for the optimised reference frame instead consists of end-members in the net rotation bounds: +### Run variants (uncertainty quantification) -1. **No-net-rotation (NNR) end-member** - zero net rotation. This rotation file is already produced by every run as "no_net_rotation_model_.rot". -2. **Best run** (`run_variant = 'best'`, the default) - NR bounded to (0.08, 0.20) deg/Myr: non-zero (as suggested by mantle flow models, eg, Becker 2006) but below the preferred geodynamic upper limit (Conrad & Behn 2010). -3. **Maximum-NR end-member** (`run_variant = 'nr_max'`) - the NR upper bound relaxed to 0.30 deg/Myr *and the NR weight halved* (inverse weight 2.0). The 0.30 deg/Myr ceiling is the top of the range permitted by asthenospheric shear / seismic anisotropy constraints (Conrad & Behn 2010 derive NR < 0.26 deg/Myr for an asthenosphere 10x less viscous than the upper mantle; anisotropy proxies allow 0.2-0.3 deg/Myr); larger values approach the Pacific hotspot (HS3) frame rates (0.33-0.44 deg/Myr), which are widely considered geodynamically implausible. The weight is halved because relaxing the bounds alone is not sufficient to produce a directional end-member: the bounds are hard penalty walls that only affect timesteps where the optimum presses against the 0.20 deg/Myr ceiling (mostly in deep time), whereas inside the walls the solution is set by the smooth trade-off between the NR cost and the other components. Halving the NR weight shifts that trade-off at every timestep. Single-timestep tests (5 Ma) confirm this: relaxing the bounds alone merely selected a different near-degenerate minimum (NR 0.123 deg/Myr), whereas bounds + halved weight produced the intended systematically stronger-rollback, higher-NR solution (NR 0.177 deg/Myr, trench retreat strengthening from +1.8 to +3.2 mm/yr, 68% of segments retreating). - -To produce the envelope, run the workflow twice (the `nr_max` variant appends its name to the model name, so the two runs write separate output files): - -``` - mpirun -np python Optimised_APM.py - OPTAPM_VARIANT=nr_max mpirun -np python Optimised_APM.py -``` +The optimisation is dominated by the NR minimisation - trench migration and continent speeds are not +independent of it (Müller et al. 2022). Perturbing the component weights barely changes the outcome, so +the uncertainty envelope is expressed as **end-members in the NR bounds** instead: -The variant NR bounds are defined in `run_variant_nr_bounds` in "Optimised_config.py" (where further variants can be added). +1. **No-net-rotation (NNR)** - the zero-NR end-member. Produced by *every* run as + `no_net_rotation_model_.rot`. +2. **Best** (`run_variant = 'best'`, the default) - NR bounded to (0.08, 0.20) deg/Myr: non-zero (as + suggested by mantle-flow models, e.g. Becker 2006) but below the preferred geodynamic upper limit + (Conrad & Behn 2010). +3. **Maximum-NR** (`run_variant = 'nr_max'`) - NR upper bound relaxed to 0.30 deg/Myr **and the NR + weight halved** (inverse weight 2.0). Relaxing the bounds alone is not enough to produce a + directional end-member (the bounds are hard walls that only bite where the optimum presses against + the ceiling); halving the weight shifts the smooth NR-vs-other trade-off at *every* timestep, + yielding the intended systematically stronger-rollback / higher-NR solution. -### Model diagnostics (summary statistics and plots) +The canonical objective-function parameters of the two runs (as printed in the run logs) are: -At the end of every run (if `generate_diagnostics = True`, the default), the workflow computes per-timestep summary statistics of the optimised model and writes them to "model_output/": +| Component | `best` run | `nr_max` run | +|---|---|---| +| Net rotation | weight 1.0, bounds 0.08-0.20 deg/Myr | **weight 2.0 (halved)**, bounds 0.08-**0.30** deg/Myr | +| Trench migration (rollback) | weight 1.0, mean bounds 0-20 mm/yr, target +10 mm/yr | weight 1.0, mean bounds 0-20 mm/yr, target +10 mm/yr | +| Plate velocity | weight 1.0, bounds 0-60 mm/yr | weight 1.0, bounds 0-60 mm/yr | +| Hotspot trails (0-80 Ma only) | weight 1.0 | weight 1.0 | -* `_diagnostics.csv` - per timestep: median and MAD (median absolute deviation) of the net rotation rate (deg/Myr, sampled at 1 Myr sub-steps within each timestep), and mean, median and MAD of the trench-normal migration velocities (mm/yr, positive = retreat) across all subduction zone segments, plus the percentage of segments retreating. -* `_net_rotation.png` - net rotation (median +/- MAD) through time, with the 0.20 deg/Myr preferred upper bound and the 0.26 deg/Myr Conrad & Behn (2010) limit marked. -* `_trench_migration.png` - trench-normal migration (mean +/- MAD) through time, with the observed present-day mean retreat range (+0.9 to +1.5 cm/yr, Schellart et al. 2008) marked. +Reference plate 701 (Africa); seed = auto-calculated paleomagnetic pole; uniform sampling; search +radius 180°. The two runs differ **only** in the net-rotation bounds and weight; everything else is +identical. -The same statistics and plots are also generated for the no-net-rotation model (outputs named `_NNR_*`). This shows how subduction zone migration behaves in an NNR world (the zero-NR end-member of the uncertainty envelope), and the NNR net rotation plot doubles as a sanity check (it should be ~zero at all times). +> **Weight convention.** Weights are *inverse* weights - each component cost is multiplied by +> `1.0 / weight` (so a higher weight means a *weaker* constraint). The values above are the 0-80 Ma +> values; for ages > 80 Ma the trench-migration and plate-velocity inverse weights become 2.0 (i.e. a +> multiplicative weight of 0.5), and the hotspot-trail constraint applies only over 0-80 Ma. The +> net-rotation weight is constant over all ages (1.0 for `best`, 2.0 for `nr_max`). -The diagnostics can be (re)generated standalone on an existing optimised model without re-running the optimisation: +Run the workflow twice to produce the envelope (the `nr_max` variant appends its name to the model name, +so the runs write separate output files): ``` - python model_diagnostics.py +python Optimised_APM.py # 'best' -> optimised + NNR +OPTAPM_VARIANT=nr_max python Optimised_APM.py # 'nr_max' -> maximum-NR optimised ``` -(This uses the current settings in "Optimised_config.py", including any `OPTAPM_*` environment overrides, eg, `OPTAPM_VARIANT=nr_max python model_diagnostics.py` to diagnose the nr_max run.) +The variant NR bounds and weights live in `run_variant_nr_bounds` / `run_variant_nr_weight` in the +config; further variants are a one-line addition. -### Quick test runs - -The age range, seed count, screening parameters and parallelisation can be temporarily overridden via environment variables without editing "Optimised_config.py" - useful for quick smoke tests: - -``` - OPTAPM_START_AGE=5 OPTAPM_END_AGE=0 OPTAPM_MODELS=49 OPTAPM_SERIAL=1 python Optimised_APM.py -``` +### Diagnostics -(`OPTAPM_SCREEN_TOP_N` and `OPTAPM_SCREEN_UNIFORM_N` override the screening parameters; a negative `OPTAPM_SCREEN_TOP_N` disables screening.) +If `generate_diagnostics = True` (the default), every run finishes by writing per-timestep summary +statistics and plots to `model_output/`: -### Plate velocity penalty +* `_diagnostics.csv` - net rotation (median ± MAD) and trench-normal migration (mean, + median ± MAD, % retreating) per timestep, +* `_net_rotation.png` - NR through time, with the 0.20 and 0.26 deg/Myr geodynamic limits + marked, +* `_trench_migration.png` - trench-normal migration through time, with the observed + present-day retreat range marked. -For the 1Ga model, the plate velocity (PV) penalty currently involves multiplying the plate velocity weight by the median velocity across all continents on the globe: +The same outputs are produced for the no-net-rotation model (`_NNR_*`), whose NR plot should +be ~zero (a built-in sanity check). Regenerate standalone without re-running the optimisation: ``` -plate_velocity_penalty = plate_velocity_weight * median_velocity_across_all_continents +python model_diagnostics.py ``` -Previously we performed experiments involving continent contouring where we also multiplied by the global fragmentation index (sum of continent areas divided sum of continent perimeters): - -``` -plate_velocity_penalty = plate_velocity_weight * median_velocity_across_all_continents * - sum(area_continent) / sum(perimeter_continent) -``` +### Quick test runs -Other experiments involved individually penalizing continents, such as: +Override the age range, seed count and parallelisation via environment variables (no config edit needed): ``` -plate_velocity_penalty = plate_velocity_weight * - mean(median_velocity_in_continent * area_continent / perimeter_continent) +OPTAPM_START_AGE=5 OPTAPM_END_AGE=0 OPTAPM_MODELS=49 OPTAPM_SERIAL=1 python Optimised_APM.py ``` -The general idea was to penalize the plate velocities of supercontinents more heavily in order to slow them down. However this tended to produce less than desirable results including conflicting with optimizing net rotation, and so was ultimately abandoned. However the code to perform contouring of continents remains and a brief overview is provided in the following section (since it may be useful in other projects). +`OPTAPM_SCREEN_TOP_N` / `OPTAPM_SCREEN_UNIFORM_N` override the screening counts (a negative +`OPTAPM_SCREEN_TOP_N` disables screening); `OPTAPM_VARIANT` selects the run variant; +`OPTAPM_TM_SCHEME` selects the trench scheme. -## Running the optimisation workflow +## Running the optimisation -The optimisation workflow can be run in serial or parallel. In parallel it can be run using `ipyparallel` or `mpi4py`. +Each run produces the optimised rotation file +`data//optimisation/optimised_rotation_model_.rot` (a single file containing the +entire optimised rotation model) plus the no-net-rotation file +`no_net_rotation_model_.rot`. -Each of these should produce the final optimised rotation file "optimised_rotation_model_.rot", -in the "data//optimisation/" directory, where ** and ** are defined -by the `data_model` and `model_name` variables in the "Optimised_config.py" script. +**Serial** (set `use_parallel = None` in the config, or `OPTAPM_SERIAL=1`): -### To run in serial +``` +python Optimised_APM.py +``` -Edit "Optimised_config.py" and change the `use_parallel` parameter to `None` and then run: +**Parallel with MPI** (set `use_parallel = MPI4PY`) - the easy way to use all the cores on a +laptop/desktop, and the way to scale across an HPC: ``` - python Optimised_APM.py +mpiexec -n python Optimised_APM.py ``` -### To run in parallel using `ipyparallel` +On an HPC, edit and submit the PBS script instead: `qsub Optimised_APM.sh`. -This is useful when running a Jupyter notebook since it supports `ipyparallel` by either: +**Parallel with ipyparallel** (set `use_parallel = IPYPARALLEL`) - useful inside a Jupyter notebook; +start the cluster from the directory containing `Optimised_APM.py` (`ipcluster start -n `) to +avoid an `ImportError: No module named objective_function`. -* Starting clusters in the Jupyter clusters tab, or -* Running `ipcluster start -n ` to start engines on *cores* number of cores manually (eg, if not using a notebook). - * **NOTE**: You should be in the directory containing "Optimised_APM.py" when you start the cluster - to avoid the error `ImportError: No module named objective_function`. +> Copy *all* directories to the HPC, even empty ones like `model_output/`, or an exception will be +> raised during execution and `mpirun` may terminate without a clear message. -Edit "Optimised_config.py" and change the `use_parallel` parameter to `IPYPARALLEL` and then run: +## Assembling all reference frames -``` - python Optimised_APM.py -``` +Once **both** optimisation runs are finished, [`add_ref_frames_to_unoptimised.py`](add_ref_frames_to_unoptimised.py) +assembles a single, comprehensive rotation model that contains every available absolute reference frame +of the plate model as a distinct (high-numbered) plate ID: -### To run in parallel using `mpi4py` +| Plate ID | Reference frame | +|---|---| +| 14 | paleomag (the model's native frame) | +| 15 | optimised mantle (the `best` run) | +| 16 | optimised mantle, maximum net rotation (the `nr_max` run) | +| 17 | no-net rotation | +| 18 | approximate true polar wander (TPW = paleomag − optimised) | -This is useful when running on a High Performance Computing (HPC) cluster since MPI is used to -spread the parallel workload across any number of nodes/cores. +> Plate ID 16 (the maximum-NR optimised frame) is new - earlier versions produced only a single +> optimised frame, so the NNR and TPW frames have shifted from 16/17 to **17/18**. -Edit "Optimised_config.py" and change the `use_parallel` parameter to `MPI4PY`. - -If you are running on a personal computer that has an MPI runtime installed then run: +The script reads the original un-optimised rotation files plus the three 005-000 reference-frame files +produced by the runs (best optimised, `nr_max` optimised, no-net rotation), and writes new rotation +files (suffix `__with_ref_frames`) into the `optimisation/` directory. Input and output +paths are taken from `Optimised_config.py`, so it normally needs no editing: ``` - mpiexec -n python Optimised_APM.py +python add_ref_frames_to_unoptimised.py ``` -...where *cores* is the number of cores to use. - -If you are running on a HPC cluster (such as NCI's Gadi) then edit the "Optimise_APM.sh" job submission script with the number of cpus, amount of memory and walltime, etc, and then run: +The resulting hierarchy is: ``` - qsub Optimise_APM.sh + 015 016 017 018 + \ \ | / + \---------/ + | + 000 <- the "default" reference frame (default_reference_frame_plate_id) + | + 014 <- paleomag (the un-optimised model hangs off here) + | + --------- + | | + 101 701 ... ``` -...to submit to the job queue. When the job is finished you should have the final optimised rotation file mentioned above. +Anchor (in GPlates / pyGPlates) to any of 014-018 to view the model in that reference frame; anchoring +000 uses whichever frame `default_reference_frame_plate_id` is set to (paleomag by default). The +reference-frame sequences are bundled into the first output rotation file (not a separate file) so that +plate 000 is always defined wherever the model is loaded. -> Note: The "Optimise_APM.sh" job submission script was used on NCI's Gadi HPC and may require modifications for other HPC systems. +> **The script requires the un-optimised input model to be in its paleomagnetic frame** (anchoring +> plate 000 of the input gives PMAG) - plate 14 is set to that native input frame and the optimised / +> NNR frames are expressed relative to it. This is exposed as the clearly-flagged +> `unoptimised_input_is_in_pmag_frame` switch at the top of the script's input parameters; the script +> **aborts** if it is left `False` rather than silently mislabelling the frames. If a model's native +> 000 frame is not paleomagnetic (e.g. the Zahirovic et al. 2022 model, whose PMAG frame is reached via +> anchor plate 701701), re-express the input so that anchor 000 = PMAG first. -> Note: Make sure to copy all directories over to the HPC (even empty directories like "model_output") otherwise an exception -will get raised during execution and mpirun (or mpiexec) will get terminated abruptly (possibly without an error message). +## Results and post-processing -## Results +The optimised rotation file is the entire optimised rotation model - it is the only rotation file you +need to load into GPlates. The optimised rotations are also written back into copies of the original +files (in `optimisation/`); only rotations that referenced plate 000 are modified. -After running the workflow, the optimised rotations will be in a file with a name like "optimisation/optimised_rotation_model_run1.rot" -(depending on the model name in "Optimised_config.py"). Note that this one rotation file contains the entire optimised rotation model -(in other words, it includes all rotations and hence is the only rotation file you need to load into GPlates). - -The optimised rotations will also be saved back into the original files (or copies of them in "optimisation/"). -All these rotation files also comprise the entire optimised rotation model. -Note that only those rotations that referenced 000 as their fixed plate will be modified (to include the optimised absolute plate motion). - -You can also optionally remove plates in the plate hierarchy to make it simpler (eg, plates below 701 are typically removed so that 701 directly references 000). -This can be done using the 'remove_plate_rotations' module in GPlately ( https://github.com/GPlates/gplately ) For example: +You can optionally simplify the plate hierarchy (e.g. remove plates below 701 so that 701 references 000 +directly) with GPlately's `remove_plate_rotations`: ``` - python -m gplately.ptt.remove_plate_rotations -u -a 0.1 5 -p 70 4 3 2 1 -o removed_70_4_3_2_1_ -- ... +python -m gplately.ptt.remove_plate_rotations -u -a 0.1 5 -p 70 4 3 2 1 -o removed_70_4_3_2_1_ -- ``` -...where you replace `...` with all the rotation files in the optimised rotation model. - -## Plotting results - -After running the workflow and post-processing (although you don't need to run `gplately.ptt.remove_plate_rotations` for this), you can plot the -trench migration stats and net rotation curves for the non-optimised and any optimised models using the Jupyter notebooks in the `figures/` directory. +## Plotting +The Jupyter notebooks in [`figures/`](figures/) plot trench-migration statistics and net-rotation +curves for the un-optimised and any optimised models (you do not need to run `remove_plate_rotations` +first). ## References for parameter choices -* Atkins, S., & Coltice, N. (2021). Constraining the range and variation of lithospheric net rotation using geodynamic modeling. *Journal of Geophysical Research: Solid Earth*, 126, e2021JB022057. https://doi.org/10.1029/2021JB022057 - convection models produce Earth-like net rotation mostly well below 0.2 deg/Myr. -* Becker, T. W. (2006). On the effect of temperature and strain-rate dependent viscosity on global mantle flow, net rotation, and plate-driving forces. *Geophysical Journal International*, 167, 943-957. https://doi.org/10.1111/j.1365-246X.2006.03172.x - mantle flow models suggest net rotation should be positive and non-zero (basis for the 0.08 deg/Myr lower bound). -* Conrad, C. P., & Behn, M. D. (2010). Constraints on lithosphere net rotation and asthenospheric viscosity from global mantle flow models and seismic anisotropy. *Geochemistry, Geophysics, Geosystems*, 11, Q05W05. https://doi.org/10.1029/2009GC002970 - net rotation < 0.26 deg/Myr for an asthenosphere 10x less viscous than the upper mantle (basis for the 0.20 deg/Myr preferred and 0.30 deg/Myr maximum upper bounds). -* Muller, R. D., et al. (2022). A tectonic-rules-based mantle reference frame since 1 billion years ago - implications for supercontinent cycles and plate-mantle system evolution. *Solid Earth*, 13, 1127-1159. https://doi.org/10.5194/se-13-1127-2022 - shows that net rotation minimisation dominates the optimisation and that trench migration is not an independent constraint under the original ('minimise') scheme (basis for the uncertainty-envelope design and the 'rollback' scheme). -* Schellart, W. P., Stegman, D. R., & Freeman, J. (2008). Global trench migration velocities and slab migration induced upper mantle volume fluxes: Constraints to find an Earth reference frame based on minimizing viscous dissipation. *Earth-Science Reviews*, 88, 118-144. https://doi.org/10.1016/j.earscirev.2008.01.005 - 62-78% of trench segments retreat across reference frames, with mean trench-normal retreat of +1.3-1.5 cm/yr and medians of +0.9-1.3 cm/yr (basis for the +10 mm/yr rollback target and the (0, 20) mm/yr mean bounds). -* Tetley, M. G., Williams, S. E., Gurnis, M., Flament, N., & Muller, R. D. (2019). Constraining absolute plate motions since the Triassic. *Journal of Geophysical Research: Solid Earth*, 124, 7231-7258. https://doi.org/10.1029/2019JB017442 - the original optAPM optimisation framework. -* Williams, S., Flament, N., Muller, R. D., & Butterworth, N. (2015). Absolute plate motions since 130 Ma constrained by subduction zone kinematics. *Earth and Planetary Science Letters*, 418, 66-77. https://doi.org/10.1016/j.epsl.2015.02.026 - trench kinematic plausibility (dominantly slow retreat, minimal advance) as a primary constraint on absolute plate motion. +* Atkins, S., & Coltice, N. (2021). Constraining the range and variation of lithospheric net rotation + using geodynamic modeling. *JGR Solid Earth*, 126, e2021JB022057. + https://doi.org/10.1029/2021JB022057 +* Becker, T. W. (2006). On the effect of temperature- and strain-rate-dependent viscosity on global + mantle flow, net rotation, and plate-driving forces. *GJI*, 167, 943-957. + https://doi.org/10.1111/j.1365-246X.2006.03172.x - basis for the 0.08 deg/Myr NR lower bound. +* Conrad, C. P., & Behn, M. D. (2010). Constraints on lithosphere net rotation and asthenospheric + viscosity from global mantle flow models and seismic anisotropy. *G-cubed*, 11, Q05W05. + https://doi.org/10.1029/2009GC002970 - basis for the 0.20 (preferred) and 0.30 (maximum) deg/Myr NR + upper bounds. +* Müller, R. D., et al. (2022). A tectonic-rules-based mantle reference frame since 1 billion years ago. + *Solid Earth*, 13, 1127-1159. https://doi.org/10.5194/se-13-1127-2022 - NR minimisation dominates the + optimisation; basis for the uncertainty-envelope design and the rollback scheme. +* Schellart, W. P., Stegman, D. R., & Freeman, J. (2008). Global trench migration velocities and slab + migration induced upper mantle volume fluxes. *Earth-Science Reviews*, 88, 118-144. + https://doi.org/10.1016/j.earscirev.2008.01.005 - basis for the +10 mm/yr rollback target and the + (0, 20) mm/yr mean bounds. +* Tetley, M. G., Williams, S. E., Gurnis, M., Flament, N., & Müller, R. D. (2019). Constraining absolute + plate motions since the Triassic. *JGR Solid Earth*, 124, 7231-7258. + https://doi.org/10.1029/2019JB017442 - the original optAPM framework. +* Williams, S., Flament, N., Müller, R. D., & Butterworth, N. (2015). Absolute plate motions since 130 + Ma constrained by subduction zone kinematics. *EPSL*, 418, 66-77. + https://doi.org/10.1016/j.epsl.2015.02.026 - trench kinematic plausibility as an APM constraint. diff --git a/add_ref_frames_to_unoptimised.py b/add_ref_frames_to_unoptimised.py index 20e507f..c7e40ae 100644 --- a/add_ref_frames_to_unoptimised.py +++ b/add_ref_frames_to_unoptimised.py @@ -1,103 +1,183 @@ -import glob +""" +add_ref_frames_to_unoptimised.py +================================ + +Post-processing step of the optAPM workflow. + +After the optimisation has been run, this script reads the original un-optimised rotation +files together with the 005-000 reference-frame rotations produced by the workflow, and +writes a new set of rotation files that contain ALL the available absolute reference frames +of a plate model as distinct (high-numbered) plate IDs: + + 14 pmag - the model's native paleomagnetic reference frame + 15 optimised - the 'best' optimised mantle reference frame + 16 optimised_max_NNR - the 'nr_max' optimised mantle reference frame + (the maximum net-rotation end-member of the uncertainty envelope) + 17 no_net_rotation - the no-net-rotation reference frame (the zero net-rotation end-member) + 18 true_polar_wander - an approximation of true polar wander, TPW = pmag - optimised(best) + +This is the step that "assembles a comprehensive rotation file with all the available +reference frames for a given plate model". Each optimisation run produces its own optimised +005-000 reference frame, so this script requires both optimisation runs (the default 'best' +run and the 'nr_max' run) plus the no-net-rotation file that every run also produces: + + 1) python Optimised_APM.py # 'best' run -> optimised + no-net-rotation + 2) OPTAPM_VARIANT=nr_max python Optimised_APM.py # 'nr_max' run -> optimised_max_NNR + 3) python add_ref_frames_to_unoptimised.py # this script -> comprehensive rotation file(s) + +The resulting rotation hierarchy looks like: + + 015 016 017 018 + \ \ | / + \------------/ + | + 000 <- the "default" reference frame (one of 014-018) + | + 014 <- pmag (the un-optimised model is attached here) + | + --------- + | | + 101 701 ... + +Anchoring (in GPlates/pyGPlates) to any of 014/015/016/017/018 reconstructs the model in +that reference frame. Anchoring 000 uses whichever frame `default_reference_frame_plate_id` +is set to below. + +IMPORTANT - reference frame of the un-optimised input model. This script requires the +un-optimised input model to be in its paleomagnetic frame (anchoring plate 000 of the input +gives PMAG), because plate 014 is set to that native input frame and the optimised / no-net +rotation 005-000 frames are all expressed relative to it. This requirement is exposed as the +clearly-flagged `unoptimised_input_is_in_pmag_frame` switch in the "Input parameters" section +below; the script aborts if it is not set True. If a model's native 000 frame is NOT +paleomagnetic (for example the Zahirovic et al. 2022 model, whose native frame is a mantle +frame and whose PMAG frame is reached by anchoring plate 701701), the input must first be +re-expressed so that anchor 000 = PMAG before running this script. + +By default all input/output filenames are derived from "Optimised_config.py" (the same config +the optimisation used), so this script normally needs no editing. Every parameter can still be +overridden in the "Input parameters" section below. +""" + import os.path import pygplates - -# -# This script reads the original un-optimised rotation files and the optimised and -# the no-net rotation reference frames, and generates a new set of rotation files -# with various reference frames added. -# -# Currently the added reference frames are the original paleomag (PMAG), the optimised mantle, -# the no-net rotation, and an approximation of true polar wander (TPW) calculated by replacing -# the PMAG reference frame with the optimised mantle reference frame. -# +import Optimised_config as config #### Input parameters #### -# Original (un-optimised) rotation files. -input_rotation_filenames = [ - '1000_0_rotfile.rot', - '1800_1000_rotfile.rot', -] -#input_rotation_filenames = glob.glob('*.rot') - -# This file contains the 005-000 optimised reference frame. -input_optimised_rotation_file = 'optimisation/optimised_rotation_model_20240725.rot' - -# This file contains the 005-000 no-net-rotation reference frame. -input_no_net_rotation_file = 'optimisation/no_net_rotation_model_20240725.rot' - -# The suffix to append to the basenames of the input rotation filenames to get the output rotation filenames. +# ===================================================================================== +# IMPORTANT - reference frame of the un-optimised INPUT model (read this!) +# ===================================================================================== +# The un-optimised input rotation files MUST be in their paleomagnetic (PMAG) frame: +# anchoring plate 000 of the input must give PMAG. Plate (014) +# in the output is set to this native input frame, and ALL of the optimised / no-net +# rotation 005-000 reference frames are expressed relative to it. The whole script (and +# the optimisation that produced the 005-000 files) depends on this assumption. # -# Note: The output rotation files also end up in the 'optimisation/' sub-directory. -output_rotation_filenames_suffix = '_20240725' - -# Plate IDs of reference frames. -# -# These are PMAG, optimised, no-net rotation and true polar wander (which is optimised with PMAG removed). +# Set this True ONLY once you have confirmed the input is in its PMAG frame. +# +# If your model's native 000 frame is NOT paleomagnetic - e.g. the Zahirovic et al. (2022) +# model, whose native frame is a mantle frame and whose PMAG frame is reached by anchoring +# plate 701701 - you must FIRST re-express the input rotations so that anchor 000 = PMAG, +# otherwise plate 014 will not be the paleomagnetic frame. (The script aborts below if this +# is left False, rather than silently mislabelling the reference frames.) +# ===================================================================================== +unoptimised_input_is_in_pmag_frame = True + +# Directory containing the data model (the 'data//' directory). All input and +# output paths below are resolved relative to this directory, so the script can be run from +# anywhere (it does not need to be run from inside the data model directory). +data_model_dir = os.path.join(config.datadir, config.data_model) + +# The 'best' (default variant) model name and the 'nr_max' variant model name. +# +# "Optimised_config.py" appends the run variant to the model name for any variant other than +# 'best', so reconstruct the base ('best') name robustly regardless of any OPTAPM_VARIANT in +# the current environment. +_base_model_name = config.model_name +if config.run_variant != 'best' and _base_model_name.endswith('_' + config.run_variant): + _base_model_name = _base_model_name[:-len('_' + config.run_variant)] +best_model_name = _base_model_name +nr_max_model_name = _base_model_name + '_nr_max' + +# Original (un-optimised) rotation files (absolute paths). +# +# Taken from the optimisation config so the same input files are used. These are the files +# whose fixed plate 000 (the PMAG frame) becomes plate 014 in the output. +input_rotation_filenames = [os.path.join(config.datadir, rotation_filename) + for rotation_filename in config.original_rotation_filenames] + +# The optimised / no-net-rotation 005-000 reference-frame files produced by the two runs. +input_optimised_rotation_file = os.path.join( + data_model_dir, 'optimisation', 'optimised_rotation_model_' + best_model_name + '.rot') +input_optimised_max_NNR_rotation_file = os.path.join( + data_model_dir, 'optimisation', 'optimised_rotation_model_' + nr_max_model_name + '.rot') +input_no_net_rotation_file = os.path.join( + data_model_dir, 'optimisation', 'no_net_rotation_model_' + best_model_name + '.rot') + +# The suffix appended to the basename of each input rotation file to form the output filename. +# The output files are written to the 'optimisation/' sub-directory of the data model directory. +output_rotation_filenames_suffix = '_' + best_model_name + '_with_ref_frames' + +# Plate IDs of the reference frames. pmag_ref_frame_plate_id = 14 optimised_ref_frame_plate_id = 15 -no_net_rotation_ref_frame_plate_id = 16 -true_polar_wander_ref_frame_plate_id = 17 +optimised_max_NNR_ref_frame_plate_id = 16 +no_net_rotation_ref_frame_plate_id = 17 +true_polar_wander_ref_frame_plate_id = 18 -# The default reference frame that plate 000 should correspond to. -# -# It should be one of the above reference frame plate IDs. -default_reference_frame_plate_id = optimised_ref_frame_plate_id +# The reference frame that plate 000 should correspond to (must be one of the IDs above). +default_reference_frame_plate_id = pmag_ref_frame_plate_id -# -# Note that the final rotation hierarchy will look something like... -# -# 015 016 017 -# | | | -# ----------- -# | -# 000 -# | -# 014 -# | -# ----- -# | | -# 101 701 -# | | -# -# ...where if 000 represents: -# 1. PMAG – then 014-000 is zero, and 000-015 (and 014-015) represents what 005-000 previously represented. -# 2. Optimised – then 000-015 is zero, and 014-000 (and 014-015) represents what 005-000 previously represented. -# -# That can be one version of the rotation model. -# And another can be a flattened version that removes the reference frame plate IDs (eg, 014, 015, 016, 017). -# The flattened version can be what's included with GPlately. -# +# The plate ID used to express true polar wander (TPW = pmag - optimised). TPW is the residual +# rotation of pmag not captured by the optimised mantle frame; following Tetley et al. (2019) +# it is evaluated on Africa (701). +true_polar_wander_reference_plate_id = 701 ########################## -# Each input rotation file will get a corresponding output rotation file. +# A short human-readable name for each reference frame plate ID (used in rotation descriptions). +reference_frame_names = { + pmag_ref_frame_plate_id: 'paleomag', + optimised_ref_frame_plate_id: 'optimised mantle', + optimised_max_NNR_ref_frame_plate_id: 'optimised mantle (max net-rotation)', + no_net_rotation_ref_frame_plate_id: 'no-net rotation', + true_polar_wander_ref_frame_plate_id: 'approx true polar wander', +} +if default_reference_frame_plate_id not in reference_frame_names: + raise ValueError('default_reference_frame_plate_id is not one of the reference frame plate IDs') + +# Refuse to run unless the user has confirmed the input is in its paleomagnetic frame +# (see the prominent note in the "Input parameters" section above). +if not unoptimised_input_is_in_pmag_frame: + raise RuntimeError( + "'unoptimised_input_is_in_pmag_frame' is False: this script requires the un-optimised " + "input model to be in its paleomagnetic frame (anchoring plate 000 of the input gives " + "PMAG), because plate {0} is set to that frame and the optimised/NNR frames are expressed " + "relative to it. Re-express the input so that anchor 000 = PMAG (e.g. the Zahirovic et al. " + "2022 model's PMAG frame is reached via anchor plate 701701), then set this switch True." + .format(pmag_ref_frame_plate_id)) + + +# +# Load the original (un-optimised) rotation files and relabel their PMAG frame (fixed plate 000) +# as plate (014). Each input file gets a corresponding output file. +# output_rotation_feature_collections = [pygplates.FeatureCollection(filename) for filename in input_rotation_filenames] -# The output rotation filenames associated with the input rotation files. output_rotation_filenames = [] for filename in input_rotation_filenames: - # Insert suffix into output rotation filenames. - basename, ext = os.path.splitext(filename) - output_filename = basename + output_rotation_filenames_suffix + ext - # Place in 'optimisation/' sub-directory. - output_filename = os.path.join('optimisation', output_filename) - output_rotation_filenames.append(output_filename) - -########################## - + basename, ext = os.path.splitext(os.path.basename(filename)) + output_rotation_filenames.append( + os.path.join(data_model_dir, 'optimisation', basename + output_rotation_filenames_suffix + ext)) -# Change fixed plate 000 to -# (and remove any existing leftover 005-000 rotation features). +# Change fixed plate 000 to (and drop any leftover 005-000 features). max_time_of_sequences_with_fixed_000 = 0.0 for rotation_feature_collection in output_rotation_feature_collections: - # Exclude any leftover 005-000 rotation features. def is_feature_005_000(feature): total_reconstruction_pole = feature.get_total_reconstruction_pole() if total_reconstruction_pole: @@ -105,344 +185,177 @@ def is_feature_005_000(feature): return fixed_plate_id == 0 and moving_plate_id == 5 return False rotation_feature_collection.remove(is_feature_005_000) - - # Change fixed plate 000 to . + for rotation_feature in rotation_feature_collection: total_reconstruction_pole = rotation_feature.get_total_reconstruction_pole() if total_reconstruction_pole: fixed_plate_id, moving_plate_id, rotation_sequence = total_reconstruction_pole - - # Change fixed plate ID 000 to . - # We want everything that references 000 to now reference . if fixed_plate_id == 0: - # Change the fixed plate ID from 000 to . + # Everything that referenced 000 (the PMAG frame) now references plate 014. rotation_feature.set_total_reconstruction_pole( pmag_ref_frame_plate_id, moving_plate_id, rotation_sequence) - # The oldest time in this sequence. - max_time = max(time_sample.get_time() for time_sample in rotation_sequence.get_time_samples()) - # The oldest time in all sequences (with fixed plate 000). + max_time = max(ts.get_time() for ts in rotation_sequence.get_time_samples()) max_time_of_sequences_with_fixed_000 = max(max_time_of_sequences_with_fixed_000, max_time) -# Original (un-optimised) rotation model. -# -# Note: We use the modified input rotation features (with fixed plate 000 changed to 014). -# -# Note: It has no 000 plate ID. Instead it's anchor plate is the PMAG reference frame plate ID. -original_rotation_model_anchored_014 = pygplates.RotationModel( +# The un-optimised model anchored to the PMAG frame (plate 014). Used for the TPW calculation. +original_rotation_model_anchored_pmag = pygplates.RotationModel( output_rotation_feature_collections, default_anchor_plate_id=pmag_ref_frame_plate_id) +# +# Helpers. +# def find_005_000_rotation_sequence(rotation_filename): + """Return the (moving 005, fixed 000) rotation sequence stored in an optimised/NNR file.""" rotation_sequence_005_000 = None - for rotation_feature in pygplates.FeatureCollection(rotation_filename): total_reconstruction_pole = rotation_feature.get_total_reconstruction_pole() if total_reconstruction_pole: fixed_plate_id, moving_plate_id, rotation_sequence = total_reconstruction_pole - if fixed_plate_id == 0 and moving_plate_id == 5: if rotation_sequence_005_000 is not None: - raise ValueError("Found multiple rotation sequences with moving plate 5 and fixed plate 0") + raise ValueError('Found multiple 005-000 sequences in {0}'.format(rotation_filename)) rotation_sequence_005_000 = rotation_sequence - if rotation_sequence_005_000 is None: - raise ValueError("Found zero rotation sequences with moving plate 5 and fixed plate 0") - + raise ValueError('Found no 005-000 sequence in {0}'.format(rotation_filename)) return rotation_sequence_005_000 -# Find the optimised (005-000) rotation sequence. -optimised_rotation_sequence = find_005_000_rotation_sequence(input_optimised_rotation_file) -# Find the no-net rotation (005-000) rotation sequence. +def mantle_frame_samples_relative_to_pmag(rotation_sequence_005_000): + """ + Convert an optimised/NNR 005-000 sequence into the rotation of that mantle frame relative + to PMAG, R(pmag -> frame), as a list of (time, finite_rotation) pairs. + + In an optimised output file the 005-000 feature stores R(000->005) where, by construction, + 000 is the mantle frame and 005 is PMAG. Hence R(pmag->frame) = R(005->000) = inverse[R(000->005)]. + """ + return [(ts.get_time(), ts.get_value().get_finite_rotation().get_inverse(), ts.is_enabled()) + for ts in rotation_sequence_005_000] + + +def true_polar_wander_samples_relative_to_pmag(optimised_rotation_sequence_005_000): + """ + R(pmag -> TPW) as a list of (time, finite_rotation) pairs. + + TPW is the part of the paleomagnetic frame not captured by the optimised mantle frame + (TPW = pmag - optimised), evaluated on Africa (701): + + R(pmag->TPW) = R(pmag->701) * R(000->005)_optimised + + where R(pmag->701) is the un-optimised Africa rotation in the PMAG frame and R(000->005) + is the stored optimised reference-frame rotation. + """ + samples = [] + for ts in optimised_rotation_sequence_005_000: + time = ts.get_time() + optimised_finite_rotation = ts.get_value().get_finite_rotation() # R(000->005) + africa_in_pmag = original_rotation_model_anchored_pmag.get_rotation( # R(pmag->701) + time, true_polar_wander_reference_plate_id) + samples.append((time, africa_in_pmag * optimised_finite_rotation, ts.is_enabled())) + return samples + + +def identity_samples(begin_time, end_time=0.0): + """Two identity samples spanning [end_time, begin_time].""" + return [(end_time, pygplates.FiniteRotation(), True), + (begin_time, pygplates.FiniteRotation(), True)] + + +# Read the three optimised/NNR 005-000 reference-frame sequences. +optimised_sequence = find_005_000_rotation_sequence(input_optimised_rotation_file) +optimised_max_NNR_sequence = find_005_000_rotation_sequence(input_optimised_max_NNR_rotation_file) no_net_rotation_sequence = find_005_000_rotation_sequence(input_no_net_rotation_file) +# Each reference frame expressed RELATIVE TO PMAG (014): R(pmag -> frame) sample lists. +ref_frame_samples_relative_to_pmag = { + pmag_ref_frame_plate_id: identity_samples(max_time_of_sequences_with_fixed_000), + optimised_ref_frame_plate_id: mantle_frame_samples_relative_to_pmag(optimised_sequence), + optimised_max_NNR_ref_frame_plate_id: mantle_frame_samples_relative_to_pmag(optimised_max_NNR_sequence), + no_net_rotation_ref_frame_plate_id: mantle_frame_samples_relative_to_pmag(no_net_rotation_sequence), + true_polar_wander_ref_frame_plate_id: true_polar_wander_samples_relative_to_pmag(optimised_sequence), +} + +# A rotation model of the reference frames relative to PMAG, used to re-express them relative to +# whichever frame the default (000) is set to: R(default->frame) = inverse[R(pmag->default)] * R(pmag->frame). +ref_frames_relative_to_pmag_features = [] +for ref_frame_plate_id, samples in ref_frame_samples_relative_to_pmag.items(): + if ref_frame_plate_id == pmag_ref_frame_plate_id: + continue # pmag is the anchor (identity relative to itself); no feature needed. + ref_frames_relative_to_pmag_features.append( + pygplates.Feature.create_total_reconstruction_sequence( + fixed_plate_id=pmag_ref_frame_plate_id, + moving_plate_id=ref_frame_plate_id, + total_reconstruction_pole=pygplates.GpmlIrregularSampling( + [pygplates.GpmlTimeSample(pygplates.GpmlFiniteRotation(rot), time, '', enabled) + for time, rot, enabled in samples]))) +ref_frames_relative_to_pmag_model = pygplates.RotationModel( + ref_frames_relative_to_pmag_features, default_anchor_plate_id=pmag_ref_frame_plate_id) + -def create_identity_rotation_sequence(rotation_description, begin_time, end_time=0.0): - rotation_time_samples = [] - identity_rotation = pygplates.FiniteRotation() - # Only need two time samples at start and end times. - for time in (end_time, begin_time): - rotation_time_samples.append( - pygplates.GpmlTimeSample( - pygplates.GpmlFiniteRotation(identity_rotation), - time, - rotation_description)) - return rotation_time_samples - -def clone_rotation_sequence(rotation_sequence, rotation_description): - rotation_time_samples = [] - for rotation_sample in rotation_sequence: - rotation_time_samples.append( - pygplates.GpmlTimeSample( - pygplates.GpmlFiniteRotation(rotation_sample.get_value().get_finite_rotation()), - rotation_sample.get_time(), - rotation_description, - rotation_sample.is_enabled())) - return rotation_time_samples - -def calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, rotation_description): - rotation_time_samples = [] - for optimised_rotation_sample in optimised_rotation_sequence: - time = optimised_rotation_sample.get_time() - # R(000->005). - optimised_finite_rotation = optimised_rotation_sample.get_value().get_finite_rotation() - # R(014->701). - finite_rotation_701_014 = original_rotation_model_anchored_014.get_rotation(time, 701) - # - # PMAG contain a mantle reference frame and true polar wander (TPW). - # Whereas OPT only contains a mantle reference frame. - # - # TPW + OPT = PMAG - # TPW = PMAG - OPT - # - # So looking at Africa this becomes... - # - # R(PMAG->701) - R(OPT->701) = R(PMAG->701) * R(701->OPT) - # = R(PMAG->OPT) - # - # ...and, using the old optimised rotation R(000->005) (where 000 was OPT and 005 was PMAG)... - # - # R(PMAG->OPT) = R(005->000) - # = inverse[R(000->005)] - # - # ...therefore... - # - # R(PMAG->701) - R(OPT->701) = R(PMAG->OPT) - # = inverse[R(000->005)] - # - # So we want 017->701 to be equivalent to the *inverse* of the optimised rotation 000->005: - # - # TPW = PMAG - OPT - # - # R(017->701) = R(PMAG->701) - R(OPT->701) - # = R(PMAG->OPT) - # = inverse[R(000->005)] - # - # And expressing this relative to 014 becomes... - # - # R(017->701) = inverse[R(000->005)] - # R(017->014) * R(014->701) = inverse[R(000->005)] - # R(017->014) = inverse[R(000->005)] * inverse[R(014->701)] - # - finite_rotation_true_polar_wander = optimised_finite_rotation.get_inverse() * finite_rotation_701_014.get_inverse() - rotation_time_samples.append( - pygplates.GpmlTimeSample( - pygplates.GpmlFiniteRotation(finite_rotation_true_polar_wander), - time, - rotation_description, - optimised_rotation_sample.is_enabled())) - return rotation_time_samples - -# Generate the rotation sequence for the PMAG reference frame. -# -# This depends on which reference frame 000 is assigned to. -if default_reference_frame_plate_id == pmag_ref_frame_plate_id: - # Both 000 and 014 are PMAG. - # - # So 000->014 is zero (identity). - rotation_time_samples_pmag = create_identity_rotation_sequence( - f' Reference frames: paleomag ({pmag_ref_frame_plate_id:03d} and 000)', - max_time_of_sequences_with_fixed_000) -elif default_reference_frame_plate_id == optimised_ref_frame_plate_id: - # 000 is optimised (and 014 is PMAG). - # - # So 000->014 is the existing optimised rotation (000->005). - rotation_time_samples_pmag = clone_rotation_sequence( - optimised_rotation_sequence, - f' Reference frames: paleomag ({pmag_ref_frame_plate_id:03d}) and optimised mantle (000)') -elif default_reference_frame_plate_id == no_net_rotation_ref_frame_plate_id: - # 000 is no-net rotation (and 014 is PMAG). - # - # So 000->014 is the existing no-net rotation (000->005). - rotation_time_samples_pmag = clone_rotation_sequence( - no_net_rotation_sequence, - f' Reference frames: paleomag ({pmag_ref_frame_plate_id:03d}) and no-net rotation (000)') -elif default_reference_frame_plate_id == true_polar_wander_ref_frame_plate_id: - # 000 is TPW (and 014 is PMAG). - # - # So 000->014 is the calculated TPW rotation sequence. - rotation_time_samples_pmag = calculate_true_polar_wander_rotation_sequence( - optimised_rotation_sequence, - f' Reference frames: paleomag ({pmag_ref_frame_plate_id:03d}) and approx true polar wander (000)') -else: - raise ValueError('Default reference frame plate ID is not one of the reference frame plate IDs') - -# Create a new PMAG rotation sequence R(000->014). -pmag_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( - fixed_plate_id=0, - moving_plate_id=pmag_ref_frame_plate_id, - total_reconstruction_pole=pygplates.GpmlIrregularSampling(rotation_time_samples_pmag)) - -# Default rotation model (anchor plate 000). -default_rotation_model = pygplates.RotationModel(output_rotation_feature_collections + [pmag_rotation_feature]) - - -def remove_pmag_from_rotation_sequence(rotation_sequence, rotation_description): - rotation_time_samples = [] - for rotation_sample in rotation_sequence: - time = rotation_sample.get_time() - finite_rotation = rotation_sample.get_value().get_finite_rotation() - # R(000->014). - finite_rotation_014_000 = default_rotation_model.get_rotation(time, pmag_ref_frame_plate_id) - # pid->014 = pid->000 * 000->014 - # pid->000 = pid->014 * inverse[000->014] - finite_rotation = finite_rotation * finite_rotation_014_000.get_inverse() - rotation_time_samples.append( - pygplates.GpmlTimeSample( - pygplates.GpmlFiniteRotation(finite_rotation), - time, - rotation_description, - rotation_sample.is_enabled())) - return rotation_time_samples - -# Generate the rotation sequences for the remaining reference frames (opt, NNR, TPW). -# -# This depends on which reference frame 000 is assigned to. -if default_reference_frame_plate_id == pmag_ref_frame_plate_id: - # Both 000 and 014 are PMAG. - # - # 015->000 is the existing optimised rotation (000->005). - rotation_time_samples_optimised = clone_rotation_sequence( - optimised_rotation_sequence, - f' Reference frames: optimised mantle ({optimised_ref_frame_plate_id:03d}) and paleomag (000)') - # 016->000 is the existing no-net rotation (000->005). - rotation_time_samples_no_net_rotation = clone_rotation_sequence( - no_net_rotation_sequence, - f' Reference frames: no-net rotation ({no_net_rotation_ref_frame_plate_id:03d}) and paleomag (000)') - # 017->000 is the calculated TPW rotation sequence. - rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence( - optimised_rotation_sequence, - f' Reference frames: approx true polar wander ({true_polar_wander_ref_frame_plate_id:03d}) and paleomag (000)') -elif default_reference_frame_plate_id == optimised_ref_frame_plate_id: - # 000 is optimised (and 014 is PMAG). - # - # 015->000 is zero (identity). - rotation_time_samples_optimised = create_identity_rotation_sequence( - f' Reference frames: optimised mantle ({optimised_ref_frame_plate_id:03d} and 000)', - max_time_of_sequences_with_fixed_000) - # 016->014 is the existing no-net rotation (000->005). - # 016->014 = 016->000 * 000->014 - # 016->000 = 016->014 * inverse[000->014] - rotation_time_samples_no_net_rotation = remove_pmag_from_rotation_sequence( - no_net_rotation_sequence, - f' Reference frames: no-net rotation ({no_net_rotation_ref_frame_plate_id:03d}) and optimised mantle (000)') - # 017->014 is the calculated TPW rotation sequence. - # 017->014 = 017->000 * 000->014 - # 017->000 = 017->014 * inverse[000->014] - rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, '') - rotation_time_samples_true_polar_wander = remove_pmag_from_rotation_sequence( - rotation_time_samples_true_polar_wander, - f' Reference frames: approx true polar wander ({true_polar_wander_ref_frame_plate_id:03d}) and optimised mantle (000)') -elif default_reference_frame_plate_id == no_net_rotation_ref_frame_plate_id: - # 000 is no-net rotation (and 014 is PMAG). - # - # 015->014 is the existing optimised rotation (000->005). - # 015->014 = 015->000 * 000->014 - # 015->000 = 015->014 * inverse[000->014] - rotation_time_samples_optimised = remove_pmag_from_rotation_sequence( - optimised_rotation_sequence, - f' Reference frames: optimised mantle ({optimised_ref_frame_plate_id:03d}) and no-net rotation (000)') - # 016->000 is zero (identity). - rotation_time_samples_no_net_rotation = create_identity_rotation_sequence( - f' Reference frames: no-net rotation ({no_net_rotation_ref_frame_plate_id:03d} and 000)', - max_time_of_sequences_with_fixed_000) - # 017->014 is the calculated TPW rotation sequence. - # 017->014 = 017->000 * 000->014 - # 017->000 = 017->014 * inverse[000->014] - rotation_time_samples_true_polar_wander = calculate_true_polar_wander_rotation_sequence(optimised_rotation_sequence, '') - rotation_time_samples_true_polar_wander = remove_pmag_from_rotation_sequence( - rotation_time_samples_true_polar_wander, - f' Reference frames: approx true polar wander ({true_polar_wander_ref_frame_plate_id:03d}) and no-net rotation (000)') -elif default_reference_frame_plate_id == true_polar_wander_ref_frame_plate_id: - # 000 is TPW (and 014 is PMAG). - # - # 015->014 is the existing optimised rotation (000->005). - # 015->014 = 015->000 * 000->014 - # 015->000 = 015->014 * inverse[000->014] - rotation_time_samples_optimised = remove_pmag_from_rotation_sequence( - optimised_rotation_sequence, - f' Reference frames: optimised mantle ({optimised_ref_frame_plate_id:03d}) and approx true polar wander (000)') - # 016->014 is the existing no-net rotation (000->005). - # 016->014 = 016->000 * 000->014 - # 016->000 = 016->014 * inverse[000->014] - rotation_time_samples_no_net_rotation = remove_pmag_from_rotation_sequence( - no_net_rotation_sequence, - f' Reference frames: no-net rotation ({no_net_rotation_ref_frame_plate_id:03d}) and approx true polar wander (000)') - # 017->000 is zero (identity). - rotation_time_samples_true_polar_wander = create_identity_rotation_sequence( - f' Reference frames: approx true polar wander ({true_polar_wander_ref_frame_plate_id:03d} and 000)', - max_time_of_sequences_with_fixed_000) -else: - raise ValueError('Default reference frame plate ID is not one of the reference frame plate IDs') - - -def invert_rotation_sequence(rotation_sequence): - rotation_time_samples = [] - for rotation_sample in rotation_sequence: - rotation_time_samples.append( - pygplates.GpmlTimeSample( - pygplates.GpmlFiniteRotation(rotation_sample.get_value().get_finite_rotation().get_inverse()), - rotation_sample.get_time(), - rotation_sample.get_description(), - rotation_sample.is_enabled())) - return rotation_time_samples - -# Create a new optimised rotation sequence R(015->000). # -# But we write it out as R(000->015), which is inverse of R(015->000), otherwise 000 is the moving plate ID and it -# will overlap in time with other sequences where 000 is also the moving plate ID (when loading into GPlates/pyGPlates). -optimised_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( - fixed_plate_id=0, - moving_plate_id=optimised_ref_frame_plate_id, - total_reconstruction_pole=pygplates.GpmlIrregularSampling( - invert_rotation_sequence(rotation_time_samples_optimised))) - -# Create a new no-net rotation sequence R(016->000). +# Build the output reference-frame features, each stored as a (fixed 000, moving ) sequence +# so that anchoring 000 gives the chosen default frame and anchoring gives that frame. # -# But we write it out as R(000->016), which is inverse of R(016->000), otherwise 000 is the moving plate ID and it -# will overlap in time with other sequences where 000 is also the moving plate ID (when loading into GPlates/pyGPlates). -no_net_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( - fixed_plate_id=0, - moving_plate_id=no_net_rotation_ref_frame_plate_id, - total_reconstruction_pole=pygplates.GpmlIrregularSampling( - invert_rotation_sequence(rotation_time_samples_no_net_rotation))) - -# Create a new true polar wander rotation sequence R(017->000). +# Stored finite rotation is R(000->frame) = R(default->frame) = inverse[R(pmag->default)] * R(pmag->frame). # -# But we write it out as R(000->016), which is inverse of R(016->000), otherwise 000 is the moving plate ID and it -# will overlap in time with other sequences where 000 is also the moving plate ID (when loading into GPlates/pyGPlates). -true_polar_wander_rotation_feature = pygplates.Feature.create_total_reconstruction_sequence( - fixed_plate_id=0, - moving_plate_id=true_polar_wander_ref_frame_plate_id, - total_reconstruction_pole=pygplates.GpmlIrregularSampling( - invert_rotation_sequence(rotation_time_samples_true_polar_wander))) - -# Add the reference frame features to a feature collection. -output_reference_frames_feature_collection = pygplates.FeatureCollection() -output_reference_frames_feature_collection.add(pmag_rotation_feature) -output_reference_frames_feature_collection.add(optimised_rotation_feature) -output_reference_frames_feature_collection.add(no_net_rotation_feature) -output_reference_frames_feature_collection.add(true_polar_wander_rotation_feature) - -# Store the reference frames in one of the output feature collections/files (we arbitrarily choose the first). +output_reference_frame_features = [] +for ref_frame_plate_id, samples_relative_to_pmag in ref_frame_samples_relative_to_pmag.items(): + + name = reference_frame_names[ref_frame_plate_id] + + if ref_frame_plate_id == pmag_ref_frame_plate_id: + # The PMAG feature carries R(000->014) = R(default->pmag) = inverse[R(pmag->default)]. + # Sample it at the default frame's native times (or identity span if default is PMAG). + if default_reference_frame_plate_id == pmag_ref_frame_plate_id: + sample_times = [t for t, _, _ in ref_frame_samples_relative_to_pmag[pmag_ref_frame_plate_id]] + else: + sample_times = [t for t, _, _ in ref_frame_samples_relative_to_pmag[default_reference_frame_plate_id]] + out_samples = [] + for time in sample_times: + pmag_to_default = ref_frames_relative_to_pmag_model.get_rotation( + time, default_reference_frame_plate_id) + out_samples.append((time, pmag_to_default.get_inverse(), True)) + else: + out_samples = [] + for time, rot_pmag_to_frame, enabled in samples_relative_to_pmag: + pmag_to_default = ref_frames_relative_to_pmag_model.get_rotation( + time, default_reference_frame_plate_id) + out_samples.append((time, pmag_to_default.get_inverse() * rot_pmag_to_frame, enabled)) + + if ref_frame_plate_id == default_reference_frame_plate_id: + description = ' Reference frames: {0} ({1:03d} and 000)'.format(name, ref_frame_plate_id) + else: + default_name = reference_frame_names[default_reference_frame_plate_id] + description = ' Reference frames: {0} ({1:03d}) and {2} (000)'.format( + name, ref_frame_plate_id, default_name) + + output_reference_frame_features.append( + pygplates.Feature.create_total_reconstruction_sequence( + fixed_plate_id=0, + moving_plate_id=ref_frame_plate_id, + total_reconstruction_pole=pygplates.GpmlIrregularSampling( + [pygplates.GpmlTimeSample(pygplates.GpmlFiniteRotation(rot), time, description, enabled) + for time, rot, enabled in out_samples]))) + + # -# Previously we added to a separate file called 'reference_frames.rot'. -# But that can be confusing for novice users who might think they can exclude that file if -# they don't want any of the reference frames (and hence just rely on the default 000 plate ID). -# But that would have been erroneous because 000 only existed in 'reference_frames.rot'. -# So now the reference frames are included within the normal output rotation files (associated with the input files). +# Write the output. The reference-frame features are bundled into the first output rotation file +# (rather than a separate 'reference_frames.rot') so that 000 is always defined wherever the model +# is loaded - a novice user cannot accidentally drop the file that defines plate 000. # -# Note: We could split the rotation sequences of the reference frames across the output rotation files -# (instead of dumping them all in just one of the files). For example, 0-1000 Ma into '1000_0_rotfile.rot' -# and 1000-1800 Ma into '1800_1000_rotfile.rot'. But we only have the filename to determine the time range -# of each rotation file. And that information isn't always in the filename. -# Besides, the original rotation sequences aren't always split up exactly according to the filename anyway. -output_reference_frames_feature_collection.add(output_rotation_feature_collections[0]) +output_reference_frames_feature_collection = pygplates.FeatureCollection(output_reference_frame_features) +for rotation_feature in output_rotation_feature_collections[0]: + output_reference_frames_feature_collection.add(rotation_feature) output_rotation_feature_collections[0] = output_reference_frames_feature_collection -# Write the output rotation filenames associated with the input rotation files. -for file_index in range(len(output_rotation_feature_collections)): - output_rotation_feature_collection = output_rotation_feature_collections[file_index] - output_rotation_filename = output_rotation_filenames[file_index] +for output_rotation_feature_collection, output_rotation_filename in zip( + output_rotation_feature_collections, output_rotation_filenames): output_rotation_feature_collection.write(output_rotation_filename) + print('Wrote {0}'.format(output_rotation_filename)) + +print('Reference frames added: pmag={0}, optimised={1}, optimised_max_NNR={2}, ' + 'no_net_rotation={3}, true_polar_wander={4} (default 000 = {5})'.format( + pmag_ref_frame_plate_id, optimised_ref_frame_plate_id, optimised_max_NNR_ref_frame_plate_id, + no_net_rotation_ref_frame_plate_id, true_polar_wander_ref_frame_plate_id, + reference_frame_names[default_reference_frame_plate_id]))