226 unify the geometries of COBs in Muller2016

gplately/utils/convert_geometries.py

@@ -0,0 +1,103 @@
+#
+#    Copyright (C) 2024 The University of Sydney, Australia
+#
+#    This program is free software; you can redistribute it and/or modify it under
+#    the terms of the GNU General Public License, version 2, as published by
+#    the Free Software Foundation.
+#
+#    This program is distributed in the hope that it will be useful, but WITHOUT
+#    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+#    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+#    for more details.
+#
+#    You should have received a copy of the GNU General Public License along
+#    with this program; if not, write to Free Software Foundation, Inc.,
+#    51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+#
+
+import pygplates
+
+# work in progress!!
+
+
+def convert_polylines_to_polygons_within_feature(
+    feature: pygplates.Feature,
+    only_convert_closed_polylines=True,
+    verify_information_model=pygplates.VerifyInformationModel.yes,
+):
+    """TODO: describe what this funciton do and explain the parameters"""
+    # First step: gather all geometry property names (ignore other properties).
+    geometry_property_names = []
+    for property in feature:
+        # If it's a geometry property then add the property name to the list.
+        if property.get_value().get_geometry():
+            geometry_property_names.append(property.get_name())
+
+    # Second step: convert polylines to polygons (don't convert other geometry types).
+    for geometry_property_name in geometry_property_names:
+        # Get all geometries with the current property name.
+        #
+        # NOTE: I'm specifying the property name rather than relying on the default.
+        geometries_with_property_name = feature.get_geometries(geometry_property_name)
+
+        # Among the current geometries, convert any polylines to polygons (leave the other types alone).
+        for geometry_index, geometry in enumerate(geometries_with_property_name):
+            if isinstance(geometry, pygplates.PolylineOnSphere):
+                # As a bonus, only convert polyline if it is closed (ie, its first and last vertices are equal).
+                # Otherwise converting to a polygon will be dodgy (as Nicky pointed out).
+                # if only_convert_closed_polylines has been set to False by users, we assume the users know what they are doing and respect the users
+                if not only_convert_closed_polylines or geometry[0] == geometry[-1]:
+                    polygon = pygplates.PolygonOnSphere(geometry)
+                    geometries_with_property_name[geometry_index] = polygon
+
+        # Set all geometries with the current property name.
+        #
+        # NOTE: I'm specifying the property name rather than relying on the default.
+        feature.set_geometry(
+            geometries_with_property_name,
+            geometry_property_name,
+            verify_information_model=verify_information_model,
+        )
+
+
+def convert_polylines_to_polygons(
+    feature_collection: pygplates.FeatureCollection,
+    only_convert_closed_polylines=True,
+    verify_information_model=pygplates.VerifyInformationModel.yes,
+):
+    """convert all polylines in a given feature collection to polygons
+
+    Parameters
+    ----------
+    feature_collection: pygplates.FeatureCollection
+        all the polylines in this feature collection will be replaced with polygons
+    """
+    for feature in feature_collection:
+        convert_polylines_to_polygons_within_feature(
+            feature,
+            only_convert_closed_polylines=only_convert_closed_polylines,
+            verify_information_model=verify_information_model,
+        )
+
+
+def convert_polygons_to_polylines(feature_collection: pygplates.FeatureCollection):
+    """similar to the function above"""
+    # TODO:
+    # get the exterior ring of the polygon and convert it to polyline?
+    # Bianca, let's try to do this function
+    # step 1: create a new function "convert_polygons_to_polylines_within_feature"
+    #         the new function should be very similar to "convert_polylines_to_polygons_within_feature"
+    #         You may copy the code from "convert_polylines_to_polygons_within_feature" to the new function and make changes to the copy
+    # step 2: get the exterior ring and interior rings from the polygon and convert them to PolylineOnSphere
+    #         John may have better ideas on this?
+    #         see the links below for relavant pygplates doc
+    # https://www.gplates.org/docs/pygplates/generated/pygplates.polygononsphere#pygplates.PolygonOnSphere.get_exterior_ring_points
+    # https://www.gplates.org/docs/pygplates/generated/pygplates.polygononsphere#pygplates.PolygonOnSphere.get_interior_ring_points
+    # https://www.gplates.org/docs/pygplates/generated/pygplates.polylineonsphere#pygplates.PolylineOnSphere.__init__
+    # step 3: try to call this function in test_convert_geometries.py
+
+    #
+    # Please feel free to ask me or John if you got stuck at some point
+    # mostly, this is a training task. John and I are obliged to provide support.
+    #
+    pass

tests-dir/debug_utils/github_issue_214.py

@@ -0,0 +1,61 @@
+#!/usr/bin/env python3
+
+import cartopy.crs as ccrs
+import matplotlib.pyplot as plt
+from plate_model_manager import PlateModelManager
+
+import gplately
+
+print(gplately.__file__)
+
+# This file is used to debug github Issue #214.
+# https://github.com/GPlates/gplately/issues/214
+# Once the problem has been fixed, the exception should go away and the plot should be successful.
+
+
+def main():
+    # Call GPlately's PlateModelManager object and request data from the Müller et al. 2016 study
+    pm_manager = PlateModelManager()
+    muller2016_model = pm_manager.get_model("Muller2016", data_dir=".")
+    rotation_model = muller2016_model.get_rotation_model()
+    topology_features = muller2016_model.get_topologies()
+    static_polygons = muller2016_model.get_static_polygons()
+
+    model = gplately.PlateReconstruction(
+        rotation_model, topology_features, static_polygons
+    )
+
+    # Obtain features for the PlotTopologies object with PlateModelManager
+    coastlines = muller2016_model.get_layer("Coastlines")
+    continents = None
+    COBs = muller2016_model.get_layer("COBs")
+
+    # Call the PlotTopologies object
+    gplot = gplately.plot.PlotTopologies(
+        model, coastlines=coastlines, continents=continents, COBs=COBs
+    )
+
+    gplot.time = 10  # Ma
+
+    fig = plt.figure(figsize=(8, 4))
+    ax = fig.add_subplot(111, projection=ccrs.Mollweide(190))
+
+    ax.set_global()
+    ax.gridlines(
+        crs=ccrs.PlateCarree(),
+        draw_labels=True,
+        linewidth=2,
+        color="gray",
+        alpha=0.5,
+        linestyle="--",
+    )
+
+    gplot.plot_continent_ocean_boundaries(ax, color="red")
+
+    plt.title(f"{gplot.time} Ma")
+
+    plt.show()
+
+
+if __name__ == "__main__":
+    main()

tests-dir/unittest/test_convert_geometries.py

@@ -0,0 +1,75 @@
+#!/usr/bin/env python3
+
+import pygplates
+
+import gplately
+
+print(gplately.__file__)
+
+from gplately.utils import convert_geometries
+
+# Bianca, finish this file to test the new functions
+
+if __name__ == "__main__":
+
+    filename = "Global_230-0Ma_GK07_AREPS_COB_Terranes.gpml"
+    feature_collection = pygplates.FeatureCollection(filename)
+    geometry_types = []
+    for feature in feature_collection:
+        geoms = feature.get_all_geometries()
+        for geometry in geoms:
+            type_name = type(geometry).__name__
+            if type_name not in geometry_types:
+                geometry_types.append(type_name)
+
+    print("feature types before conversion: ", geometry_types)
+
+    # Bianca, let's create a .gpml file and put all PolylineOnSphere geometries into the .gpml file
+    # we would like to see what these PolylineOnSphere geometries look like
+    # step 1: create a new FeatureCollection, something like "new_feature_collection = pygplates.FeatureCollection()"
+    # step 2: loop through the original feature collection, something like "for feature in feature_collection:"
+    # step 3: for each feature, get all geometries, something like "geoms = feature.get_all_geometries()"
+    # step 4: loop through all geometries with an inner loop, something like "for geometry in geoms:"
+    # step 5: check if the geometry is a PolylineOnSphere, something like "if isinstance(geometry, pygplates.PolylineOnSphere):"
+    # step 6: for each PolylineOnSphere geometries, create a new feature, something like "new_feature  = pygplates.Feature()"
+    # step 7: put the PolylineOnSphere geometry into the new feature, something like "new_feature.set_geometry(geometry)"
+    # step 8: add the new feature into the new feature collection, something like "new_feature_collection.add(new_feature)"
+    # step 9: save the new feature collection to a .gpml file
+    # step 10: open the .gpml file in GPlates desktop software. https://www.earthbyte.org/download-gplates-2-5/
+
+    #############
+    new_feature_collection = pygplates.FeatureCollection()
+
+    for feature in feature_collection:
+        geoms = feature.get_all_geometries()
+        for geometry in geoms:
+            if isinstance(geometry, pygplates.PolylineOnSphere):
+                new_feature = pygplates.Feature()
+                new_feature.set_geometry(geometry)
+                new_feature_collection.add(new_feature)
+
+    new_feature_collection.write("polylines-in-muller2016.gpml")
+    #############
+
+    convert_geometries.convert_polylines_to_polygons(
+        feature_collection,
+        only_convert_closed_polylines=False,
+        verify_information_model=pygplates.VerifyInformationModel.no,
+    )
+
+    # now let's see if the geometries have been converted
+    geometry_types = []
+    for feature in feature_collection:
+        geoms = feature.get_all_geometries()
+        for geometry in geoms:
+            type_name = type(geometry).__name__
+            # print(type_name)
+            # if type_name == "PolylineOnSphere":
+            #    print(geoms)
+            if type_name not in geometry_types:
+                geometry_types.append(type_name)
+
+    print("feature types after conversion: ", geometry_types)
+
+    # Write the modified feature collection to output file.
+    feature_collection.write("converted-cob-muller2016.gpmlz")