Reorganise plotting module

gplately/geometry.py

@@ -79,6 +79,9 @@ class has the following derived GeometryOnSphere classes:
     "PolygonOnSphere",
     "PolylineOnSphere",
     "pygplates_to_shapely",
+    "shapelify_features",
+    "shapelify_feature_lines",
+    "shapelify_feature_polygons",
     "shapely_to_pygplates",
     "wrap_geometries",
 ]
@@ -165,6 +168,113 @@ def to_shapely(self, central_meridian=0.0, tessellate_degrees=None):
         )
 
 
+def shapelify_features(features, central_meridian=0.0, tessellate_degrees=None):
+    """Generate Shapely `MultiPolygon` or `MultiLineString` geometries
+    from reconstructed feature polygons.
+
+    Notes
+    -----
+    Some Shapely polygons generated by `shapelify_features` cut longitudes of 180
+    or -180 degrees. These features may appear unclosed at the dateline, so Shapely
+    "closes" these polygons by connecting any of their open ends with lines. These
+    lines may manifest on GeoAxes plots as horizontal lines that span the entire
+    global extent. To prevent this, `shapelify_features` uses pyGPlates'
+    [DateLineWrapper](https://www.gplates.org/docs/pygplates/generated/pygplates.datelinewrapper)
+    to split a feature polygon into multiple closed polygons if it happens to cut the
+    antimeridian.
+    Another measure taken to ensure features are valid is to order exterior coordinates
+    of Shapely polygons anti-clockwise.
+
+    Parameters
+    ----------
+    features : iterable of <pygplates.Feature>, <ReconstructedFeatureGeometry> or <GeometryOnSphere>
+        Iterable containing reconstructed polygon features.
+    central_meridian : float
+        Central meridian around which to perform wrapping; default: 0.0.
+    tessellate_degrees : float or None
+        If provided, geometries will be tessellated to this resolution prior
+        to wrapping.
+
+    Returns
+    -------
+    all_geometries : list of `shapely.geometry.BaseGeometry`
+        Shapely geometries converted from the given reconstructed features. Any
+        geometries at the dateline are split.
+
+    See Also
+    --------
+    geometry.pygplates_to_shapely : convert PyGPlates geometry objects to
+    Shapely geometries.
+    """
+    if isinstance(
+        features,
+        (
+            pygplates.Feature,
+            pygplates.ReconstructedFeatureGeometry,
+            pygplates.GeometryOnSphere,
+        ),
+    ):
+        features = [features]
+
+    geometries = []
+    for feature in features:
+        if isinstance(feature, pygplates.Feature):
+            geometries.extend(feature.get_all_geometries())
+        elif isinstance(feature, pygplates.ReconstructedFeatureGeometry):
+            geometries.append(feature.get_reconstructed_geometry())
+        elif isinstance(feature, (pygplates.GeometryOnSphere, pygplates.LatLonPoint)):
+            geometries.append(feature)
+        elif isinstance(feature, pygplates.DateLineWrapper.LatLonMultiPoint):
+            geometries.append(
+                pygplates.MultiPointOnSphere(
+                    [i.to_lat_lon() for i in feature.get_points()]
+                )
+            )
+        elif isinstance(feature, pygplates.DateLineWrapper.LatLonPolyline):
+            geometries.append(pygplates.PolylineOnSphere(feature.get_points()))
+        elif isinstance(feature, pygplates.DateLineWrapper.LatLonPolygon):
+            geometries.append(
+                pygplates.PolygonOnSphere(
+                    [i.to_lat_lon() for i in feature.get_exterior_points()]
+                )
+            )
+
+    return [
+        pygplates_to_shapely(
+            i,
+            force_ccw=True,
+            validate=True,
+            central_meridian=central_meridian,
+            tessellate_degrees=tessellate_degrees,
+            explode=False,
+        )
+        for i in geometries
+    ]
+
+
+shapelify_feature_lines = shapelify_features
+shapelify_feature_polygons = shapelify_features
+
+
+def explode_geometries(geometries):
+    if isinstance(geometries, _BaseGeometry):
+        geometries = [geometries]
+    out = []
+    for geometry in geometries:
+        if geometry.is_empty:
+            continue
+        if isinstance(geometry, _BaseMultipartGeometry):
+            out.extend(
+                [
+                    i for i in list(geometry.geoms)
+                    if isinstance(i, _BaseGeometry) and not i.is_empty
+                ]
+            )
+        else:
+            out.append(geometry)
+    return out
+
+
 def pygplates_to_shapely(
     geometry,
     central_meridian=0.0,

gplately/gpml.py

@@ -271,4 +271,13 @@ def _load_FeatureCollection(geometry):
     elif geometry is None:
         return None
     else:
-        raise ValueError("geometry is an invalid type", type(geometry))
+        try:
+            fc = pygplates.FeatureCollection(
+                pygplates.FeaturesFunctionArgument(geometry).get_features()
+            )
+            fc.filenames = []
+            return fc
+        except Exception as err:
+            raise TypeError(
+                "geometry is an invalid type", type(geometry)
+            ) from err

gplately/plot/__init__.py

@@ -0,0 +1,38 @@
+"""Tools for reconstructing and plotting geological features and feature data through time.
+
+Methods in this module reconstruct geological features using
+[pyGPlates' `reconstruct` function](https://www.gplates.org/docs/pygplates/generated/pygplates.reconstruct.html),
+turn them into plottable Shapely geometries, and plot them onto
+Cartopy GeoAxes using Shapely and GeoPandas.
+
+Classes
+-------
+ * `PlotTopologies`
+ * `SubductionTeeth`
+
+Functions
+---------
+ * `plot_subduction_teeth`
+ * `shapelify_features` and its aliases:
+   -  `shapelify_feature_lines`
+   -  `shapelify_feature_polygons`
+"""
+from ..geometry import (
+    shapelify_features,
+    shapelify_feature_lines,
+    shapelify_feature_polygons,
+)
+from .plot_topologies import PlotTopologies
+from .subduction_teeth import (
+    SubductionTeeth,
+    plot_subduction_teeth,
+)
+
+__all__ = [
+    "PlotTopologies",
+    "SubductionTeeth",
+    "plot_subduction_teeth",
+    "shapelify_features",
+    "shapelify_feature_lines",
+    "shapelify_feature_polygons",
+]

gplately/plot/_axes_tools.py

@@ -0,0 +1,34 @@
+import cartopy.crs as ccrs
+import numpy as np
+from matplotlib.axes import Axes
+
+
+def meridian_from_ax(ax: Axes) -> float:
+    if hasattr(ax, "projection") and isinstance(ax.projection, ccrs.Projection):
+        proj = ax.projection
+        return meridian_from_projection(projection=proj)
+    return 0.0
+
+
+def meridian_from_projection(projection: ccrs.Projection) -> float:
+    x = np.mean(projection.x_limits)
+    y = np.mean(projection.y_limits)
+    return ccrs.PlateCarree().transform_point(x, y, projection)[0]
+
+
+def transform_distance_axes(d: float, ax: Axes, inverse=False) -> float:
+    axes_bbox = ax.get_position()
+    fig_bbox = ax.figure.bbox_inches  # display units (inches)
+
+    axes_width = axes_bbox.width * fig_bbox.width
+    axes_height = axes_bbox.height * fig_bbox.height
+
+    # Take mean in case they're different somehow
+    xlim = ax.get_xlim()
+    x_factor = np.abs(xlim[1] - xlim[0]) / axes_width  # map units per display unit
+    ylim = ax.get_ylim()
+    y_factor = np.abs(ylim[1] - ylim[0]) / axes_height
+    factor = 0.5 * (x_factor + y_factor)
+    if inverse:
+        return d / factor
+    return d * factor

gplately/plot/_clean_polygons.py

@@ -0,0 +1,168 @@
+import cartopy.crs as ccrs
+import geopandas as gpd
+import numpy as np
+from shapely.geometry import (
+    LineString,
+    MultiPolygon,
+    Point,
+    Polygon,
+    box,
+)
+from shapely.geometry.base import BaseMultipartGeometry
+from shapely.ops import linemerge, substring
+
+from ._axes_tools import meridian_from_projection
+
+
+def clean_polygons(data, projection):
+    data = gpd.GeoDataFrame(data)
+    data = data.explode(ignore_index=True)
+
+    if data.crs is None:
+        data.crs = ccrs.PlateCarree()
+
+    if isinstance(
+        projection,
+        (
+            ccrs._RectangularProjection,
+            ccrs._WarpedRectangularProjection,
+        ),
+    ):
+        central_longitude = meridian_from_projection(projection)
+        dx = 1.0e-3
+        dy = 5.0e-2
+        rects = (
+            box(
+                central_longitude - 180,
+                -90,
+                central_longitude - 180 + dx,
+                90,
+            ),
+            box(
+                central_longitude + 180 - dx,
+                -90,
+                central_longitude + 180,
+                90,
+            ),
+            box(
+                central_longitude - 180,
+                -90 - dy * 0.5,
+                central_longitude + 180,
+                -90 + dy * 0.5,
+            ),
+            box(
+                central_longitude - 180,
+                90 - dy * 0.5,
+                central_longitude + 180,
+                90 + dy * 0.5,
+            ),
+        )
+        rects = gpd.GeoDataFrame(
+            {"geometry": rects},
+            geometry="geometry",
+            crs=ccrs.PlateCarree(),
+        )
+        data = data.overlay(rects, how="difference")
+
+    projected = data.to_crs(projection)
+
+    # If no [Multi]Polygons, return projected data
+    for geom in projected.geometry:
+        if isinstance(geom, (Polygon, MultiPolygon)):
+            break
+    else:
+        return projected
+
+    proj_width = np.abs(projection.x_limits[1] - projection.x_limits[0])
+    proj_height = np.abs(projection.y_limits[1] - projection.y_limits[0])
+    min_distance = np.mean((proj_width, proj_height)) * 1.0e-4
+
+    boundary = projection.boundary
+    if np.array(boundary.coords).shape[1] == 3:
+        boundary = type(boundary)(np.array(boundary.coords)[:, :2])
+    return fill_all_edges(projected, boundary, min_distance=min_distance)
+
+
+def fill_all_edges(data, boundary, min_distance=None):
+    data = gpd.GeoDataFrame(data).explode(ignore_index=True)
+
+    def drop_func(geom):
+        if hasattr(geom, "exterior"):
+            geom = geom.exterior
+        coords = np.array(geom.coords)
+        return np.all(np.abs(coords) == np.inf)
+
+    to_drop = data.geometry.apply(drop_func)
+    data = (data[~to_drop]).copy()
+
+    def filt_func(geom):
+        if hasattr(geom, "exterior"):
+            geom = geom.exterior
+        coords = np.array(geom.coords)
+        return np.any(np.abs(coords) == np.inf) or (
+            min_distance is not None and geom.distance(boundary) < min_distance
+        )
+
+    to_fix = data.index[data.geometry.apply(filt_func)]
+    for index in to_fix:
+        fixed = fill_edge_polygon(
+            data.geometry.at[index],
+            boundary,
+            min_distance=min_distance,
+        )
+        data.geometry.at[index] = fixed
+    return data
+
+
+def fill_edge_polygon(geometry, boundary, min_distance=None):
+    if isinstance(geometry, BaseMultipartGeometry):
+        return type(geometry)(
+            [fill_edge_polygon(i, boundary, min_distance) for i in geometry.geoms]
+        )
+    if not isinstance(geometry, Polygon):
+        geometry = Polygon(geometry)
+    coords = np.array(geometry.exterior.coords)
+
+    segments_list = []
+    segment = []
+    for x, y in coords:
+        if (np.abs(x) == np.inf or np.abs(y) == np.inf) or (
+            min_distance is not None and boundary.distance(Point(x, y)) <= min_distance
+        ):
+            if len(segment) > 1:
+                segments_list.append(segment)
+                segment = []
+            continue
+        segment.append((x, y))
+    if len(segments_list) == 0:
+        return geometry
+    segments_list = [LineString(i) for i in segments_list]
+
+    out = []
+    for i in range(-1, len(segments_list) - 1):
+        segment_before = segments_list[i]
+        point_before = Point(segment_before.coords[-1])
+
+        segment_after = segments_list[i + 1]
+        point_after = Point(segment_after.coords[0])
+
+        d0 = boundary.project(point_before, normalized=True)
+        d1 = boundary.project(point_after, normalized=True)
+        boundary_segment = substring(boundary, d0, d1, normalized=True)
+
+        if boundary_segment.length > 0.5 * boundary.length:
+            if d1 > d0:
+                seg0 = substring(boundary, d0, 0, normalized=True)
+                seg1 = substring(boundary, 1, d1, normalized=True)
+            else:
+                seg0 = substring(boundary, d0, 1, normalized=True)
+                seg1 = substring(boundary, 0, d1, normalized=True)
+            boundary_segment = linemerge([seg0, seg1])
+
+        if i == -1:
+            out.append(segment_before)
+        out.append(boundary_segment)
+        if i != len(segments_list) - 2:
+            out.append(segment_after)
+
+    return Polygon(np.vstack([i.coords for i in out]))