Calc-Subduction-Stats-Volume-Area-notebook2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=orange>Subduction Zone Kinematics: calculating area and volume flux</font>\n",
    "#### <font color=blue>Notebook 2</font>\n",
    "\n",
    "This notebook illustrates how the subduction zone kinematics genereated in the first notebook ('Make-Subduction-Stats-Table-notebook1.ipynb') can be used to calculate subducting plate area and slab flux. It also illustrates how this data can be visualised and plotted.\n",
    "\n",
    "This notebook therefore uses data precomputed and stored in a csv file, and loaded into a pandas data table for manipulation.\n",
    "\n",
    "\n",
    "This was the second in a series of notebooks used to generate the data presened in the following manuscript:\n",
    "\n",
    "*'Subduction history reveals Cretaceous slab superflux as a possible cause for the mid-Cretaceous plume pulse and superswell events'* by Madison East, R Dietmar Müller, Simon Williams, Sabin Zahirovic and Christian Heine\n",
    "\n",
    "\n",
    "(Notebook written by Simon Williams and modified made by Madison East)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell loads in the required files for the chosen plate model, as well as the precomputed data #\n",
    "#####################################################################################################\n",
    "\n",
    "import pygplates\n",
    "import numpy as np\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import gridspec\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "from topology_plotting import *\n",
    "from matplotlib.patches import Polygon\n",
    "from scipy.special import erfinv\n",
    "import pandas as pd\n",
    "import math\n",
    "import statistics\n",
    "from lithosphere_thickness_plate_model import * # this loads in the coded which is used later to convert seafloor age to lithospheric thickness\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "#######################################################\n",
    "# Define Input files for Muller 2016 AREPS model\n",
    "#######################################################\n",
    "# Specify the name of the GPlates rotation file\n",
    "rotation_filename = '../Data_AREPS_clean/Global_EarthByte_230-0Ma_GK07_AREPS.rot'\n",
    "\n",
    "# input topologies to be used with the subduction_convergence script\n",
    "input_topology_filename = ['../Data_AREPS_clean/Global_EarthByte_230-0Ma_GK07_AREPS_PlateBoundaries.gpml',\\\n",
    "             '../Data_AREPS_clean/Global_EarthByte_230-0Ma_GK07_AREPS_Topology_BuildingBlocks.gpml']\n",
    "\n",
    "# Static polygons to determine whether subduction segment is adjacent to continent or not\n",
    "static_polygon_filename = '../Data_AREPS_clean/Shapefiles/StaticPolygons/Global_EarthByte_GPlates_PresentDay_StaticPlatePolygons_2015_v1.shp'\n",
    "static_polygon_features = pygplates.FeatureCollection(static_polygon_filename)\n",
    "continental_polygon_features = []\n",
    "for feature in static_polygon_features:\n",
    "    if feature.get_feature_type() == pygplates.FeatureType.gpml_closed_continental_boundary:\n",
    "        continental_polygon_features.append(feature)\n",
    "\n",
    "\n",
    "######################################################\n",
    "rotation_model = pygplates.RotationModel(rotation_filename)\n",
    "\n",
    "# specify time range and resolution for plots\n",
    "threshold_sampling_distance_radians = np.radians(0.5)\n",
    "\n",
    "# Define the time snapshots at which to get the subduction zone properties\n",
    "min_time = 0.\n",
    "max_time = 230.\n",
    "time_step = 1.\n",
    "\n",
    "# Set the delta time for velocity calculations\n",
    "velocity_delta_time = 1.\n",
    "\n",
    "# Typically the achor plate id should be 0\n",
    "anchor_plate_id = 0\n",
    "\n",
    "############################################################################################################\n",
    "# Previously calculated results are loaded here from a .csv file, and stored in a pandas data table.\n",
    "# This is the data that was generated in the first notebook ('Make-Subduction-Stats-Table-notebook1.ipynb')\n",
    "df = pd.read_csv('SubductionTable_clean_0_230Ma.csv')\n",
    "\n",
    "# This takes care of any minor construction errors in the plate model where small sections of convergence zones may momentarily become divergent\n",
    "# this converts any negative convergence rate to zero (these were very minor occurances in the AREPS model, but if using a new model, check these instances first)\n",
    "df.ix[df.conv_rate<0,'conv_rate'] = 0\n",
    "\n",
    "############################################################################################################\n",
    "# Below we calculate the orthogonal migration and convergence rates, which are important to get true estimates of area and volume flux\n",
    "\n",
    "mr = np.asarray(df['migr_rate'])\n",
    "mo = np.asarray(df['migr_obliq'])\n",
    "df['ortho_migr_rate'] = pd.Series(mr*np.sin(np.radians(np.abs(mo))), index=df.index)\n",
    "\n",
    "cr = np.asarray(df['conv_rate'])\n",
    "co = np.asarray(df['conv_obliq'])\n",
    "df['ortho_conv_rate'] = pd.Series(cr*np.abs(np.cos(np.radians(co))), index=df.index)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time for the calculations..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "########### Lithospheric Thickness and Subduction Volume Calculations  ###############\n",
    "\n",
    "## Thickness (in metres): Equation and constants from 'The Solid Earth' (2005), Sandwell (2001) and Grose et al. (2012) \n",
    "## T-values in degrees celcius, thermal diffusivity (k) in m^2/sec\n",
    "T1 = 1150.\n",
    "To = 0.\n",
    "Tm = 1350.\n",
    "kappa = 0.804e-6\n",
    "Myr2sec=1e6*365*24*60*60\n",
    "lithosphere_thickness_boundarylayer = erfinv((T1-To)/(Tm-To))*2*np.sqrt(kappa)*np.sqrt(Myr2sec*np.asarray(df['SeafloorAge']))\n",
    "\n",
    "for i in range(0, len(lithosphere_thickness_boundarylayer)):\n",
    "    if (lithosphere_thickness_boundarylayer[i] > 125000) :\n",
    "        lithosphere_thickness_boundarylayer[i] = 125000\n",
    "        \n",
    "## Thickness using translated MATLAB script\n",
    "lithosphere_thickness_plate_model = list(lithosphere_thickness_boundarylayer)\n",
    "for i in range(0, len(lithosphere_thickness_boundarylayer)) :\n",
    "#for i in range(0, 100) :\n",
    "    lithosphere_thickness_plate_model[i] = plate_isotherm_depth(np.asarray(df['SeafloorAge'])[i], T1)\n",
    "\n",
    "## To convert arc_length from degrees on a sphere to m (using earth's radius = 6371000 m)\n",
    "arc_length_m = 2*math.pi*6371000*((np.asarray(df['arc_length']))/360)\n",
    "\n",
    "## Calculate Subduction Volume (in m^3 per year)\n",
    "subduction_volume_m3y = (np.asarray(df['ortho_conv_rate'])/100) * lithosphere_thickness_plate_model * arc_length_m\n",
    "\n",
    "## Calculate Subduction Volume (in m^3 per Ma)\n",
    "subduction_volume_m3Ma = (np.asarray(df['ortho_conv_rate'])/100) * lithosphere_thickness_plate_model * arc_length_m * 1e6\n",
    "\n",
    "## Calculate Subduciton Volume (slab flux) (in km^3 per year)\n",
    "subduction_volume_km3y = subduction_volume_m3y/1e9 \n",
    "\n",
    "print len(df)\n",
    "print np.asarray(df['SeafloorAge'])\n",
    "print lithosphere_thickness_boundarylayer\n",
    "print arc_length_m\n",
    "print subduction_volume_m3y\n",
    "#print lithosphere_thickness_plate_model\n",
    "#print pd.Series(lithosphere_thickness)\n",
    "\n",
    "# Add the parameters calculated above to our data table\n",
    "df['arc_length_m'] = pd.Series(arc_length_m, index=df.index)\n",
    "df['lithosphere_thickness_boundarylayer'] = pd.Series(lithosphere_thickness_boundarylayer, index=df.index)\n",
    "df['lithosphere_thickness_plate_model'] = pd.Series(lithosphere_thickness_plate_model, index=df.index)\n",
    "df['subduction_volume_m3y'] = pd.Series(subduction_volume_m3y, index=df.index)\n",
    "df['subduction_volume_m3Ma'] = pd.Series(subduction_volume_m3Ma, index=df.index)\n",
    "df['subduction_volume_km3y'] = pd.Series(subduction_volume_km3y, index=df.index)\n",
    "\n",
    "\n",
    "########## Area calculations ##############\n",
    "\n",
    "#To calculate subducting plate area (in m^2/year)\n",
    "subduction_area_m2y = (np.asarray(df['ortho_conv_rate'])/100) * arc_length_m \n",
    "\n",
    "# To calculate subducitng plate area (in km^2/year)\n",
    "subduction_area_km2y = subduction_area_m2y/1e6\n",
    "\n",
    "df['subduction_area_m2y'] = pd.Series(subduction_area_m2y, index=df.index)\n",
    "df['subduction_area_km2y'] = pd.Series(subduction_area_km2y, index=df.index)\n",
    "\n",
    "\n",
    "df_AllTimes = df\n",
    "\n",
    "# Display the contents of the table containing all subduction\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save data to csv file!!\n",
    "df_AllTimes.to_csv('SubductionVolumesAreaTable_clean_%0.0f_%0.0fMa.csv' % (min_time,max_time))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualisation...\n",
    "\n",
    "The plot below shows time snapshots of the model plate boundaries and allows you to visualise your calcualted data\n",
    "- each subduction zone segment is represented by a coloured dot, and according to the colour map, depicting what you have chosen to plot (e.g. slab flux or subducting plate area)\n",
    "- bright blue lines show mid-ocean ridges\n",
    "- orange lines show 'other' plate boundaries\n",
    "\n",
    "You can choose what range of time steps you want to plot and what calculation to show. Plot features can be altered and the maps can be exported. By exporting all time steps, you have the potential to generate an animation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.patches as mpatches\n",
    "from topology_plotting import *\n",
    "\n",
    "#Define which timesteps you want to plot\n",
    "times = np.arange(225,231,1)\n",
    "for time in times:\n",
    "    subset = df_AllTimes[(df_AllTimes['time']==time)]\n",
    "\n",
    "    coastlines_file = '../Data_most_recent/Global_EarthByte_230-0Ma_GK07_AREPS_Coastlines.gpmlz'\n",
    "    pygplates.reconstruct(coastlines_file, rotation_model, 'tmp.shp', time)\n",
    "\n",
    "    pygplates.reconstruct(continental_polygon_features, rotation_model, 'tmp2.shp', time)\n",
    "\n",
    "    fig = plt.figure(figsize=(16,6),dpi=150)\n",
    "    ax_map = fig.add_axes([0,0,0.8,1.0])\n",
    "    lon0=0\n",
    "    m = Basemap(projection='moll', lon_0=lon0, resolution='c', ax=ax_map)\n",
    "    #m = Basemap(resolution='c',projection='ortho',lat_0=30.,lon_0=lon0)\n",
    "    cp = m.drawmapboundary()\n",
    "    m.drawparallels(np.arange(-90,90,30))\n",
    "    m.drawmeridians(np.arange(-180,180,30))\n",
    "\n",
    "    # Plot reconstructed coastlines\n",
    "    shp_info = m.readshapefile('tmp','shp',drawbounds=True,color='none')\n",
    "    for nshape,seg in enumerate(m.shp):\n",
    "        poly = Polygon(seg,facecolor='grey',edgecolor='none',alpha=0.8,zorder=0.5)\n",
    "        plt.gca().add_patch(poly)\n",
    "\n",
    "    # Plot reconstructed continental shelf\n",
    "    shp_info = m.readshapefile('tmp2','shp',drawbounds=True,color='none')\n",
    "    for nshape,seg in enumerate(m.shp):\n",
    "        poly = Polygon(seg,facecolor='grey',edgecolor='none',alpha=0.3,zorder=0.25)\n",
    "        plt.gca().add_patch(poly)\n",
    "\n",
    "    plot_velocities_and_topologies(m,input_topology_filename,rotation_model,time,\n",
    "                                   delta_time=5,res=10,scale=4000,lon0=lon0,clip_path=cp)\n",
    "\n",
    "    # in the below lines you define which of your calculations you want to plot e.g. subduction_volume_km3y\n",
    "    x, y = m(np.asarray(subset.lon), np.asarray(subset.lat))\n",
    "    l1 = m.scatter(x,y,c=subset['subduction_volume_km3y'],s=40,edgecolor='',zorder=5,\n",
    "                   cmap=plt.cm.Blues,vmin=0,vmax=6.5e-1)\n",
    "\n",
    "    #plt.title('Subduction Volume Flux (%s Ma)' % time, size=18)\n",
    "    plt.title('%s Ma' % time, size=38)\n",
    "\n",
    "    plt.colorbar(l1,extend='max',label='Volume ($km^3$/yr)')\n",
    "    \n",
    "    #velocity = mpatches.Patch(color='w', label='Velocity Vectors',alpha=0.4)\n",
    "    #MOR = mpatches.Patch(color='darkturquoise', label='Mid Ocean Ridge',alpha=0.5)\n",
    "    #transform = mpatches.Patch(color='darkorange', label='Fracture Zones',alpha=0.6)\n",
    "    #continents = mpatches.Patch(color='grey', label='Continental Shelf',alpha=0.3)\n",
    "    #coastlines = mpatches.Patch(color='grey', label='Continents',alpha=0.7)\n",
    "    \n",
    "    #plt.legend(handles=[velocity, MOR, transform, continents, coastlines], loc=1)\n",
    "    \n",
    "\n",
    "    #plt.savefig('Volume_flux_animation/figures_vol/Subduction_Volume_Flux_%s Ma.jpeg' % time ,bbox_inches='tight', dpi=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Total Volume (slab flux)\n",
    "To understand temporal changes, rather than spatial variation at a single time step, we will need to combine the data from our numerous subduction segments at each time step. The below cell sums the total subduction volume for each 1 Ma timestep and saves it as a csv file. A separate file is made for volume per year and per Ma."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Saving Summed Values for 'm^3 per year'\n",
    "\n",
    "times = np.arange(0,231,1)\n",
    "total_sub_volume_m3y = []\n",
    "for time in times:\n",
    "    subset1 = df_AllTimes[(df_AllTimes['time']==time)]\n",
    "    total_sub_volume_m3y.append(math.fsum(subset1.subduction_volume_m3y))    \n",
    "#print total_sub_volume_m3y\n",
    "\n",
    "data = pd.DataFrame()\n",
    "data['Time_(Ma)'] = pd.Series(times)\n",
    "data['Subduction_Volume_Flux_(m^3/year)'] = pd.Series(total_sub_volume_m3y)\n",
    "data\n",
    "data.to_csv('Summed_data_plus_stats/Subduction_Volume_Flux_Summed_Data_clean_m3y_%0.0f_%0.0fMa.csv' % (min_time,max_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Saving Summed Values for 'm^3 per Ma'\n",
    "\n",
    "times = np.arange(0,231,1)\n",
    "total_sub_volume_m3Ma = []\n",
    "for time in times:\n",
    "    subset1 = df_AllTimes[(df_AllTimes['time']==time)]\n",
    "    total_sub_volume_m3Ma.append(math.fsum(subset1.subduction_volume_m3Ma))    \n",
    "#print total_sub_volume_m3Ma\n",
    "\n",
    "data = pd.DataFrame()\n",
    "data['Time_(Ma)'] = pd.Series(times)\n",
    "data['Subduction_Volume_Flux_(m^3/Ma)'] = pd.Series(total_sub_volume_m3Ma)\n",
    "data\n",
    "data.to_csv('Summed_data_plus_stats/Subduction_Volume_Flux_Summed_Data_clean_m3Ma_%0.0f_%0.0fMa.csv' % (min_time,max_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Saving Summed Values for 'km^3 per yr'\n",
    "\n",
    "times = np.arange(0,231,1)\n",
    "total_sub_volume_km3y = []\n",
    "for time in times:\n",
    "    subset1 = df_AllTimes[(df_AllTimes['time']==time)]\n",
    "    total_sub_volume_km3y.append(math.fsum(subset1.subduction_volume_km3y))    \n",
    "#print total_sub_volume_km3y\n",
    "\n",
    "data = pd.DataFrame()\n",
    "data['Time_(Ma)'] = pd.Series(times)\n",
    "data['Subduction_Volume_Flux_(km^3/year)'] = pd.Series(total_sub_volume_km3y)\n",
    "data\n",
    "data.to_csv('Summed_data_plus_stats/Subduction_Volume_Flux_Summed_Data_clean_km3y_%0.0f_%0.0fMa.csv' % (min_time,max_time))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Total Area (subducting plate area)\n",
    "Here we do the same thing that was done for volume, but for area values instead."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Summing area values for subducting plate area (in m^2/yr)\n",
    "\n",
    "times = np.arange(0,231,1)\n",
    "total_sub_area_m2y = []\n",
    "for time in times:\n",
    "    subset1 = df_AllTimes[(df_AllTimes['time']==time)]\n",
    "    total_sub_area_m2y.append(math.fsum(subset1.subduction_area_m2y))\n",
    "#print total_sub_area_m2y \n",
    "\n",
    "data = pd.DataFrame()\n",
    "data['Time_(Ma)'] = pd.Series(times)\n",
    "data['Subduction_Area_Flux_(m^2/year)'] = pd.Series(total_sub_area_m2y)\n",
    "data\n",
    "data.to_csv('Summed_data_plus_stats/Subduction_Area_Flux_Summed_Data_clean_m2y_%0.0f_%0.0fMa.csv' % (min_time,max_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Summing area values for subducting plate area (in km^2/yr)\n",
    "\n",
    "times = np.arange(0,231,1)\n",
    "total_sub_area_km2y = []\n",
    "for time in times:\n",
    "    subset1 = df_AllTimes[(df_AllTimes['time']==time)]\n",
    "    total_sub_area_km2y.append(math.fsum(subset1.subduction_area_km2y))\n",
    "#print total_sub_area_km2y \n",
    "\n",
    "data = pd.DataFrame()\n",
    "data['Time_(Ma)'] = pd.Series(times)\n",
    "data['Subduction_Area_Flux_(km^2/year)'] = pd.Series(total_sub_area_km2y)\n",
    "data\n",
    "data.to_csv('Summed_data_plus_stats/Subduction_Area_Flux_Summed_Data_clean_km2y_%0.0f_%0.0fMa.csv' % (min_time,max_time))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}