diff --git a/Notebooks/07-WorkingWithPlateTectonicStats.ipynb b/Notebooks/07-WorkingWithPlateTectonicStats.ipynb index 78b8e5b0..4c455be1 100644 --- a/Notebooks/07-WorkingWithPlateTectonicStats.ipynb +++ b/Notebooks/07-WorkingWithPlateTectonicStats.ipynb @@ -17,6 +17,8 @@ "- Mean subduction velocity standard deviation (cm/yr)\n", "- Crustal surface area subducted by trenches over 1 yr (km^2/yr)\n", "- Length of all transform boundaries (km)\n", + "- Mean global trench velocities (cm/yr)\n", + "- Global trench velocity standard deviation (cm/yr)\n", "\n", "for a given plate model (a GPlately `PlateReconstruction` object) and a `reconstruction_time`. \n", "\n", @@ -55,18 +57,11 @@ "metadata": {}, "outputs": [], "source": [ - "pm_manager = PlateModelManager()\n", - "muller2019_model = pm_manager.get_model(\"Muller2019\", data_dir=\"plate-model-repo\")\n", - "\n", - "rotation_model = muller2019_model.get_rotation_model()\n", - "topology_features = muller2019_model.get_topologies()\n", - "static_polygons = muller2019_model.get_static_polygons()\n", - "\n", - "coastlines = muller2019_model.get_layer('Coastlines')\n", - "continents = muller2019_model.get_layer('ContinentalPolygons')\n", - "COBs = muller2019_model.get_layer('COBs')\n", - "\n", - "model = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)" + "muller2019_model = PlateModelManager().get_model(\"Muller2019\", data_dir=\"plate-model-repo\")\n", + "model = gplately.PlateReconstruction(\n", + " muller2019_model.get_rotation_model(), \n", + " muller2019_model.get_topologies(), \n", + " muller2019_model.get_static_polygons())" ] }, { @@ -169,6 +164,15 @@ " \n", " # Global mean subduction velocity and stdev of all trench segments (cm/year)\n", " subd_vel_mean, subd_vel_std = np.mean(subduction_vel), np.std(subduction_vel)\n", + "\n", + " # Absolute velocity(cm/year) of the trench’s motion in the orthogonal direction towards the overriding plate. \n", + " # Negative if moving towards the overriding plate (Trench advance)\n", + " # Positive if moving away from the overriding plate (Trench retreat)\n", + " # The purpose of \"-np.fabs()\" is to allow \"np.cos()\" to produce the correct sign of value(negative or positive).\n", + " trench_absolute_vel = -np.fabs(subduction_data[:,4]) * np.cos(np.radians(subduction_data[:,5]))\n", + "\n", + " # Global mean and standard deviation of trench velocity(cm/year) of all trench segments\n", + " trench_abs_vel_mean, trench_abs_vel_std = np.mean(trench_absolute_vel), np.std(trench_absolute_vel)\n", " \n", " # Area subducted by trenches over 1 yr (km^2/yr)\n", " # Convert subduction velocities from cm/yr to km/yr; trench lengths are already in km.\n", @@ -185,7 +189,8 @@ " # Return a set of boundary stats in a tuple\n", " data = (\n", " reconstruction_time, ridge_len, ridge_vel_mean, ridge_vel_std,\n", - " ridge_surface_area, subduction_len, subd_vel_mean, subd_vel_std, subd_surface_area\n", + " ridge_surface_area, subduction_len, subd_vel_mean, subd_vel_std, subd_surface_area,\n", + " trench_abs_vel_mean, trench_abs_vel_std\n", " )\n", " \n", " return data" @@ -212,7 +217,7 @@ "time_array = np.arange(max_time,min_time-1, -1)\n", "\n", "# Initialise data ndarray\n", - "data_array = np.zeros((9,time_array.size))\n", + "data_array = np.zeros((11,time_array.size))\n", "\n", "# Get plate tectonic stats for each reconstruction time\n", "for t, time in enumerate(time_array):\n", @@ -246,7 +251,9 @@ "'Global subduction zone lengths (km)',\n", "'Mean global subduction velocities (cm/yr)',\n", "'Global subduction velocity standard deviation (cm/yr)',\n", - "'Surface area of crust subducted at trenches (km^2/yr)']\n", + "'Surface area of crust subducted at trenches (km^2/yr)',\n", + "'Mean global trench velocities (cm/yr)',\n", + "'Global trench velocity standard deviation (cm/yr)']\n", "\n", "# Turn data array into a pandas dataframe.\n", "muller19_df = pd.DataFrame(np.column_stack(data_array), columns = headers)\n", @@ -310,10 +317,11 @@ "metadata": {}, "source": [ "### Calculate plate tectonic stats of other plate models\n", - "We repeat the workflow above to calculate the plate tectonic stats of two additional plate reconstruction models:\n", + "We repeat the workflow above to calculate the plate tectonic stats of three additional plate reconstruction models:\n", "\n", "1. Muller et al. 2016\n", - "2. Merdith 2021\n", + "2. Merdith et al. 2021\n", + "3. Zahirovic et al. 2022\n", "\n" ] }, @@ -331,12 +339,15 @@ "muller16_time_array = np.arange(max_time,min_time-1, -1)\n", "\n", "# Initialise data ndarray\n", - "muller16_data = np.zeros((9,muller16_time_array.size))\n", + "muller16_data = np.zeros((11,muller16_time_array.size))\n", "\n", "# Create the PlateReconstruction object with Muller et al 2016\n", - "gdownload = gplately.download.DataServer(\"Muller2016\")\n", - "rotation_model, topology_features, static_polygons = gdownload.get_plate_reconstruction_files()\n", - "muller2016 = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)\n", + "m = PlateModelManager().get_model(\"Muller2016\", data_dir=\"plate-model-repo\")\n", + "muller2016 = gplately.PlateReconstruction(\n", + " m.get_rotation_model(),\n", + " m.get_topologies(),\n", + " m.get_static_polygons(),\n", + ")\n", "\n", "# Get plate tectonic stats for Muller et al 2016\n", "for t, time in enumerate(muller16_time_array):\n", @@ -374,12 +385,15 @@ "merdith21_time_array = np.arange(max_time,min_time-1, -1)\n", "\n", "# Initialise data ndarray\n", - "merdith21_data = np.zeros((9,merdith21_time_array.size))\n", + "merdith21_data = np.zeros((11,merdith21_time_array.size))\n", "\n", "# Create the PlateReconstruction object with Merdith 2021\n", - "gdownload = gplately.download.DataServer(\"Merdith2021\")\n", - "rotation_model, topology_features, static_polygons = gdownload.get_plate_reconstruction_files()\n", - "merdith2021 = gplately.PlateReconstruction(rotation_model, topology_features, static_polygons)\n", + "m = PlateModelManager().get_model(\"Merdith2021\", data_dir=\"plate-model-repo\")\n", + "merdith2021 = gplately.PlateReconstruction(\n", + " m.get_rotation_model(),\n", + " m.get_topologies(),\n", + " m.get_static_polygons(),\n", + ")\n", "\n", "# Get plate tectonic stats for Merdith 2021\n", "for t, time in enumerate(merdith21_time_array):\n", @@ -403,6 +417,52 @@ "merdith21_df.to_csv(\"./NotebookFiles/Merdith2021-PlateTectonicStats.csv\", encoding='utf-8', index=False)" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "82f21113-541b-4993-bee6-2e428f6070b7", + "metadata": {}, + "outputs": [], + "source": [ + "# Zahirovic et al. 2022\n", + "# Create a time array\n", + "min_time = 0 \n", + "max_time = 250\n", + "zahirovic22_time_array = np.arange(max_time,min_time-1, -1)\n", + "\n", + "# Initialise data ndarray\n", + "zahirovic22_data = np.zeros((11, zahirovic22_time_array.size))\n", + "\n", + "# Create the PlateReconstruction object with Zahirovic et al. 2022\n", + "m = PlateModelManager().get_model(\"Zahirovic2022\", data_dir=\"plate-model-repo\")\n", + "zahirovic2022 = gplately.PlateReconstruction(\n", + " m.get_rotation_model(),\n", + " m.get_topologies(),\n", + " m.get_static_polygons(),\n", + ")\n", + "\n", + "# Get plate tectonic stats for Merdith 2021\n", + "for t, time in enumerate(zahirovic22_time_array):\n", + " zahirovic22_data[:,t] = get_plate_tectonic_stats(zahirovic2022, time)\n", + " #print(\"Calculated stats for {} Ma!\".format(time))\n", + " gplately.tools.update_progress((max_time-time)/(max_time-min_time))\n", + "print(\"Done calculating stats!\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1454cb90-79ed-4de4-8024-9ffc0b20424e", + "metadata": {}, + "outputs": [], + "source": [ + "# Turn the Zahirovic et al. 2022 plate tectoncic stats array into a pandas dataframe.\n", + "zahirovic22_df = pd.DataFrame(np.column_stack(zahirovic22_data), columns = headers)\n", + "\n", + "# Save dataframe to CSV\n", + "zahirovic22_df.to_csv(\"./NotebookFiles/Zahirovic2022-PlateTectonicStats.csv\", encoding='utf-8', index=False)" + ] + }, { "cell_type": "markdown", "id": "33a85cf2", @@ -427,19 +487,24 @@ "n_axes = len(axes)\n", "\n", "# Access subduction zone and ridge lengths from pandas dataframes\n", - "muller19_trench_lengths = muller19_df['Global subduction zone lengths (km)'].to_list()\n", - "muller16_trench_lengths = muller16_df['Global subduction zone lengths (km)'].to_list()\n", - "merdith21_trench_lengths = merdith21_df['Global subduction zone lengths (km)'].to_list()\n", - "muller19_ridge_lengths = muller19_df['Global mid-ocean ridge lengths (km)'].to_list()\n", - "muller16_ridge_lengths = muller16_df['Global mid-ocean ridge lengths (km)'].to_list()\n", - "merdith21_ridge_lengths = merdith21_df['Global mid-ocean ridge lengths (km)'].to_list()\n", - "\n", + "column_name = 'Global subduction zone lengths (km)'\n", + "muller19_trench_lengths = muller19_df[column_name].to_list()\n", + "muller16_trench_lengths = muller16_df[column_name].to_list()\n", + "merdith21_trench_lengths = merdith21_df[column_name].to_list()\n", + "zahirovic22_trench_lengths = zahirovic22_df[column_name].to_list()\n", + "\n", + "column_name = 'Global mid-ocean ridge lengths (km)'\n", + "muller19_ridge_lengths = muller19_df[column_name].to_list()\n", + "muller16_ridge_lengths = muller16_df[column_name].to_list()\n", + "merdith21_ridge_lengths = merdith21_df[column_name].to_list()\n", + "zahirovic22_ridge_lengths = zahirovic22_df[column_name].to_list()\n", "\n", "# Prepare different time arrays for each model as Muller 2016 stats start at 230 Ma, while \n", "# the others start at 250 Ma\n", "muller19_time_array = muller19_df['Time (Ma)'].to_list()\n", "muller16_time_array = muller16_df['Time (Ma)'].to_list()\n", - "muller21_time_array = merdith21_df['Time (Ma)'].to_list()\n", + "merdith21_time_array = merdith21_df['Time (Ma)'].to_list()\n", + "zahirovic22_time_array = zahirovic22_df['Time (Ma)'].to_list()\n", "\n", "\n", "# Subduction zone length plot\n", @@ -447,7 +512,8 @@ "axes[0].set_title(\"(a)\", x=-0.1, y=1, fontsize=\"12\")\n", "axes[0].plot(muller19_time_array, gaussian_filter(muller19_trench_lengths, sigma=1), c='k', linestyle=\":\", label=\"Muller et al. 2019\")\n", "axes[0].plot(muller16_time_array, gaussian_filter(muller16_trench_lengths,sigma=1), c='k', alpha = 0.7, label='Muller et al. 2016')\n", - "axes[0].plot(muller21_time_array, gaussian_filter(merdith21_trench_lengths,sigma=1), c='k', alpha = 0.3, label=\"Merdith 2021\")\n", + "axes[0].plot(merdith21_time_array, gaussian_filter(merdith21_trench_lengths,sigma=1), c='k', alpha = 0.3, label=\"Merdith 2021\")\n", + "axes[0].plot(zahirovic22_time_array, gaussian_filter(zahirovic22_trench_lengths,sigma=1), c='k', alpha = 0.7, linestyle=\"-.\", label=\"Zahirovic 2022\")\n", "axes[0].tick_params(direction=\"in\", length=5, top=True, right=True)\n", "axes[0].tick_params(direction=\"in\", which=\"minor\", length=2.5, top=True, right=True)\n", "axes[0].minorticks_on()\n", @@ -457,7 +523,8 @@ "axes[1].set_title(\"(b)\", x=-0.1, y=1, fontsize=\"12\")\n", "axes[1].plot(muller19_time_array, gaussian_filter(muller19_ridge_lengths,sigma=1), c='k', linestyle=\":\", label=\"Muller et al. 2019\")\n", "axes[1].plot(muller16_time_array, gaussian_filter(muller16_ridge_lengths,sigma=1), c='k', alpha = 0.7, label=\"Muller et al. 2016\")\n", - "axes[1].plot(muller21_time_array, gaussian_filter(merdith21_ridge_lengths,sigma=1), c='k', alpha = 0.3, label=\"Merdith 2021\")\n", + "axes[1].plot(merdith21_time_array, gaussian_filter(merdith21_ridge_lengths,sigma=1), c='k', alpha = 0.3, label=\"Merdith 2021\")\n", + "axes[1].plot(zahirovic22_time_array, gaussian_filter(zahirovic22_ridge_lengths,sigma=1), c='k', alpha = 0.7, linestyle=\"-.\", label=\"Zahirovic 2022\")\n", "axes[1].legend(bbox_to_anchor=(0.9, -0.2), ncol = 3, frameon=False)\n", "axes[1].tick_params(direction=\"in\", length=5, top=True, right=True)\n", "axes[1].tick_params(direction=\"in\", which=\"minor\", length=2.5, top=True, right=True)\n", @@ -487,28 +554,50 @@ "metadata": {}, "outputs": [], "source": [ - "fig, axes = plt.subplots(2, 1, figsize=(8,6), dpi=200)\n", + "fig, axes = plt.subplots(3, 1, figsize=(8,9), dpi=200)\n", "fig.subplots_adjust(hspace=0.05)\n", "n_axes = len(axes)\n", "\n", "# Access subduction convergence and ridge spreading rates (along with their standard deviation) from the dataframes\n", "# Turn lists into numpy arrays.\n", - "muller19_subduction_rate = np.array(muller19_df['Mean global ridge velocities (cm/yr)'].to_list())\n", - "muller16_subduction_rate = np.array(muller16_df['Mean global ridge velocities (cm/yr)'].to_list())\n", - "merdith21_subduction_rate = np.array(merdith21_df['Mean global ridge velocities (cm/yr)'].to_list())\n", - "muller19_subduction_rate_sd = np.array(muller19_df['Global ridge velocity standard deviation (cm/yr)'].to_list())\n", - "muller16_subduction_rate_sd = np.array(muller16_df['Global ridge velocity standard deviation (cm/yr)'].to_list())\n", - "merdith21_subduction_rate_sd = np.array(merdith21_df['Global ridge velocity standard deviation (cm/yr)'].to_list())\n", - "\n", - "muller19_spreading_rate = np.array(muller19_df['Mean global subduction velocities (cm/yr)'].to_list())\n", - "muller16_spreading_rate = np.array(muller16_df['Mean global subduction velocities (cm/yr)'].to_list())\n", - "merdith21_spreading_rate = np.array(merdith21_df['Mean global subduction velocities (cm/yr)'].to_list())\n", - "muller19_spreading_rate_sd = np.array(muller19_df['Global subduction velocity standard deviation (cm/yr)'].to_list())\n", - "muller16_spreading_rate_sd = np.array(muller16_df['Global subduction velocity standard deviation (cm/yr)'].to_list())\n", - "merdith21_spreading_rate_sd = np.array(merdith21_df['Global subduction velocity standard deviation (cm/yr)'].to_list())\n", + "column_name = 'Mean global ridge velocities (cm/yr)'\n", + "muller19_subduction_rate = np.array(muller19_df[column_name].to_list())\n", + "muller16_subduction_rate = np.array(muller16_df[column_name].to_list())\n", + "merdith21_subduction_rate = np.array(merdith21_df[column_name].to_list())\n", + "zahirovic22_subduction_rate = np.array(zahirovic22_df[column_name].to_list())\n", + "\n", + "column_name = 'Global ridge velocity standard deviation (cm/yr)'\n", + "muller19_subduction_rate_sd = np.array(muller19_df[column_name].to_list())\n", + "muller16_subduction_rate_sd = np.array(muller16_df[column_name].to_list())\n", + "merdith21_subduction_rate_sd = np.array(merdith21_df[column_name].to_list())\n", + "zahirovic22_subduction_rate_sd = np.array(zahirovic22_df[column_name].to_list())\n", + "\n", + "column_name = 'Mean global subduction velocities (cm/yr)'\n", + "muller19_spreading_rate = np.array(muller19_df[column_name].to_list())\n", + "muller16_spreading_rate = np.array(muller16_df[column_name].to_list())\n", + "merdith21_spreading_rate = np.array(merdith21_df[column_name].to_list())\n", + "zahirovic22_spreading_rate = np.array(zahirovic22_df[column_name].to_list())\n", + "\n", + "column_name = 'Global subduction velocity standard deviation (cm/yr)'\n", + "muller19_spreading_rate_sd = np.array(muller19_df[column_name].to_list())\n", + "muller16_spreading_rate_sd = np.array(muller16_df[column_name].to_list())\n", + "merdith21_spreading_rate_sd = np.array(merdith21_df[column_name].to_list())\n", + "zahirovic22_spreading_rate_sd = np.array(zahirovic22_df[column_name].to_list())\n", + "\n", + "column_name = 'Mean global trench velocities (cm/yr)'\n", + "muller19_trench_rate = np.array(muller19_df[column_name].to_list())\n", + "muller16_trench_rate = np.array(muller16_df[column_name].to_list())\n", + "merdith21_trench_rate = np.array(merdith21_df[column_name].to_list())\n", + "zahirovic22_trench_rate = np.array(zahirovic22_df[column_name].to_list())\n", + "\n", + "column_name = 'Global trench velocity standard deviation (cm/yr)'\n", + "muller19_trench_rate_sd = np.array(muller19_df[column_name].to_list())\n", + "muller16_trench_rate_sd = np.array(muller16_df[column_name].to_list())\n", + "merdith21_trench_rate_sd = np.array(merdith21_df[column_name].to_list())\n", + "zahirovic22_trench_rate_sd = np.array(zahirovic22_df[column_name].to_list())\n", "\n", "# Subduction convergence rates with standard deviation\n", - "axes[0].set_ylabel('Ridge \\n spreading \\n rate (cm/yr)', rotation=\"horizontal\", fontsize=\"12\")\n", + "axes[0].set_ylabel('Ridge spreading \\n rate (cm/yr)', fontsize=\"12\")\n", "axes[0].yaxis.set_label_coords(-0.125, 0.45)\n", "axes[0].set_title(\"(a)\", x=-0.2, y=0.9, fontsize=\"12\")\n", "axes[0].set_ylim([0, 18])\n", @@ -530,6 +619,14 @@ " gaussian_filter(merdith21_subduction_rate-merdith21_subduction_rate_sd, sigma=2), \n", " gaussian_filter(merdith21_subduction_rate+merdith21_subduction_rate_sd, sigma=2),\n", " edgecolor='k', color='k', alpha=0.05)\n", + "\n", + "# Zahirovic 2022\n", + "axes[0].plot(zahirovic22_time_array, gaussian_filter(zahirovic22_subduction_rate,sigma=1), c='k', alpha = 0.7, linestyle=\"-.\", label=\"Zahirovic22 et al. 2022\")\n", + "axes[0].fill_between(merdith21_time_array, \n", + " gaussian_filter(merdith21_subduction_rate-merdith21_subduction_rate_sd, sigma=2), \n", + " gaussian_filter(merdith21_subduction_rate+merdith21_subduction_rate_sd, sigma=2),\n", + " edgecolor='k', color='k', alpha=0.05)\n", + "\n", "# Tick labels\n", "axes[0].tick_params(direction=\"in\", length=5, top=True, right=True)\n", "axes[0].tick_params(direction=\"in\", length=5, top=True, right=True)\n", @@ -538,7 +635,7 @@ "\n", "\n", "# Mid-ocean ridge spreading rates with standard deviation\n", - "axes[1].set_ylabel('Subduction \\n convergence \\n rate (cm/yr)', rotation=\"horizontal\", fontsize=\"12\")\n", + "axes[1].set_ylabel('Subduction convergence \\n rate (cm/yr)', fontsize=\"12\")\n", "axes[1].yaxis.set_label_coords(-0.125, 0.45)\n", "axes[1].set_title(\"(b)\", x=-0.2, y=0.9, fontsize=\"12\")\n", "axes[1].set_ylim([0, 18])\n", @@ -560,12 +657,54 @@ " gaussian_filter(merdith21_spreading_rate+merdith21_spreading_rate_sd, sigma=2),\n", " edgecolor='k', color='k', alpha=0.05)\n", "\n", - "axes[1].legend(bbox_to_anchor=(0.9, -0.25), ncol = 3, frameon=False)\n", + "axes[1].plot(zahirovic22_time_array, gaussian_filter(zahirovic22_spreading_rate,sigma=1), c='k', alpha = 0.5, linestyle=\"-.\", label=\"Zahirovic et al. 2022\")\n", + "axes[1].fill_between(merdith21_time_array, \n", + " gaussian_filter(zahirovic22_spreading_rate - zahirovic22_spreading_rate_sd, sigma=2), \n", + " gaussian_filter(zahirovic22_spreading_rate + zahirovic22_spreading_rate_sd, sigma=2),\n", + " edgecolor='k', color='k', alpha=0.05)\n", + "\n", "axes[1].tick_params(direction=\"in\", length=5, top=True, right=True)\n", "axes[1].tick_params(axis=\"x\", which=\"major\", direction=\"inout\", length=5)\n", "axes[1].tick_params(direction=\"in\", which=\"minor\", length=2.5, top=True, right=True)\n", "axes[1].minorticks_on()\n", "\n", + " \n", + "# trench velocities with standard deviation\n", + "axes[2].set_ylabel('Trench Velocities (cm/yr)', fontsize=\"12\")\n", + "axes[2].yaxis.set_label_coords(-0.125, 0.45)\n", + "axes[2].set_title(\"(c)\", x=-0.2, y=0.9, fontsize=\"12\")\n", + "axes[2].set_ylim([-9.5, 13.5])\n", + "axes[2].plot(muller19_time_array, gaussian_filter(muller19_trench_rate,sigma=1), c='k', linestyle=\":\", label=\"Muller et al. 2019\")\n", + "axes[2].fill_between(muller19_time_array, \n", + " gaussian_filter(muller19_trench_rate - muller19_trench_rate_sd, sigma=2), \n", + " gaussian_filter(muller19_trench_rate + muller19_trench_rate_sd, sigma=2),\n", + " edgecolor='k', color='k', alpha=0.05)\n", + "\n", + "axes[2].plot(muller16_time_array, gaussian_filter(muller16_trench_rate,sigma=1), c='k', alpha = 0.7, label=\"Muller et al. 2016\")\n", + "axes[2].fill_between(muller16_time_array, \n", + " gaussian_filter(muller16_trench_rate - muller16_trench_rate_sd, sigma=2), \n", + " gaussian_filter(muller16_trench_rate + muller16_trench_rate_sd, sigma=2),\n", + " edgecolor='k', color='k', alpha=0.05)\n", + "\n", + "axes[2].plot(merdith21_time_array, gaussian_filter(merdith21_spreading_rate,sigma=1), c='k', alpha = 0.3, label=\"Merdith et al. 2021\")\n", + "axes[2].fill_between(merdith21_time_array, \n", + " gaussian_filter(merdith21_trench_rate - merdith21_trench_rate_sd, sigma=2), \n", + " gaussian_filter(merdith21_trench_rate + merdith21_trench_rate_sd, sigma=2),\n", + " edgecolor='k', color='k', alpha=0.05)\n", + "\n", + "axes[2].plot(zahirovic22_time_array, gaussian_filter(zahirovic22_trench_rate,sigma=1), c='k', alpha = 0.5, linestyle=\"-.\", label=\"Zahirovic et al. 2022\")\n", + "axes[2].fill_between(merdith21_time_array, \n", + " gaussian_filter(zahirovic22_trench_rate - zahirovic22_trench_rate_sd, sigma=2), \n", + " gaussian_filter(zahirovic22_trench_rate + zahirovic22_trench_rate_sd, sigma=2),\n", + " edgecolor='k', color='k', alpha=0.05)\n", + "\n", + "\n", + "axes[2].legend(bbox_to_anchor=(0.9, -0.25), ncol = 3, frameon=False)\n", + "axes[2].tick_params(direction=\"in\", length=5, top=True, right=True)\n", + "axes[2].tick_params(axis=\"x\", which=\"major\", direction=\"inout\", length=5)\n", + "axes[2].tick_params(direction=\"in\", which=\"minor\", length=2.5, top=True, right=True)\n", + "axes[2].minorticks_on()\n", + "\n", "# Place tick labels on the bottom plot only.\n", "for i, ax in enumerate(axes):\n", " ax.set_xlim(max_time, min_time)\n", @@ -597,23 +736,28 @@ "\n", "# Access crustal production and destruction rates from the dataframes\n", "# Turn lists into numpy arrays.\n", - "muller19_crustal_production_rate = np.array(muller19_df['Surface area of crust produced at ridges (km^2/yr)'].to_list())\n", - "muller16_crustal_production_rate = np.array(muller16_df['Surface area of crust produced at ridges (km^2/yr)'].to_list())\n", - "merdith21_crustal_production_rate = np.array(merdith21_df['Surface area of crust produced at ridges (km^2/yr)'].to_list())\n", + "column_name = 'Surface area of crust produced at ridges (km^2/yr)'\n", + "muller19_crustal_production_rate = np.array(muller19_df[column_name].to_list())\n", + "muller16_crustal_production_rate = np.array(muller16_df[column_name].to_list())\n", + "merdith21_crustal_production_rate = np.array(merdith21_df[column_name].to_list())\n", + "zahirovic22_crustal_production_rate = np.array(zahirovic22_df[column_name].to_list())\n", "\n", - "muller19_crustal_destruction_rate = np.array(muller19_df['Surface area of crust subducted at trenches (km^2/yr)'].to_list())\n", - "muller16_crustal_destruction_rate = np.array(muller16_df['Surface area of crust subducted at trenches (km^2/yr)'].to_list())\n", - "merdith21_crustal_destruction_rate = np.array(merdith21_df['Surface area of crust subducted at trenches (km^2/yr)'].to_list())\n", + "column_name = 'Surface area of crust subducted at trenches (km^2/yr)'\n", + "muller19_crustal_destruction_rate = np.array(muller19_df[column_name].to_list())\n", + "muller16_crustal_destruction_rate = np.array(muller16_df[column_name].to_list())\n", + "merdith21_crustal_destruction_rate = np.array(merdith21_df[column_name].to_list())\n", + "zahirovic22_crustal_destruction_rate = np.array(zahirovic22_df[column_name].to_list())\n", "\n", "\n", "# Crustal production rates\n", - "axes[0].set_ylabel('Crustal \\n production \\n rate (km$^2$/yr)', rotation=\"horizontal\", fontsize=\"12\")\n", + "axes[0].set_ylabel('Crustal production \\n rate (km$^2$/yr)', fontsize=\"12\")\n", "axes[0].yaxis.set_label_coords(-0.125, 0.45)\n", "axes[0].set_title(\"(a)\", x=-0.2, y=0.9, fontsize=\"12\")\n", "\n", "axes[0].plot(muller19_time_array, gaussian_filter(muller19_crustal_production_rate, sigma=1), c='k', linestyle=\":\", label=\"Muller et al. 2019\")\n", "axes[0].plot(muller16_time_array, gaussian_filter(muller16_crustal_production_rate,sigma=1), c='k', alpha = 0.7, label='Muller et al. 2016')\n", "axes[0].plot(merdith21_time_array, gaussian_filter(merdith21_crustal_production_rate,sigma=1), c='k', alpha = 0.3, label=\"Merdith et al. 2021\")\n", + "axes[0].plot(zahirovic22_time_array, gaussian_filter(zahirovic22_crustal_production_rate,sigma=1), c='k', linestyle=\"-.\", alpha = 0.5, label=\"Zahirovic et al. 2022\")\n", "\n", "# Tick labels\n", "axes[0].tick_params(direction=\"in\", length=5, top=True, right=True)\n", @@ -623,12 +767,13 @@ "\n", "\n", "# Crustal destruction rates\n", - "axes[1].set_ylabel('Crustal \\n destruction \\n rate (km$^2$/yr)', rotation=\"horizontal\", fontsize=\"12\")\n", + "axes[1].set_ylabel('Crustal destruction \\n rate (km$^2$/yr)', fontsize=\"12\")\n", "axes[1].yaxis.set_label_coords(-0.125, 0.45)\n", "axes[1].set_title(\"(b)\", x=-0.2, y=0.9, fontsize=\"12\")\n", "axes[1].plot(muller19_time_array, gaussian_filter(muller19_crustal_destruction_rate,sigma=1), c='k', linestyle=\":\", label=\"Muller et al. 2019\")\n", "axes[1].plot(muller16_time_array, gaussian_filter(muller16_crustal_destruction_rate,sigma=1), c='k', alpha = 0.7, label=\"Muller et al. 2016\")\n", "axes[1].plot(merdith21_time_array, gaussian_filter(merdith21_crustal_destruction_rate,sigma=1), c='k', alpha = 0.3, label=\"Merdith et al. 2021\")\n", + "axes[1].plot(zahirovic22_time_array, gaussian_filter(zahirovic22_crustal_destruction_rate,sigma=1), c='k', linestyle=\"-.\", alpha = 0.5, label=\"Zahirovic et al. 2022\")\n", "\n", "axes[1].legend(bbox_to_anchor=(0.9, -0.25), ncol = 3, frameon=False)\n", "axes[1].tick_params(direction=\"in\", length=5, top=True, right=True)\n", diff --git a/tests-dir/unittest/trench_velocities.py b/tests-dir/unittest/trench_velocities.py new file mode 100644 index 00000000..c3ecc55b --- /dev/null +++ b/tests-dir/unittest/trench_velocities.py @@ -0,0 +1,55 @@ +import gplately +import numpy as np +from plate_model_manager import PlateModelManager + +zahirovic2022_model = PlateModelManager().get_model( + "Zahirovic2022", data_dir="plate-model-repo" +) +assert zahirovic2022_model + +z22_reconstruction_model = gplately.PlateReconstruction( + zahirovic2022_model.get_rotation_model(), + zahirovic2022_model.get_topologies(), + zahirovic2022_model.get_static_polygons(), +) + +# change this variable for other reconstruction time +reconstruction_time = 0 + +subduction_data = z22_reconstruction_model.tessellate_subduction_zones( + reconstruction_time +) + +# Latitudes and longitudes of points along trench segments +subduction_lon = subduction_data[:, 0] +subduction_lat = subduction_data[:, 1] +print(f"The longitudes of trench sample locations({len(subduction_lon)}):") +print(subduction_lon) +print() +print(f"The latitudes of trench sample locations({len(subduction_lat)}):") +print(subduction_lat) +print() + +# Absolute velocity(cm/year) of the trench’s motion in the orthogonal direction towards the overriding plate. +# Negative if moving towards the overriding plate (Trench advance) +# Positive if moving away from the overriding plate (Trench retreat) +# The purpose of "-np.fabs()" is to allow "np.cos()" to produce the correct sign of value(negative or positive). +trench_absolute_vel = -np.fabs(subduction_data[:, 4]) * np.cos( + np.radians(subduction_data[:, 5]) +) +print( + f"The trench absolute velocities({len(trench_absolute_vel)}) at {reconstruction_time} Ma:" +) +print(trench_absolute_vel) +print() + +# Global mean and standard deviation of trench velocity(cm/year) of all trench segments +trench_abs_vel_mean = np.mean(trench_absolute_vel) +trench_abs_vel_std = np.std(trench_absolute_vel) +print(f"The mean of trench absolute velocities at {reconstruction_time} Ma:") +print(trench_abs_vel_mean) +print() +print( + f"The standard deviation of trench absolute velocities at {reconstruction_time} Ma:" +) +print(trench_abs_vel_std)