Hi all
I’m searching how to unreconstruct, manually and by Python, drawn feature for one time.
For instance : I drawn an entity at 115Ma and I want to know where is located at Present Day.
But location of my drawn entity is considered at Present Day and not at 115Ma.
No doubt I forgot several things. Entity was drawn with an other GIS tools.
Only fields informed of attribute table : PLATEID, FROMAGE, TOAGE
Thanks for your help.
Hi Fabrice,
If you have geometry at 115Ma and know which plate it’s on you can just use pygplates.RotationModel.get_rotation() and then reverse the returned rotation with pygplates.FiniteRotation.get_inverse() and then multiply your geometry by that reversed rotation (as shown in the example code documented in that FiniteRotation link).
If you don’t know the plate then it’s a little more involved but I can show you?
Also this assumes you can get your geometry out of whatever GIS format it’s in. If it’s a format supported by GPlates then it’s relatively straightforward using pyGPlates.
Regards John
Hi John
Thank you
Unfortunately our unreconstruct script doesn’t work.
For our exercice first we reconstruct our object and then we’ve tryed to unreconstruct.
Reconstruct operation it’s ok. No the opposite.
It doesn’t work but no error message
rotation_features = pygplates.FeatureCollection(rot_file)
vgp_features = pygplates.FeatureCollection(shp_to_rotate)
pygplates.reverse_reconstruct(vgp_features, rotation_features, age_r)
Regards
Hi Fabrice,
Your vgp_features will have their present day geometries modified by the call to pygplates.reverse_reconstruct(). It doesn’t generate any new output (like pygplates.reconstruct()), instead it modifies the features you pass into it as follows…
Before calling pygplates.reverse_reconstruct(), your vgp_features are assumed to have geometries representing a past time (in your case age_r). After calling that function, your vgp_features should have geometries representing present day. So perhaps double-check that they have in fact been modified.
The reason for this is features are always meant to store geometries (internally) in their present day position, however when that’s not the case then the pygplates.reverse_reconstruct() function comes in handy (that’s really the scenario it’s designed for), in which case it rectifies the situation by replacing the feature’s internal non-present-day geometries with a present-day version (by reverse reconstructing using the feature’s plate ID - in your case you have a Shapefile and so the usual PLATEID1 attribute will be used for that). Then once that’s done you are free to reconstruct that feature to any past reconstruction time.
Regards,
John
Hi John,
I work with Fabrice on the subject, i have wrote the script. It appears that the function don’t change any geometry in the input shapefile (i have tried on 60 shapefiles polyline, polygon or points).
I have tried 2 scenarii,
in first: i apply different Plate ID to a shapefile, representing a past time, and i try the reverse reconstruct function but nothing change.
So to be sure that’s plate model and rotation model works,
in a second try, i have successfully reconstructed a shapefile in past time with a rotation model and plate model then from this reconstructed shapefile (with same rotation model) i have try to make the reverse reconstruct method to come back at the initial state but as in the first scenario nothing changes.
With the pygplates.reconstruct() function, i have already seen an incompability with our arcmap python library’s version but with some adjustment i bypass it. It’s possible that’s it’s the same with the pygplates.reverse_reconstruct() function, but i have tried outside arcpy build and nothing change, so i am pretty stuck !
Regards,
Florent
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
Hi John,
No, i use directly Shapefile as initial_vgp_feature, the easiest way as we use mainly ESRI solution for our GIS data.
shp_with_plate_id = r"C:\Local\00_InputShp\Input_test.shp" # Input Shapefile with columns [PLATEID1, FROMAGE, TOAGE]
shp_rotated = r"C:\Local\01_OutputShp\Output_test.shp"
rot_file = r"C:\Local\T_ROT_GEM_jan21.rot"
age_r = 100
1 RECONSTRUCT THE SHAPEFILE
Load FC to rotate (if shapefile format, the pygplate feature collection can contain only 1 file)
vgp_features = pygplates.FeatureCollection(shp_with_plate_id)
Load rotation file as pygplate feature collection
rotation_features = pygplates.FeatureCollection(rot_file)
pygplates.reconstruct(vgp_features, rotation_features, shp_rotated, age_r)
del vgp_features, rotation_features2 REVERSE RECONSTRUCT THE RECONSTRUCTED SHAPEFILE
shp_to_derotate = shp_rotated # The shapefile to reverse reconstruct, is the shapefile reconstructed in first step
Load rotation file as pygplate feature collection
rotation_features = pygplates.FeatureCollection(rot_file)
Load FC to rotate (if shapefile format, the pygplate feature collection can contain only 1 file)
vgp_features = pygplates.FeatureCollection(shp_to_derotate)
pygplates.reverse_reconstruct(vgp_features, rotation_features, age_r)
del vgp_features, rotation_features
Hi Florent,
Thanks for the code, yeah that’s pretty much what I figured you were doing. That is, loading from Shapefile and also reconstructing directly to Shapefile.
You’ll want to use GPML instead of Shapefile (as noted above) since you’re using VGP features. Unfortunately that makes interoperability a bit more difficult.
How are you actually creating the VGP features in the first place? Ultimately it sounds like you just need to reverse reconstruct them (I assume the reconstructing was just for test purposes). So if you run pygplates.reverse_reconstruct() on a GPML file (containing VGPs) then it should work (assuming the VGPs in the GPML file are OK). Are you wanting to display the present day VGP site position or pole position or both in ArcGIS (or whichever GIS you’re using)?
Regards,
John
Hello John,
Sorry for this response period.
We just want to reverse reconstruct geometric files.
We want to display the present day VGP site position in ArcGIS. As feature collection works with shapefile for the reconstruct method, i suppose the reverse method was able to run too.
So the only way is to use GPML File but can we easely convert ESRI Shapefile format to GPML format ?
Regards,
Florent
Hi Florent,
For regular features yes, but for VGP features no.
How are you actually creating the VGP features in the first place? From this post of yours I can see you are loading a Shapefile that appears to contain VGP features (because you call it vgp_features). So I’m just wondering how you got VGP features into the Shapefile in the first place.
Regards,
John
Re,
In first place, i have intersect my original files with plate model to create the shapefile “shp_with_plate_id” from a spatial join then i run the reconstruct method. To do the reverse reconstruct, i do the same but before i reconstruct the plate model in the age position then try the reverse.reconstruct() at this age.
Florent
Hi Florent,
I think the feature types in your shp_with_plate_id file are VirtualGeomagneticPole (this will be the value of the GPGIM_TYPE shapefile attribute).
Probably the easiest workaround at this stage is to change the feature type to something else (like UnclassifiedFeature). Then you can still use Shapefiles and everything should work fine (reconstructing and reverse reconstructing).
The problem is really due to the VirtualGeomagnetic feature type. Usually it has both a site position and a pole position. Since you’re only interested in the VGP site position, your original files might have just contained the site position (and not the pole position). However when GPlates/pyGPlates reads the Shapefile it sees a single position (in each feature) and assumes it’s the pole position. So by changing the feature type to something other than a VGP feature I think you can bypass this issue.
Regards,
John
Hello John,
Yes, it’s works !
Thank you very much, i has thinked i could’nt directly to reference the shapefile path as input parameter for reverse_reconstruct() function because few time ago, with reconstruct() function, i have an issue it’s why i use the featurecollection() function before.
Thanks, for your help !
Florent Montarou