Hi Florent,
I wrote a small script to test this out and can see a couple of potential issues.
Firstly, you’ll probably want to use GPML files when storing VirtualGeomagneticPole (VGP) features. We don’t really have a Shapefile attribute mapping for VGP properties (it’s a bit of a special case feature).
Secondly, GPML format is not an export format and so cannot be used as the output of pygplates.reconstruct()
.
So this first script shows how you would normally use pygplates.reconstruct()
to reconstruct (not to an exported file) and then use pygplates.reverse_reconstruct()
to reverse reconstruct that result. However this does not work well with VGP features for the reasons explained in the code comments below (essentially it’s because there’s two types of geometry in a VGP, the pole position and the site position, and currently only the former is reconstructed)…
import pygplates
initial_vgp_features = pygplates.FeatureCollection(r'C:\Users\John\Development\Usyd\gplates\sample_data\2.2\SampleData\FeatureCollections\Palaeomagnetism\gpml\Europe2004_RM_Npoles.vgp.gpml')
rotation_model = pygplates.RotationModel(r'C:\Users\John\Development\Usyd\gplates\sample_data\2.2\SampleData\FeatureCollections\Rotations\Matthews_etal_GPC_2016_410-0Ma_GK07.rot')
reconstruction_time = 100
# pygplates.reconstruct currently only reconstructs VGP *poles* (not the average sample site position).
#
# TODO: pygplates needs to add a new reconstructed geometry type called pygplates.ReconstructedVirtualGeomagneticPole
# (derived from pygplates.ReconstructedFeatureGeometry) that exposes the extra information (such as site position, inclination, declination, etc).
reconstructed_vgps = []
pygplates.reconstruct(initial_vgp_features, rotation_model, reconstructed_vgps, reconstruction_time)
vgp_features = []
for reconstructed_vgp in reconstructed_vgps:
# Clone the initial VGP feature - we're going to modify it by storing *reconstructed* geometries in it.
initial_vgp_feature = reconstructed_vgp.get_feature()
vgp_feature = initial_vgp_feature.clone()
# Set the reconstructed VGP pole.
vgp_feature.set_geometry(reconstructed_vgp.get_reconstructed_geometry(), pygplates.PropertyName.gpml_pole_position)
# ...the VGP pole position is the default geometry in a VGP feature, so the above line is the same as...
# vgp_feature.set_geometry(reconstructed_vgp.get_reconstructed_geometry())
# Set the reconstructed VGP site position (it's not the default geometry, so we need to extract it from
# pygplates.ReconstructedVirtualGeomagneticPole, but pygplates doesn't yet support that so the following call
# to 'reconstructed_vgp.get_average_sample_site_position()' isn't yet implemented and this won't work.
#
# TODO: Support pygplates.ReconstructedVirtualGeomagneticPole so that the following works.
#
# vgp_feature.set_geometry(reconstructed_vgp.get_average_sample_site_position(), pygplates.PropertyName.gpml_average_sample_site_position)
vgp_features.append(vgp_feature)
# We store *reconstructed* positions in the VGP feature, but features should store present day geometries, so rectify that.
# Note that 'pygplates.reverse_reconstruct()' correctly reverse reconstructs a VGP's pole position and site position.
# However the VGP site position was not actually reconstructed (or not actually stored in VGP feature above).
pygplates.reverse_reconstruct(vgp_features, rotation_model, reconstruction_time)
pygplates.FeatureCollection(vgp_features).write('vgp_features.gpml')
This second script works around that by explicitly reconstructing both the pole position and site position of each VGP. If you then load both the initial VGP file and the vgp_features.gpml file generated by this script into GPlates you’ll see they match up…
import pygplates
initial_vgp_features = pygplates.FeatureCollection(r'C:\Users\John\Development\Usyd\gplates\sample_data\2.2\SampleData\FeatureCollections\Palaeomagnetism\gpml\Europe2004_RM_Npoles.vgp.gpml')
rotation_model = pygplates.RotationModel(r'C:\Users\John\Development\Usyd\gplates\sample_data\2.2\SampleData\FeatureCollections\Rotations\Matthews_etal_GPC_2016_410-0Ma_GK07.rot')
reconstruction_time = 100
vgp_features = []
for initial_vgp_feature in initial_vgp_features:
# Clone the initial VGP feature - we're going to modify it by storing *reconstructed* geometries in it.
vgp_feature = initial_vgp_feature.clone()
# The rotation used to reconstructed both the VGP pole position and site position.
reconstruction_rotation = rotation_model.get_rotation(reconstruction_time, vgp_feature.get_reconstruction_plate_id())
# Reconstruct the VGP pole position and store in the VGP feature.
pole_position = vgp_feature.get_geometry(pygplates.PropertyName.gpml_pole_position)
reconstructed_pole_position = reconstruction_rotation * pole_position
vgp_feature.set_geometry(reconstructed_pole_position, pygplates.PropertyName.gpml_pole_position)
# Reconstruct the VGP site position and store in the VGP feature.
site_position = vgp_feature.get_geometry(pygplates.PropertyName.gpml_average_sample_site_position)
reconstructed_site_position = reconstruction_rotation * site_position
vgp_feature.set_geometry(reconstructed_site_position, pygplates.PropertyName.gpml_average_sample_site_position)
vgp_features.append(vgp_feature)
# We store *reconstructed* positions in the VGP feature, but features should store present day geometries, so rectify that.
# Note that 'pygplates.reverse_reconstruct()' correctly reverse reconstructs a VGP's pole position and site position.
pygplates.reverse_reconstruct(vgp_features, rotation_model, reconstruction_time)
pygplates.FeatureCollection(vgp_features).write('vgp_features.gpml')
Regards,
John