244 experiment pygmt plotting

docker/env.yaml

@@ -7,3 +7,5 @@ dependencies:
   - gplately 
   - jupyter
   - moviepy
+  - gmt
+  - pygmt

gplately/auxiliary.py

@@ -0,0 +1,56 @@
+from typing import Union
+
+from plate_model_manager import PlateModel, PlateModelManager
+
+from .mapping.plot_engine import PlotEngineType
+from .plot import PlotTopologies
+from .reconstruction import PlateReconstruction
+
+
+def get_gplot(
+    model_name: str,
+    model_repo_dir: str,
+    age: Union[int, float],
+    plot_engine: PlotEngineType = PlotEngineType.CARTOPY,
+) -> PlotTopologies:
+    """auxiliary function to get gplot object"""
+    try:
+        model = PlateModelManager().get_model(model_name, data_dir=model_repo_dir)
+    except:
+        model = PlateModel(model_name, data_dir=model_repo_dir, readonly=True)
+
+    if model is None:
+        raise Exception(f"Unable to get model ({model_name})")
+
+    topology_features = None
+    static_polygons = None
+    coastlines = None
+    COBs = None
+    continents = None
+
+    all_layers = model.get_avail_layers()
+
+    if "Topologies" in all_layers:
+        topology_features = model.get_layer("Topologies")
+    if "StaticPolygons" in all_layers:
+        static_polygons = model.get_layer("StaticPolygons")
+    if "Coastlines" in all_layers:
+        coastlines = model.get_layer("Coastlines")
+    if "COBs" in all_layers:
+        COBs = model.get_layer("COBs")
+    if "ContinentalPolygons" in all_layers:
+        continents = model.get_layer("ContinentalPolygons")
+
+    m = PlateReconstruction(
+        model.get_rotation_model(),
+        topology_features=topology_features,
+        static_polygons=static_polygons,
+    )
+    return PlotTopologies(
+        m,
+        coastlines=coastlines,
+        COBs=COBs,
+        continents=continents,
+        time=age,
+        plot_engine=plot_engine,
+    )

gplately/mapping/__init__.py

@@ -0,0 +1 @@
+# This submodule contains code to plot maps.

gplately/mapping/plot_engine.py

@@ -0,0 +1,46 @@
+#
+#    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.
+#
+
+from enum import Enum
+from abc import ABC, abstractmethod
+
+from geopandas.geodataframe import GeoDataFrame
+
+
+class PlotEngineType(Enum):
+    CARTOPY = 1
+    PYGMT = 2
+
+
+class PlotEngine(ABC):
+    @abstractmethod
+    def plot_geo_data_frame(self, gdf: GeoDataFrame, **kwargs):
+        pass  # This is an abstract method, no implementation here.
+
+    @abstractmethod
+    def plot_pygplates_features(self, features, **kwargs):
+        pass  # This is an abstract method, no implementation here.
+
+    @abstractmethod
+    def plot_subduction_zones(
+        self,
+        gdf_subduction_left: GeoDataFrame,
+        gdf_subduction_right: GeoDataFrame,
+        color="blue",
+        **kwargs,
+    ):
+        pass  # This is an abstract method, no implementation here.

gplately/mapping/pygmt_plot.py

@@ -0,0 +1,160 @@
+#
+#    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.
+#
+from geopandas.geodataframe import GeoDataFrame
+import pygmt
+from plot_engine import PlotEngine
+
+pygmt.config(
+    FONT_ANNOT=8,
+    FONT_LABEL=8,
+    FONT=8,
+    MAP_TICK_PEN="0.75p",
+    MAP_FRAME_PEN="0.75p",
+    MAP_TICK_LENGTH_PRIMARY="4p",
+)
+
+# ----- parameters for plot
+region = "d"
+width = 10
+projection = "N180/"
+x_offset = width + 2
+
+# plate boundary stuff
+plateboundary_width = "0.5p"
+age_font = "12p,Helvetica,black"
+label_font = "12p,Helvetica,black"
+label_offset = "j0/-0.5c"
+label_position = "TC"
+
+
+class PygmtPlotEngine(PlotEngine):
+    def __init__(self, projection="N180/10c", region="d"):
+        self.fig = pygmt.Figure()
+        self.fig.basemap(region=region, projection=projection, frame="lrtb")
+
+    def plot_geo_data_frame(self, gdf: GeoDataFrame, **kwargs):
+        plot_geo_data_frame(self.fig, gdf, **kwargs)
+
+    def plot_pygplates_features(self, features, **kwargs):
+        pass
+
+    def plot_subduction_zones(
+        self,
+        gdf_subduction_left: GeoDataFrame,
+        gdf_subduction_right: GeoDataFrame,
+        color="blue",
+        **kwargs,
+    ):
+        plot_subduction_zones(
+            self.fig, gdf_subduction_left, gdf_subduction_right, color=color, **kwargs
+        )
+
+
+def get_pygmt_basemap_figure(projection="N180/10c", region="d"):
+    fig = pygmt.Figure()
+    fig.basemap(region=region, projection=projection, frame="lrtb")
+    return fig
+
+
+def plot_subduction_zones(
+    fig: pygmt.Figure,
+    gdf_subduction_left: GeoDataFrame,
+    gdf_subduction_right: GeoDataFrame,
+    color="blue",
+    **kwargs,
+):
+    fig.plot(
+        data=gdf_subduction_left, pen=f"0.5p,{color}", fill=color, style="f0.2/0.08+l+t"
+    )
+    fig.plot(
+        data=gdf_subduction_right,
+        pen=f"0.5p,{color}",
+        fill=color,
+        style="f0.2/0.08+r+t",
+    )
+
+
+def plot_geo_data_frame(fig: pygmt.Figure, gdf: GeoDataFrame, **kwargs):
+    line_width = "0.1p"
+    line_color = "blue"
+
+    if "edgecolor" in kwargs.keys():
+        if isinstance(kwargs["edgecolor"], str):
+            line_color = kwargs["edgecolor"]
+        else:
+            raise Exception(
+                "The edgecolor parameter is not string. Currently, the pygmt plot engine only supports colour name."
+            )
+
+    if "linewidth" in kwargs.keys():
+        line_width = f"{kwargs['linewidth']}p"
+
+    fill = None
+    if "facecolor" in kwargs.keys() and kwargs["facecolor"].lower() != "none":
+        fill = f"{kwargs['facecolor']}"
+
+    if line_color.lower() == "none":
+        line_width = "0"
+        line_color = fill
+
+    if "fill" in kwargs.keys():
+        fill = kwargs["fill"]
+
+    if "pen" in kwargs.keys():
+        pen = kwargs["pen"]
+    else:
+        pen = f"{line_width},{line_color}"
+
+    style = None
+    if "style" in kwargs.keys():
+        style = kwargs["style"]
+
+    label = None
+    if "gmtlabel" in kwargs.keys():
+        label = kwargs["gmtlabel"]
+
+    fig.plot(
+        data=gdf.geometry, pen=pen, fill=fill, style=style, transparency=0, label=label
+    )
+
+    """
+    fig.plot(data=gdf_coastlines, fill=coastline_color, frame=["xa0", "ya0"],  transparency=0)
+
+    fig.plot(data=gdf_topo_plates.geometry, pen='%s,%s' % (plateboundary_width, plate_colour), frame="lrtb")
+    fig.plot(data=gdf_subduction_left, pen='%s,%s' % (plateboundary_width, subduction_zone_colour), fill=subduction_zone_colour, style='f0.2/0.08+l+t')
+    fig.plot(data=gdf_subduction_right, pen='%s,%s' % (plateboundary_width, subduction_zone_colour), fill=subduction_zone_colour, style='f0.2/0.08+r+t')
+    fig.plot(data=gdf_ridges_transforms, pen='%s,%s' % (plateboundary_width, ridge_colour))
+    fig.plot(data=gplot.get_transforms(), pen='%s,%s' % (plateboundary_width, transform_color))
+
+    fig.text(text='gplot.get_transforms(): %s Ma' % age, position=label_position, no_clip=True, font=label_font, offset=label_offset)
+
+    fig.shift_origin(xshift=x_offset)
+    fig.basemap(region=region, projection="%s%sc" % (projection, width), frame="lrtb")
+    fig.plot(data=gdf_cobs, fill=COB_color, transparency=0, )
+    fig.plot(data=gdf_coastlines, fill=coastline_color, frame=["xa0", "ya0"],  transparency=0)
+
+    fig.plot(data=gdf_topo_plates.geometry, pen='%s,%s' % (plateboundary_width, plate_colour), frame="lrtb", label='other plate boundary types')
+    fig.plot(data=gdf_subduction_left, pen='%s,%s' % (plateboundary_width, subduction_zone_colour), fill=subduction_zone_colour, style='f0.2/0.08+l+t', label='subduction zones')
+    fig.plot(data=gdf_subduction_right, pen='%s,%s' % (plateboundary_width, subduction_zone_colour), fill=subduction_zone_colour, style='f0.2/0.08+r+t')
+    fig.plot(data=gdf_ridges_transforms, pen='%s,%s' % (plateboundary_width, ridge_colour), label='ridges and transforms')
+
+    # from gpml: transforms
+    fig.plot(data=gdf_topo_transforms, pen='%s,%s' % (plateboundary_width, transform_color), label = 'transforms')
+    fig.text(text='FeatureType.gpml_transform: %s Ma' % age, position=label_position, no_clip=True, font=label_font, offset=label_offset)
+
+    fig.legend(position='jBL+o-2.7/0', box="+gwhite+p0.5p")
+    """

gplately/plot.py

@@ -45,6 +45,11 @@
     validate_topology_availability,
 )
 from .gpml import _load_FeatureCollection
+from .mapping.plot_engine import PlotEngineType
+from .mapping.pygmt_plot import (
+    plot_geo_data_frame,
+    plot_subduction_zones as pygmt_plot_subduction_zones,
+)
 from .pygplates import FeatureCollection as _FeatureCollection
 from .reconstruction import PlateReconstruction as _PlateReconstruction
 from .tools import EARTH_RADIUS
@@ -298,12 +303,13 @@ def __init__(
         COBs=None,
         time=None,
         anchor_plate_id=0,
+        plot_engine: PlotEngineType = PlotEngineType.CARTOPY,
     ):
         self.plate_reconstruction = plate_reconstruction
 
         if self.plate_reconstruction.topology_features is None:
             self.plate_reconstruction.topology_features = []
-            logger.warn("Plate model does not have topology features.")
+            logger.warning("Plate model does not have topology features.")
 
         self.base_projection = ccrs.PlateCarree()
 
@@ -321,6 +327,7 @@ def __init__(
         self._topologies = None
 
         self._anchor_plate_id = self._check_anchor_plate_id(anchor_plate_id)
+        self._plot_engine = plot_engine
 
         # store topologies for easy access
         # setting time runs the update_time routine
@@ -722,16 +729,24 @@ def _plot_feature(self, ax, get_feature_func, **kwargs):
             tessellate_degrees=tessellate_degrees,
         )
 
+        if not isinstance(gdf, gpd.GeoDataFrame):
+            raise Exception(
+                f"Expecting a GeoDataFrame object, but the gdf is {type(gdf)}"
+            )
+
         if len(gdf) == 0:
             logger.warning("No feature found for plotting. Do nothing and return.")
             return ax
 
-        if hasattr(ax, "projection"):
-            gdf = _clean_polygons(data=gdf, projection=ax.projection)
+        if self._plot_engine == PlotEngineType.PYGMT:
+            return plot_geo_data_frame(fig=ax, gdf=gdf, **kwargs)
         else:
-            kwargs["transform"] = self.base_projection
+            if hasattr(ax, "projection"):
+                gdf = _clean_polygons(data=gdf, projection=ax.projection)
+            else:
+                kwargs["transform"] = self.base_projection
 
-        return gdf.plot(ax=ax, **kwargs)
+            return gdf.plot(ax=ax, **kwargs)
 
     @validate_reconstruction_time
     @append_docstring(GET_DATE_DOCSTRING.format("coastlines"))
@@ -1071,10 +1086,10 @@ def get_subduction_direction(self):
         trench_right_features = shapelify_feature_lines(self.trench_right)
 
         gdf_left = gpd.GeoDataFrame(
-            {"geometry": trench_left_features}, geometry="geometry"
+            {"geometry": trench_left_features}, geometry="geometry", crs="EPSG:4326"
         )
         gdf_right = gpd.GeoDataFrame(
-            {"geometry": trench_right_features}, geometry="geometry"
+            {"geometry": trench_right_features}, geometry="geometry", crs="EPSG:4326"
         )
 
         return gdf_left, gdf_right
@@ -1116,6 +1131,13 @@ def plot_subduction_teeth(
             See `Matplotlib` keyword arguments
             [here](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html).
         """
+        if self._plot_engine == PlotEngineType.PYGMT:
+            gdf_subduction_left, gdf_subduction_right = self.get_subduction_direction()
+            pygmt_plot_subduction_zones(
+                ax, gdf_subduction_left, gdf_subduction_right, color=color, **kwargs
+            )
+            return
+
         if not self.plate_reconstruction.topology_features:
             logger.warn(
                 "Plate model does not have topology features. Unable to plot_subduction_teeth."

tests-dir/unittest/test_plot.py

@@ -1,4 +1,7 @@
 #!/usr/bin/env python3
+# import matplotlib
+
+# matplotlib.use("QtAgg")
 
 import sys
 
@@ -11,6 +14,7 @@
 import gplately
 from gplately import PlateReconstruction, PlotTopologies
 
+
 print(gplately.__file__)
 
 # test the plot function with the new PlateModel class