Variable Initial Crustal Thickness + Deforming Mesh Bugs

Hi. I am looking at crustal thickness changes, but I think currently you can only assign a constant initial crustal thickness for the entire deforming region? Is there any way you can have a variable initial crustal thickness and/or load in a grd/nc file with thickness values?

Also sometimes when I try to create a deforming mesh, I can add all of my lines to the boundary, but when I click create, no mesh appears. Is this a common bug, and are there specific types of features that must be added, or a certain way of saving the file, which can avoid this?


Hi Edward,

Currently you can only assign a constant initial thickness in GPlates. We plan to implement the importing of an initial thickness grid as you noted.

However in pyGPlates you can now do this and it will be available in the next release (after GPlates 2.3 is released, which is soon). You would still need to provide the point locations and initial thicknesses yourself (eg, sample a thickness grd/nc grid). Although before the release we will expose the ability to generate uniformly spaced points on the globe (ie, not just uniform lat/lon) inside a chosen deforming region (ie, polygon boundary). Which means you would just need to take those points and sample the initial thickness grid yourself.

It sounds like a problem with feature IDs of some sort. Have a look at this post, but that might not be the answer since you’re not saving and reloading. Even though the mesh does not appear, can you save the mesh to GPML and send to me along with the file(s) containing the features it references (ie, the ones you added to the mesh)?

You can add a regular feature (reconstructed by its plate ID) or reconstructed by half stage rotation (eg, mid-ocean ridge). But you cannot add a VirtualGeomagneticPole, SmallCircle, motion path or flowline.

Hi John,

Thanks for your help. For pygplates, are you referring to the MultiPointOnSphere function? I’m able to provide the point locations and thickness values, but how would I then connect this to a deforming mesh and reconstruct the thickness through time? Is there a way to create a reconstructed scalar coverage layer in pyGPlates and then load it into GPlates?

Okay, I’ve saved the files and emailed to you.

Best wishes,

Thanks for the files.

I noticed your feature IDs are all zero, which I’ve never seen before. That would explain the issue you are seeing. This is in your blocs.gpml file. Did you create it outside of GPlates/pyGPlates (as in create the XML some other way) ?

It’ll be a function that returns a MultiPointOnSphere, or a list of PointOnSphere, that is the uniform points inside the deforming network. But it won’t be available until the next release.

Have a look at this code from a version of pyBacktrack currently in development. The part that loads the initial scalars (in this case initial stretching factors) is:

initial_scalars = {pygplates.ScalarType.gpml_crustal_stretching_factor : reconstructed_crustal_stretching_factors}

…but this ability won’t be available until the next pyGPlates release.

Yes, I created the files in ArcGIS, loaded a .shp into GPlates and then saved it as a .gpml file.

And I will look forward to the next pyGPlates release! Any idea of the timescale for this?

Can you send me the Shapefile you created in ArcGIS? Thanks. And I guess the gpml file you saved was the blocs.gpml that I now have (if not could you send the GPML also).

I’m not really sure, but will know better once GPlates 2.3 has been released - even my estimate for that has not been very accurate :wink:

Thanks for sending the Shapefile.

I can see that you’ve mapped feature IDs to your Line_Numb Shapefile attribute which is zero for all your features, but GPlates requires feature IDs to be unique. This is because when the topology looks for a boundary section with feature ID 0 it will find more than one. This results in ambiguity so it chooses neither (this way it becomes more obvious that the topology is broken, rather than taking the risk of arbitrarily choosing a boundary section which could the wrong one).

I would suggest mapping your Blocs.shp again (you might need to delete the Blocs.shp.gplates.xml file to clear the current mapping) and load Blocs.shp back into GPlates again and map like this…

…noting that Feature ID: is left as <none>.

Then save it as Blocs.gpml. I know you’re not doing this, but for others if you keep loading the Shapefile (instead of GPML) then GPlates will create new feature IDs each time you load, and the topologies will not be able to find those features. This is something we need to look into when we overhaul our non-native file IO (eg, maybe we’ll create a new Shapefile attribute to store the feature ID in).

Then you can create topologies (unfortunately you’ll need to create them again).

Alternatively you could map Feature ID: to a Shapefile attribute that you know will be unique across all your features. But generally the feature ID should follow the string format that GPlates gives it (which starts with GPlates- and has alphanumeric characters and -, _ and . only). It’s best to allow GPlates to provide the feature ID (which it will do if you leave it mapped as <none>). And, as you have correctly done, we always recommend using GPML format when working with topologies (both the topologies and their referenced sections).

Great, thank you! Everything is working now.