From 1df3c448106780c0735a625f6291069873dd4aea Mon Sep 17 00:00:00 2001 From: John Cannon Date: Fri, 25 Oct 2024 15:27:04 +1100 Subject: [PATCH] 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 --- gplately/ptt/subduction_convergence.py | 5 +- tests-dir/debug_utils/github_issue_269.py | 140 ++++++++++++++++++ .../pytestcases/test_1_reconstructions.py | 9 +- tests-dir/pytestcases/test_2_points.py | 23 ++- 4 files changed, 164 insertions(+), 13 deletions(-) create mode 100644 tests-dir/debug_utils/github_issue_269.py diff --git a/gplately/ptt/subduction_convergence.py b/gplately/ptt/subduction_convergence.py index b6806fc3..c5853a4a 100644 --- a/gplately/ptt/subduction_convergence.py +++ b/gplately/ptt/subduction_convergence.py @@ -422,7 +422,6 @@ 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(): @@ -430,8 +429,6 @@ def subduction_convergence( 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, @@ -439,7 +436,7 @@ def subduction_convergence( 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, diff --git a/tests-dir/debug_utils/github_issue_269.py b/tests-dir/debug_utils/github_issue_269.py new file mode 100644 index 00000000..300df340 --- /dev/null +++ b/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() diff --git a/tests-dir/pytestcases/test_1_reconstructions.py b/tests-dir/pytestcases/test_1_reconstructions.py index 160a6502..bd257298 100644 --- a/tests-dir/pytestcases/test_1_reconstructions.py +++ b/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): diff --git a/tests-dir/pytestcases/test_2_points.py b/tests-dir/pytestcases/test_2_points.py index d5ae3641..ed0f550a 100644 --- a/tests-dir/pytestcases/test_2_points.py +++ b/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 + )