[GPlates-discuss] Reconstructing plate boundaries with PyGPlates

[Simon Williams]

Hi Adam,

Based on your sample code, a likely issue is related to the input files that you define.

Where you define the ‘features_name’, there is only one file name given, 'Global_EarthByte_230-0Ma_GK07_AREPS_PlateBoundaries.gpml'

However, for the topologies to be completely resolved from the AREPS file bundle there is another file called something like 'Global_EarthByte_230-0Ma_GK07_AREPS_Topology_BuildingBlocks.gpml’ which includes some extra elements (specifically topological line segments, which are kept in a separate file from the topological polygons for various reasons). If you don’t also point to this file when you tell your script to resolve the topologies, then the polygons that rely on the topological lines in the second file won’t get constructed correctly. And most of the topological lines in the current model are located in the Mediterranean region, so it makes sense that this is your issue. 
(I’d guess that you had this file loaded into GPlates when you made your screenshot - however, it is not entirely obvious sometimes how different layers in the GPlates layer manager relate to different source files, they often share the same names giving the illusion of a one-to-one relationship between layers and files, but this is not always the case)


Suggest you try the following:

define the features as a python list of filenames, 

features_name= ['../../reconstructions/Muller_etal_AREPS_Supplement/Global_EarthByte_230-0Ma_GK07_AREPS_PlateBoundaries.gpml’,
                            '../../reconstructions/Muller_etal_AREPS_Supplement/Global_EarthByte_230-0Ma_GK07_AREPS_Topology_BuildingBlocks.gpml’]

in theory you can pass this list of files idrectly to resolve_topologies:

pygplates.resolve_topologies(features_name, poles, resolved_topologies, reconstruction_t,shared_boundary_sections)


Alternatively if you want to avoid repetitive file io you could create a feature collection from the list of files with something like 

topology_features = pygplates.FeatureCollection()
for filename in features_name:
    topology_feature = pygplates.FeatureCollection(filename)
    topology_features.add(topology_feature)

then call resolve_topologies as you did before.

Try that and see if it makes the Mediterranean look better, if not it must be something else. Note also that you will likely get some ERROR messages anyway, since there are many small and hard-to-see topological errors all around the model at different times.

cheers,
Simon




> On 28 Oct 2017, at 6:12 AM, Adam Holt <adamholt at mit.edu> wrote:
> 
> Hi all,
> 
> I have a question about reconstructing plate boundaries using PyGPlates: I am using the Muller 2016 reconstruction to resolve plate topologies at 80 Ma, output boundaries, and then plot. My issue is that when I "resolve_topologies" in PyGPlates, I get errors for a few of the plates (error format: "ERROR: Failed to create a ResolvedTopologicalBoundary"), so that when I plot the boundaries there are gaps (e.g. in Mediterranean: See attached "PyGPlates" figure). When I do the same in GPlates (See attached "GPlates" figure), the plate boundaries are complete. Thus, it seems like I am doing something wrong with my simple PyGPlates script (copied below), and wondered if anyone could tell me what this is? Any help would be greatly appreciated!   
> 
> Thanks in advance!
> Adam Holt
> 
> 
> 
> --------------------------------------------------------------------------------------------------
> features_name='../../reconstructions/Muller_etal_AREPS_Supplement/Global_EarthByte_230-0Ma_GK07_AREPS_PlateBoundaries.gpml'
> poles_name='../../reconstructions/Muller_etal_AREPS_Supplement/Global_EarthByte_230-0Ma_GK07_AREPS.rot'
> ages_grd_name='/home/aholt/Dropbox/datasets/Muller_etal_2016_AREPS_Agegrids/netCDF_0-230Ma/agegrid_80.nc'
> reconstruction_t=80
> 
> topology_features = pygplates.FeatureCollection(features_name)
> poles = pygplates.RotationModel(poles_name);
> resolved_topologies = []; shared_boundary_sections = []
> pygplates.resolve_topologies(topology_features, poles, resolved_topologies, reconstruction_t,shared_boundary_sections)
> 
> plot_name='temp/temp_plot80Ma'
> subprocess.check_call(['./plot_plate_boundaries.sh',str(plot_name),ages_grd_name,'temp/points.txt','50000'])  # initialize gmt plot
> 
> for shared_boundary_section in shared_boundary_sections:
>    for shared_sub_segment in shared_boundary_section.get_shared_sub_segments():
>        points = shared_sub_segment.get_resolved_geometry().to_lat_lon_array()
>        np.savetxt('temp/points.txt',points[:,[1,0]],fmt='%.3f')
>        subprocess.check_call(['./plot_plate_boundaries.sh',str(plot_name),ages_grd_name,'temp/points.txt','1','blue']) # plot segments
> 
> subprocess.check_call(['./plot_plate_boundaries.sh',str(plot_name),ages_grd_name,'temp/points.txt','10000']) # finish plot
> ------------------------------------------------------------------------------------------------------
> <80Ma_plate-bounds_GPlates.png><80Ma_plate-bounds_PyGPlates.png>_______________________________________________
> GPlates-discuss mailing list
> GPlates-discuss at mailman.sydney.edu.au
> https://mailman.sydney.edu.au/mailman/listinfo/gplates-discuss