Make-Subduction-Stats-Table-notebook1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color=orange>Extracting Subduction Zone Kinematics & Parameters (pyGPlates method)</font>\n",
    "\n",
    "#### <font color=blue>Notebook 1</font>\n",
    "\n",
    "**This notebook can be used to extract subduction zone kinematics and parameters from a chosen plate model. These include:**\n",
    "- overriding and subducting plate ID's\n",
    "- arc length\n",
    "- lat and lon\n",
    "- time\n",
    "- convergence rate and obliquity\n",
    "- migration rate and obliquity \n",
    "- age of the subducting seafloor\n",
    "- distance to nearest continent\n",
    "\n",
    "**It illustrates usage of the updated subduction zone kinematics pygplates method, using the Müller et al AREPS model:**\n",
    "\n",
    "- implemented to work directly on the plate model gpml files (whereas before, you had to export the resolved topologies at 1Myr increments first)\n",
    "- this notebook also interpolates the age of the subducting plate age from age grids if provided\n",
    "- within the notebook data are stored in pandas dataframes, and there is an illustration of how to export them to both csv or hdf5 files\n",
    "\n",
    "This was the first 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 with minor modifications made by Madison East)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell loads in the required files for the chosen plate model and defines some required functions #\n",
    "#########################################################################################################\n",
    "\n",
    "import pygplates\n",
    "import subduction_convergence as sc #(required code stored in a .py file)\n",
    "import numpy as np\n",
    "import os\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import gridspec\n",
    "from netCDF4 import Dataset\n",
    "import scipy.interpolate as spi\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "import inpaint\n",
    "import pandas as pd\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "###################################################################\n",
    "# Define Input files for Muller 2016 AREPS model (or chosen 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",
    "# The list array here specifies the age grid files, \n",
    "# assuming that the age grids are all together in a folder, with common filenames that\n",
    "# differ only in the integer number somewhere in the filename that corresponds to the age in Ma\n",
    "# The two parts of the list specify the parts of the filename before and after the number\n",
    "# For example:\n",
    "# ['/path_to_age_grids/age_grids_at_time_','_Ma.grd]\n",
    "agegrid_filename = ['../AgeGrids_most_recent/agegrid_final_nomask_','.grd']\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",
    "#################################################\n",
    "# Included functions\n",
    "#################################################\n",
    "def sample_grid_using_scipy(x,y,grdfile):\n",
    "    \n",
    "    with Dataset(grdfile, 'r', format='NETCDF4') as data:\n",
    "        try:\n",
    "            lon = np.asarray(data['x'])\n",
    "            lat = np.asarray(data['y'])\n",
    "        except KeyError:\n",
    "            lon = np.asarray(data['lon'])\n",
    "            lat = np.asarray(data['lat'])\n",
    "    \n",
    "        Zg = np.asarray(data['z'])\n",
    "    \n",
    "    test = inpaint.fill_ndimage(Zg)\n",
    "    \n",
    "    lut=spi.RectBivariateSpline(lon,lat,test.T)\n",
    "    result = []\n",
    "    for xi,yi in zip(x,y):\n",
    "        result.append(lut(xi, yi)[0][0])\n",
    "            \n",
    "    return result\n",
    "\n",
    "\n",
    "def get_nearest_continental_polygon(subduction_data,continental_polygons):\n",
    "\n",
    "    nearest_continent_plate_id = []\n",
    "    distance_to_nearest_continent = []\n",
    "    \n",
    "    for index,row in subduction_data.iterrows():\n",
    "        \n",
    "        min_distance_to_all_features = np.radians(180)\n",
    "        nearest_continent = None        \n",
    "        \n",
    "        seed_point = pygplates.PointOnSphere(row['lat'], row['lon'])\n",
    "        \n",
    "        for polygon in continental_polygons:\n",
    "            if polygon is not None:\n",
    "                min_distance_to_feature = pygplates.GeometryOnSphere.distance(\n",
    "                    polygon.get_reconstructed_geometry(),\n",
    "                    seed_point,\n",
    "                    min_distance_to_all_features,\n",
    "                    geometry1_is_solid=True)\n",
    "\n",
    "                # If the current geometry is nearer than all previous geometries then\n",
    "                # its associated feature is the nearest feature so far.\n",
    "                if min_distance_to_feature is not None:\n",
    "                    min_distance_to_all_features = min_distance_to_feature\n",
    "                    nearest_continent = polygon.get_feature().get_reconstruction_plate_id()\n",
    "                    \n",
    "        nearest_continent_plate_id.append(nearest_continent)\n",
    "        distance_to_nearest_continent.append(min_distance_to_all_features*pygplates.Earth.mean_radius_in_kms)\n",
    "\n",
    "    return nearest_continent_plate_id,distance_to_nearest_continent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Number crunching section"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a dataframe that contains all subduction data for full sequence of time steps\n",
    "\n",
    "# Data frame template defining the column names\n",
    "DataFrameTemplate = ('lon','lat','conv_rate','conv_obliq','migr_rate',\n",
    "                     'migr_obliq','arc_length','arc_azimuth',\n",
    "                     'subducting_plate','overriding_plate','time')\n",
    "\n",
    "# Create an empty dataframe to concatenate results to\n",
    "df_AllTimes = pd.DataFrame(columns=DataFrameTemplate)\n",
    "\n",
    "# list of reconstruction times\n",
    "reconstruction_times = np.arange(min_time,max_time+1,time_step)\n",
    "\n",
    "# Iterate over time steps\n",
    "for reconstruction_time in reconstruction_times:\n",
    "        \n",
    "    # Get the subduction kinematics database \n",
    "    subduction_data = []\n",
    "    output_data = sc.subduction_convergence(\n",
    "                rotation_filename,\n",
    "                input_topology_filename,\n",
    "                threshold_sampling_distance_radians,\n",
    "                reconstruction_time,\n",
    "                velocity_delta_time,\n",
    "                anchor_plate_id)\n",
    "\n",
    "    # Make a flat list of subduction stats to input into the proximity test\n",
    "    for data in output_data:\n",
    "        subduction_data.append(data+(reconstruction_time,))\n",
    "    \n",
    "    # convert list array to dataframe\n",
    "    df = pd.DataFrame(subduction_data, columns = DataFrameTemplate)\n",
    "    \n",
    "    ######################\n",
    "    # Age grid sampling\n",
    "    grdfile = '%s%0.0f%s' % (agegrid_filename[0],reconstruction_time,agegrid_filename[1])\n",
    "    print 'working on file %s' % grdfile\n",
    "    sample_result = sample_grid_using_scipy(df['lon'],df['lat'],grdfile)\n",
    "    \n",
    "    # Append interpolate seafloor ages to dataframe\n",
    "    df['SeafloorAge'] = pd.Series(sample_result)\n",
    "    \n",
    "    \n",
    "    ######################\n",
    "    # Find nearest continental polygon\n",
    "    reconstructed_continental_polygons = []\n",
    "    pygplates.reconstruct(continental_polygon_features,\n",
    "                          rotation_model,\n",
    "                          reconstructed_continental_polygons,\n",
    "                          reconstruction_time,\n",
    "                          anchor_plate_id=anchor_plate_id)\n",
    "    pid,dist = get_nearest_continental_polygon(df,reconstructed_continental_polygons)\n",
    "    df['NearestContinentPID'] = pd.Series(pid)\n",
    "    df['DistanceToContinent'] = pd.Series(dist)\n",
    "\n",
    "    # append dataframe \n",
    "    df_AllTimes = df_AllTimes.append(df)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Options for saving/loading results to/from file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save to csv file (prefered file format)\n",
    "df_AllTimes.to_csv('SubductionTable_clean_%0.0f_%0.0fMa.csv' % (min_time,max_time)) \n",
    "\n",
    "# Save to an hdf5 file\n",
    "#df_AllTimes.to_hdf('SubductionTable_%0.0f_%0.0fMa.h5' % (min_time,max_time),'SubductionTable') \n",
    "\n",
    "# Example of loading previously calculated results from a hdf5 file as created using the line above\n",
    "#df_AllTimes = pd.read_hdf('SubductionTable_0_230Ma.h5')\n",
    "\n",
    "#plt.scatter(df.lon,df.lat,c=df.DistanceToContinent,edgecolors='')\n",
    "#plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display a preveiw of the data we just created\n",
    "df_AllTimes"
   ]
  },
  {
   "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
}