I have an issue with linking a reconstructed raster for the South China Sea region, that gives the wrong location in a polygon. The rotation file and raster data do not have an issue, i think, and the polygon position is at the right position, but the filled reconstruction raster does not show the suitable location, Please guide me if someone can understand.
Hi wajid,
Can you message me the files and I’ll take a look?
It could be something like a non-zero rotation at present day (usually the rotations in a rotation file should be zero at 0Ma).
Thank you for your quick and positive response. i am trying to attach the rotation file and gpml polygon with this message, but its not working. however i take a screenshot of the rotation file and its result. although you are right about the rotation file to be 0 at 0Ma, but when i try to change its value in rotation file then its misplaced. i can share my rotation file with you in your email, please send me your email address. here is mine, wajid.ouc@gmail.com
Thank you, John, for your email. I am sorry for the late response, I have attached the rotation and the necessary gpml files to your email, please find it.
waiting for your response.
Best regards,
Wajid
Good morning, professor, Until now, the email could not be delivered to you because of some issue with your email. can you guide me to send you those files via another email or any other way?
Hi Wajid, Thanks for sending the files. The issue is indeed the non-zero rotations at present day as we suspected. And I’ve sent you modified rotation and GPML files as a workaround.
I’ll cover the workaround procedure here in case others run into this issue in the future…
The issue arises during raster reconstruction when there are non-zero rotations at present day (0Ma). The static polygons, that are used to reconstruct the raster, will then get rotated to a new position, even at 0 Ma (ie, a different position than what’s stored in the static polygons file).
And raster reconstruction does not take this difference into account. Normally this is fine because it is highly recommended to only use rotation models that have zero rotations at present day. Ideally raster reconstruction would allow non-zero present-day rotations (perhaps in a future version of GPlates).
So the workaround involves modifying the rotations so that present-day rotations are zero, and also modifying the geometries in any reconstructable features (GPML files) using the rotation model to compensate.
The following pyGPlates code does that. You can copy it into a file, such as convert.py, then specify your own rotation files, and any reconstructable files that use those rotations, and then run python convert.py (after having installed pyGPlates). The output will be a modified copy of those rotation/reconstructable files but with a “converted_” prefix. Loading those files into GPlates (instead of the originals) will then result in correct raster reconstruction.
I should probably add this code somewhere, like PlateTectonicTools, but there are still a few limitations - the issue of crossovers (rotation sequences with same moving plate but different fixed plates) and the issue of a moving/fixed rotation sequence split across files (eg, 0-250Ma in one file and 250-410Ma in another) are not handled in the following code.
import os.path
import pygplates
###########################################################################################################
# Convert rotation features to have a zero rotation at present day for all plates.
#
# Also convert reconstructable features (ie, containing geometry) to use the converted rotation model.
###########################################################################################################
# List of input rotation files to convert.
#
# Output files will have "converted_" prefix.
rotation_files = [
'Global-rotation.rot',
]
# List of input reconstructable files to convert.
#
# Output files will have "converted_" prefix.
reconstructable_files = [
'900005-polygon-01.gpml',
'ERSCS-900007-ER-301-15-0Ma-copy.gpml',
'ET-900005-PhSP-608-14-11Ma-copy.gpml',
'EU-900007-EU-301-15-0Ma.gpml',
'SC-900001-ERU-301-15-0Ma-copy.gpml',
'SESC-900003-EU-301-15-0Ma-copy.gpml',
'SPhSP-900004-PhSP-608-15-0Ma-copy.gpml',
'WT-900001-ERU-301-15-0Ma-04.gpml',
]
def convert_rotation_features(
rotation_feature_collections):
"""
Function to make all moving/fixed plate rotation sequences have zero rotation at present day.
Returns a new list of modified feature collections.
"""
output_rotation_feature_collections = []
for rotation_feature_collection in rotation_feature_collections:
# Create an empty output rotation feature collection for each input rotation feature collection.
output_rotation_feature_collection = []
output_rotation_feature_collections.append(output_rotation_feature_collection)
# Iterate over rotation features in current collection.
for rotation_feature in rotation_feature_collection:
# Make sure it's a rotation feature (should be).
if not rotation_feature.get_total_reconstruction_pole():
continue
# Clone the input feature before we start making modifications.
# Otherwise we'll be modifying the input rotation feature (which the caller might not expect).
output_rotation_feature = rotation_feature.clone()
# Extract the moving/fixed rotation sequence from the rotation feature.
_, _, output_rotation_sequence = output_rotation_feature.get_total_reconstruction_pole()
# Get the rotation samples - enabled and disabled samples.
output_rotation_samples = output_rotation_sequence.get_time_samples()
if not output_rotation_samples:
# No time samples (shouldn't happen).
continue
# Iterate over time samples in current moving/fixed rotation sequence.
#
# First find the non-zero present-day rotation.
# Note that it's possible the sequence starts at a non-zero time,
# or that it does start at present day but already has a zero rotation.
# In both cases we won't modify the rotation sequence.
#
# Then (if there was a non-zero present-day rotation) adjust all finite rotations
# in the sequence at non-zero times (including disabled time samples).
present_day_rotation_adjustment = None
for time_sample in output_rotation_samples:
if not present_day_rotation_adjustment:
# See if it's the first enabled sample at present day (if current sequence has any).
if time_sample.is_enabled() and time_sample.get_time() == 0.0:
present_day_rotation = time_sample.get_value().get_finite_rotation()
# See if present-day rotation is non-zero.
if not present_day_rotation.represents_identity_rotation():
# Calculate the rotation adjustment from the non-zero present-day rotation.
present_day_rotation_adjustment = present_day_rotation.get_inverse()
# Set the present day rotation to the identity rotation.
time_sample.get_value().set_finite_rotation(pygplates.FiniteRotation())
else:
# Get the current rotation, adjust it for present day, and
# set the adjusted rotation back in the time sample.
#
# The finite rotation R(t) at current time t is composed of non-zero present-day rotation R(0)
# and the stage rotation between them R(0->t):
# R(t) = R(0->t) * R(0)
# Post-multiplying both sides by "inverse[R(0)]" gives us:
# R(t) * inverse[R(0)] = R(0->t) * R(0) * inverse[R(0)]
# = R(0->t)
# But we have set R(0) to the zero rotation, so our new finite rotation R'(t) is just:
# R'(t) = R(0->t)
# = R(t) * inverse[R(0)]
finite_rotation = time_sample.get_value().get_finite_rotation()
new_finite_rotation = finite_rotation * present_day_rotation_adjustment
time_sample.get_value().set_finite_rotation(new_finite_rotation)
output_rotation_feature_collection.append(output_rotation_feature)
# Return our output feature collections as a list of pygplates.FeatureCollection.
return [pygplates.FeatureCollection(rotation_feature_collection)
for rotation_feature_collection in output_rotation_feature_collections]
#
# Make all moving/fixed plate rotation sequences have zero rotation at present day.
#
# Includes propagating any non-zero rotation (at present day) into finite rotations at non-zero times.
#
# Read the input rotation feature collections.
input_rotation_feature_collections = [pygplates.FeatureCollection(input_rotation_filename)
for input_rotation_filename in rotation_files]
# Convert rotations.
output_rotation_feature_collections = convert_rotation_features(input_rotation_feature_collections)
# Write the modified rotation feature collections to disk.
for rotation_feature_collection_index in range(len(output_rotation_feature_collections)):
output_rotation_feature_collection = output_rotation_feature_collections[rotation_feature_collection_index]
# Each output filename is the input filename with a "converted_" prefix prepended.
input_rotation_filename = rotation_files[rotation_feature_collection_index]
dir, file_basename = os.path.split(input_rotation_filename)
output_rotation_filename = os.path.join(dir, '{0}{1}'.format('converted_', file_basename))
output_rotation_feature_collection.write(output_rotation_filename)
def convert_reconstructable_features(
reconstructable_feature_collections,
rotation_model_or_features):
"""
Function to reconstruct all reconstructable features to present day and store reconstructed geometry back into feature.
Returns a new list of modified feature collections.
"""
rotation_model = pygplates.RotationModel(rotation_model_or_features)
output_reconstructable_feature_collections = []
for reconstructable_feature_collection in reconstructable_feature_collections:
# Create an empty output reconstructable feature collection for each input reconstructable feature collection.
output_reconstructable_feature_collection = []
output_reconstructable_feature_collections.append(output_reconstructable_feature_collection)
# Iterate over reconstructable features in current collection.
for reconstructable_feature in reconstructable_feature_collection:
# Clone the input feature before we start making modifications.
output_reconstructable_feature = reconstructable_feature.clone()
# Get the feature's valid time period.
feature_begin_time, feature_end_time = output_reconstructable_feature.get_valid_time()
# Temporarily set the valid time period to all time in case feature is not valid at present day.
# We want to be able to reconstruct to present day.
output_reconstructable_feature.set_valid_time(pygplates.GeoTimeInstant.create_distant_past(), 0)
# Reconstruct feature to present day.
reconstructed_feature_geometries = []
pygplates.reconstruct(output_reconstructable_feature, rotation_model, reconstructed_feature_geometries, 0.0)
# Get present-day reconstructed geometries.
geometry_property_name = reconstructed_feature_geometries[0].get_property().get_name() # assume all geometries have same property name
feature_geometries = [reconstructed_feature_geometry.get_reconstructed_geometry()
for reconstructed_feature_geometry in reconstructed_feature_geometries]
# Set the new present-day geometries in the feature.
output_reconstructable_feature.set_geometry(feature_geometries, geometry_property_name)
# Restore the original valid time period.
output_reconstructable_feature.set_valid_time(feature_begin_time, feature_end_time)
output_reconstructable_feature_collection.append(output_reconstructable_feature)
# Return our output feature collections as a list of pygplates.FeatureCollection.
return [pygplates.FeatureCollection(reconstructable_feature_collection)
for reconstructable_feature_collection in output_reconstructable_feature_collections]
#
# Reconstruct all reconstructable features to present day and store reconstructed geometry back into feature.
#
# This ensures that the geometry stored in feature is the actual geometry at present day.
# Previously it would have been reconstructed with a non-zero present-day rotation, but after conversion will use zero present-day rotation.
#
# Read the input reconstructable feature collections.
input_reconstructable_feature_collections = [pygplates.FeatureCollection(input_reconstructable_filename)
for input_reconstructable_filename in reconstructable_files]
# Convert the reconstructable features.
output_reconstructable_feature_collections = convert_reconstructable_features(input_reconstructable_feature_collections, input_rotation_feature_collections)
# Write the modified reconstructable feature collections to disk.
for reconstructable_feature_collection_index in range(len(output_reconstructable_feature_collections)):
output_reconstructable_feature_collection = output_reconstructable_feature_collections[reconstructable_feature_collection_index]
# Each output filename is the input filename with a "converted_" prefix prepended.
input_reconstructable_filename = reconstructable_files[reconstructable_feature_collection_index]
dir, file_basename = os.path.split(input_reconstructable_filename)
output_reconstructable_filename = os.path.join(dir, '{0}{1}'.format('converted_', file_basename))
output_reconstructable_feature_collection.write(output_reconstructable_filename)
Thank you very much, John, for the workaround and step-by-step guidance. I am trying to follow your instructions and hope to build accurate rasters.
Best Regards
Wajid