269 incorrect trench normal angle calculation (#296)

* add debug file for issue 269

* Fix incorrect double reversal of topological line sub-segment normals.

The topo line's sub-segment geometries are already getting reversed
(when necessary) which already reverses its great circle arc normals.

* increase the error magin from 1000 to 3000 in test_cont_arc_length() function

* fix test_2_points.py. assert the two numpy arrays are not emtpy.

* change error margin back to 1000 and change the approx_lengths at 100ma from 44000 to 46000 in test_cont_arc_length() function

---------

Co-authored-by: michaelchin <michaelchin78@hotmail.com>

gplately/ptt/subduction_convergence.py

@@ -422,24 +422,21 @@ def subduction_convergence(
                     )
 
                     sub_sub_segment_geometry = sub_sub_segment.get_resolved_geometry()
-                    sub_sub_segment_trench_normal_reversal = trench_normal_reversal
                     # If sub-sub-segment was reversed when it contributed to the topological line shared sub-segment then
                     # we need to use that reversed geometry so that it has the same order of points as the topological line.
                     if sub_sub_segment.was_geometry_reversed_in_topology():
                         # Create a new sub-sub-segment polyline with points in reverse order.
                         sub_sub_segment_geometry = pygplates.PolylineOnSphere(
                             sub_sub_segment_geometry[::-1]
                         )
-                        #  The trench normal direction is also reversed.
-                        sub_sub_segment_trench_normal_reversal = -trench_normal_reversal
 
                     _sub_segment_subduction_convergence(
                         output_data,
                         time,
                         sub_sub_segment_geometry,
                         trench_plate_id,
                         subducting_plate_id,
-                        sub_sub_segment_trench_normal_reversal,
+                        trench_normal_reversal,
                         trench_length_radians,
                         distance_along_trench_radians,
                         threshold_sampling_distance_radians,

tests-dir/debug_utils/github_issue_269.py

@@ -0,0 +1,140 @@
+import cartopy.crs as ccrs
+import cartopy.mpl.gridliner as grd
+import geopandas as gpd
+import matplotlib.cm as cm
+import matplotlib.pyplot as plt
+import numpy as np
+from plate_model_manager import PlateModelManager
+from shapely.geometry import LineString, Point
+
+import gplately
+
+pm_manager = PlateModelManager()
+muller2019_model = pm_manager.get_model("Muller2019", data_dir="plate-model-repo")
+rotation_model = muller2019_model.get_rotation_model()
+topology_features = muller2019_model.get_topologies()
+static_polygons = muller2019_model.get_static_polygons()
+coastlines = muller2019_model.get_coastlines()
+continents = muller2019_model.get_continental_polygons()
+COBs = muller2019_model.get_COBs()
+
+## loading reconstruction
+model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)
+time = 100  # Ma
+subduction_data = model.tessellate_subduction_zones(time, return_geodataframe=True)
+
+
+profiles = []
+trench_lats = []
+trench_lons = []
+n_steps = 20
+for i in range(0, len(subduction_data) - 1):
+    dlon1 = n_steps * np.sin(
+        np.radians(subduction_data.iloc[i]["trench normal angle (degrees)"])
+    )
+    dlat1 = n_steps * np.cos(
+        np.radians(subduction_data.iloc[i]["trench normal angle (degrees)"])
+    )
+    x1 = subduction_data.iloc[i]["geometry"].x
+    y1 = subduction_data.iloc[i]["geometry"].y
+    ilon1 = x1 + dlon1
+    ilat1 = y1 + dlat1
+    start_point = Point(x1 + 0.1 * dlon1, y1 + 0.1 * dlat1)
+    end_point = Point(ilon1, ilat1)
+    profile = LineString([start_point, end_point])
+    profiles.append(profile)
+
+central_longitude = 180
+Profile = gpd.GeoDataFrame(geometry=profiles)
+mollweide_proj = f"+proj=moll +lon_0={central_longitude} +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
+
+Profile = Profile.set_crs(mollweide_proj)
+
+
+cmap = "tab20c"
+fig, axes = plt.subplots(
+    2,
+    1,
+    figsize=(16, 12),
+    dpi=100,
+    subplot_kw={"projection": ccrs.Mollweide(central_longitude=central_longitude)},
+)
+
+# First subplot (left side)
+ax1 = axes[0]
+ax1.gridlines(
+    color="0.7",
+    linestyle="--",
+    xlocs=np.arange(-180, 180, 15),
+    ylocs=np.arange(-90, 90, 15),
+)
+
+gplot1 = gplately.PlotTopologies(
+    model, coastlines=coastlines, continents=continents, COBs=COBs, time=time
+)
+# ax1.set_title('Subduction zones and mid-ocean ridges reconstructed to %i Ma' % (time))
+
+# Plot shapefile features, subduction zones, and MOR boundaries at the specified time
+gplot1.plot_continent_ocean_boundaries(ax1, color="b", alpha=0.05)
+gplot1.plot_continents(ax1, facecolor="palegoldenrod", alpha=0.2)
+gplot1.plot_coastlines(ax1, color="DarkKhaki")
+# gplot1.plot_ridges_and_transforms(ax1, color="red")
+gplot1.plot_trenches(ax1, color="k", alpha=0.3)
+gplot1.plot_subduction_teeth(ax1, color="k", alpha=0.7)
+
+# Add subduction data plot (with a colorbar)
+cbar1 = subduction_data.plot(
+    ax=ax1,
+    column="trench normal angle (degrees)",
+    transform=ccrs.PlateCarree(),
+    alpha=0.2,
+    cmap=cmap,
+    legend=False,
+)
+
+sm = cm.ScalarMappable(cmap=cmap)
+sm.set_array(subduction_data["trench normal angle (degrees)"])
+sm.set_clim(0, 360)
+
+# Add colorbar using the same Axes object used for plotting
+colorbar = plt.colorbar(
+    sm,
+    ax=ax1,
+    orientation="vertical",
+    shrink=0.6,
+    label="trench normal angle (degrees)",
+)
+
+# Second subplot (right side)
+ax2 = axes[1]
+ax2.gridlines(
+    color="0.7",
+    linestyle="--",
+    xlocs=np.arange(-180, 180, 15),
+    ylocs=np.arange(-90, 90, 15),
+)
+
+gplot2 = gplately.PlotTopologies(
+    model, coastlines=coastlines, continents=continents, COBs=COBs, time=time
+)
+ax2.set_title("Profile View")
+
+# Plot shapefile features
+gplot2.plot_continent_ocean_boundaries(ax2, color="b", alpha=0.05)
+gplot2.plot_continents(ax2, facecolor="palegoldenrod", alpha=0.2)
+gplot2.plot_coastlines(ax2, color="DarkKhaki")
+# gplot2.plot_ridges_and_transforms(ax2, color="red")
+gplot2.plot_trenches(ax2, color="k", alpha=0.2)
+gplot2.plot_subduction_teeth(ax2, color="k")
+
+# Add profile plot (with colorbar)
+cbar2 = Profile.plot(ax=ax2, transform=ccrs.PlateCarree(), alpha=0.2)
+# fig.colorbar(cbar2, ax=ax2, fraction=0.04, pad=0.04)
+
+# Set the global extents for both subplots
+ax1.set_global()
+ax2.set_global()
+
+# Adjust layout
+plt.tight_layout()
+plt.show()

tests-dir/pytestcases/test_1_reconstructions.py

@@ -153,7 +153,7 @@ def test_cont_arc_length(time, model, muller_2019_model):
     )
     approx_lengths = {
         0: 52200,
-        100: 44000,
+        100: 46000,
     }
     err_msg = (
         "Could not calculate total continental arc lengths for Muller et "
@@ -198,10 +198,9 @@ def test_rotation_model_copy(model):
 
 
 def test_reconstruction_filenames(model):
-    assert (
-        "Muller_etal_2019_PlateBoundaries_DeformingNetworks.gpmlz"
-        in [os.path.basename(i) for i in model.topology_features.filenames]
-    )
+    assert "Muller_etal_2019_PlateBoundaries_DeformingNetworks.gpmlz" in [
+        os.path.basename(i) for i in model.topology_features.filenames
+    ]
 
 
 def test_pickle_reconstruction(model):

tests-dir/pytestcases/test_2_points.py

@@ -24,25 +24,37 @@
 
 """
 
+
 # TESTING THE POINTS OBJECT
 @pytest.mark.parametrize("time", reconstruction_times)
 def test_point_reconstruction(time, gpts):
     rlons, rlats = gpts.reconstruct(time, return_array=True, anchor_plate_id=0)
-    assert (rlons, rlats), "Unable to reconstruct point data to {} Ma with Muller et al. (2019).".format(time)
+    assert (
+        rlons.size and rlats.size
+    ), "Unable to reconstruct point data to {} Ma with Muller et al. (2019).".format(
+        time
+    )
+
 
-
 # TESTING PLATE VELOCITY CALCULATIONS
 @pytest.mark.parametrize("time", reconstruction_times)
 def test_plate_velocity(time, gpts):
     plate_vel = gpts.plate_velocity(time, delta_time=1)
-    assert plate_vel, "Unable to calculate plate velocities of point data at {} Ma with Muller et al. (2019).".format(time)
+    assert (
+        plate_vel
+    ), "Unable to calculate plate velocities of point data at {} Ma with Muller et al. (2019).".format(
+        time
+    )
+
 
 def test_point_attributes(gpts):
     attr = np.arange(0, gpts.size)
     gpts.add_attributes(FROMAGE=attr, TOAGE=attr)
 
+
 def test_pickle_Points(gpts):
     import pickle
+
     gpts_dump = pickle.dumps(gpts)
     gpts_load = pickle.loads(gpts_dump)
 
@@ -51,5 +63,8 @@ def test_pickle_Points(gpts):
     gpts_dump = pickle.dumps(gpts)
     gpts_load = pickle.loads(gpts_dump)
 
+
 def test_change_ancbor_plate(gpts):
-    gpts.rotate_reference_frames(50, from_rotation_reference_plate=0, to_rotation_reference_plate=101)
+    gpts.rotate_reference_frames(
+        50, from_rotation_reference_plate=0, to_rotation_reference_plate=101
+    )