Inconsistent result when reconstructing a feature with rotaion angle and the function "reconstruct"

[hujs]

Hi all,

I have met a problem when reconstructing a feature using the Muller et al. 2019 model.

I did the reconstruction using the “reconstruct” function as below:
reconstructed_geometry = []
pygplates.reconstruct(feature_to_be_reconstructed, rotation_model, reconstructed_geometry,
reconstruction_time=time)
print reconstructed_geometry[0].get_reconstructed_geometry().to_lat_lon()

I also tried a different way:
plateid = feature_to_be_reconstructed.get_reconstruction_plate_id()
rotation_angle = rotation_model.get_rotation(time, plateid, 0, 0)
print (rotation_angle*(assigned_point_features_10[i].get_geometry())).to_lat_lon()

I got two different results. This is confusing. I thought the two methods are equivalent. Am I missing something?

Let me know if anyone can help me with this. Thanks!

[john.cannon]

How is assigned_point_features_10 related to feature_to_be_reconstructed? Are they the same?

[hujs]

Yes. They are the same.

The initial purpose is to partition a point feature into plates to get the plate id. The point was at (4.322765776364536, -172.00856066713223) at 130 Ma. And the partitioning function gives a plate id of 926 and its present-day location at (45.601579512218066, 30.569660725405726).

When I reconstruct the point back to 130 Ma using the pygplates.reconstruct function, I get the same location as the initial point. But we I compute the rotation angle using the time 130 Ma and the plate id 926 and then rotate it back to 130 Ma, I get a different result.

[john.cannon]

@hujs

Yes. They are the same.

I suspect they might have different geometries (even though they represent the same feature).

@hujs

The point was at (4.322765776364536, -172.00856066713223) at 130 Ma. And the partitioning function gives a plate id of 926 and its present-day location at (45.601579512218066, 30.569660725405726).

It sounds like you’re using a non-zero (ie, 130Ma) reconstruction time in the partitioning function. In which case I’m assuming feature_to_be_reconstructed contains that point at 130Ma (4.322765776364536, -172.00856066713223) and assigned_point_features_10 is the output of the partitioning function which contains the present day point location (45.601579512218066, 30.569660725405726).

Is that right? If so, then they contain different geometries. If not, then can you post some more code (such as how you partition)?

And maybe confirm whether the feature geometries are the same with:

print feature_to_be_reconstructed[0].get_geometry().to_lat_lon()
print assigned_point_features_10[i].get_geometry().to_lat_lon()

[hujs]

Blockquote
It sounds like you’re using a non-zero (ie, 130Ma) reconstruction time in the partitioning function. In which case I’m assuming feature_to_be_reconstructed contains that point at 130Ma (4.322765776364536, -172.00856066713223 ) and assigned_point_features_10 is the output of the partitioning function which contains the present day point location (45.601579512218066, 30.569660725405726 ).
Blockquote
Yes, exactly. But I could not understand why they contain different geometries. Below is a clearer demo.

[john.cannon]

Thanks for the clarification, I see what you mean now.

So the problem, as I now understand it, is you’re reconstructing assigned_point_feature[0] two different ways and getting two different results. I had previously thought you were comparing reconstructions of assigned_point_feature[0] (containing present day geometry) with reconstructions of point_feature (containing geometry already reconstructed to 130Ma). But that’s not the case here.

In that case, I think the problem is due to plate 926 having a non-zero rotation at present day (relative to anchor plate 0). And since you are explicitly specifying a from_time of 0 in:

rotation_angle = rotation_model.get_rotation(time, pid, 0, 0)

…(ie, the first 0) then the non-zero present day rotation is excluded from the returned rotation. In other words, it’s returning the rotation of plate 926 at 130Ma relative to 0Ma (which excludes the fact that the geometry stored in your feature also undergoes a non-zero rotation just to get to 0Ma).

Instead try using:

rotation_angle = rotation_model.get_rotation(time, pid)

…and see if you now get the same reconstructed results.

The reason for this distinction is explained at the bottom of RotationModel.get_rotation() in the following note:

Explicitly setting from_time to zero can give a different result than not specifying from_time at all if the moving plate (or fixed plate) has a non-zero finite rotation at present day (relative to the anchor plate). However all present-day finite rotations should ideally be zero (identity), so typically there should not be a difference.

…and so it’s probably better to avoid non-zero present day rotations whenever possible as it usually leads to confusing results.

The above note assumes you’re using the latest pyGPlates public release (revision 28) which introduced that distinction to avoid inconsistencies that started appearing whenever there were non-zero present day rotations.

[hujs]

Thank you, John.

I can understand the issue now. I just made a test to delete the zeros in function call. The results are consistent now.