pds_travel_time_plot.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "###\n",
    "# Plot P and S wave epicentral distance vs arrival time for M5.0+ earthquakes from 1981-present\n",
    "# Data sources : s3://secedc-pds/earthquake_catalogs\n",
    "#                s3://scedc-pds/event_phases \n",
    "#\n",
    "# Author : Aparna Bhaskaran (aparnab at caltech.edu)\n",
    "###\n",
    "\n",
    "import pandas as pd\n",
    "import boto3\n",
    "import os\n",
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#S3 bucket, SCEDC public data set\n",
    "scedc_bucket = \"scedc-pds\" \n",
    "#prefix for index of the catalog in csv format\n",
    "catalog_prefix = \"earthquake_catalogs/index/csv/\" \n",
    "#prefix for phases\n",
    "phase_prefix = \"event_phases\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/aparnabhaskaran/anaconda3/lib/python3.6/site-packages/boto3/compat.py:88: PythonDeprecationWarning: Boto3 will no longer support Python 3.6 starting May 30, 2022. To continue receiving service updates, bug fixes, and security updates please upgrade to Python 3.7 or later. More information can be found here: https://aws.amazon.com/blogs/developer/python-support-policy-updates-for-aws-sdks-and-tools/\n",
      "  warnings.warn(warning, PythonDeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "s3 = boto3.resource('s3')\n",
    "bucket = s3.Bucket(scedc_bucket)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# retrieve the list of index files at s3://scedc-pds/earthquake_catalogs/index/csv/\n",
    "key = catalog_prefix \n",
    "index_list = list(bucket.objects.filter(Prefix=key))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# sample index csv \n",
    "#PREFIX,MS_FILENAME,PHASE_FILENAME,ORIGIN_DATETIME,ET,GT,MAG,M,LAT,LON,DEPTH,Q,EVID,NPH,NGRM\n",
    "#2022/2022_003/,39901471.ms,39901471.phase,2022-01-03 01:17:08.92,eq,l,0.80,l,32.864,-115.990,8.7,A,39901471,23,1321\n",
    "#2022/2022_003/,39901479.ms,39901479.phase,2022-01-03 01:33:27.31,eq,l,0.80,l,33.783,-117.844,2.9,A,39901479,18,1924\n",
    "#2022/2022_003/,39901487.ms,39901487.phase,2022-01-03 01:52:27.61,eq,l,0.84,l,36.151,-118.070,3.6,A,39901487,24,448\n",
    "#2022/2022_003/,39901495.ms,39901495.phase,2022-01-03 02:03:34.59,eq,l,0.45,l,33.482,-116.443,13.1,A,39901495,28,1544\n",
    "#2022/2022_003/,39901519.ms,39901519.phase,2022-01-03 02:34:48.20,eq,l,3.92,l,35.284,-119.391,17.7,A,39901519,109,3233       \n",
    "\n",
    "# from the list of indexes, only download the indexes for 1981-present\n",
    "for obj in index_list:\n",
    "    index_file = os.path.basename(obj.key)\n",
    "    if int(index_file.split('_')[0]) >= 1981:\n",
    "        #delete the csv if it exists\n",
    "        Path(index_file).unlink()\n",
    "        #download\n",
    "        bucket.download_file(obj.key, index_file)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#read the csv files and create a catalog\n",
    "catalog = pd.DataFrame()\n",
    "for csv in sorted(Path('.').glob(\"*.csv\")):\n",
    "    if catalog.empty:\n",
    "        catalog = pd.read_csv(csv)\n",
    "    else:\n",
    "        catalog = pd.concat([catalog, pd.read_csv(csv)])\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "57"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#query the catalog for local events (gtype = 'l') with magnitude >= 5.0 (MAG >= 5.0)\n",
    "temp = catalog.loc[(catalog['MAG'] >= 5.0) & (catalog['GT'] == 'l') , [\"PREFIX\",\"PHASE_FILENAME\",\"EVID\"]]\n",
    "\n",
    "# collect only those records that have a non-null PHASE_FILENAME entry in the catalog\n",
    "phases = temp[temp[\"PHASE_FILENAME\"].notna()]\n",
    "len(phases)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# download the phase files; s3://scedc-pds/event_phases/year=YYYY/YYYY_JJJ/evid.phase\n",
    "for i in phases.itertuples():\n",
    "    key = phase_prefix + \"/\" + i.PREFIX + i.PHASE_FILENAME\n",
    "    # delete the PHASE_FILENAME if it exists\n",
    "    Path(i.PHASE_FILENAME).unlink()\n",
    "    bucket.download_file(key, i.PHASE_FILENAME)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data structures to store tuples of (distance, arrival time) for P and S waves\n",
    "p = []\n",
    "s = []\n",
    "\n",
    "for phase_file in sorted(Path('.').glob(\"*.phase\")):\n",
    "    with open(phase_file,'r') as fp:\n",
    "        alllines = fp.readlines()\n",
    "    for line in alllines[1:]:\n",
    "        net, sta, seedchan, location, lat, lon, elev, phase, first_motion, sig_quality, pick_quality, distance, time = [item for item in line.strip().split(' ') if item != '']\n",
    "        if phase.lower() == 'p':\n",
    "            p.append((float(distance), float(time)))\n",
    "        if phase.lower() == 's':\n",
    "            s.append((float(distance), float(time)))\n",
    "            \n",
    "# sort in ascending order by distance            \n",
    "p.sort()\n",
    "s.sort()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import (MultipleLocator, AutoMinorLocator)\n",
    "import numpy as np\n",
    "#convert the tuples to numpy arrays for plotting\n",
    "# p-wave\n",
    "xp = np.asarray([item[0] for item in p])\n",
    "yp = np.asarray([item[1] for item in p])\n",
    "\n",
    "# s-wave\n",
    "xs = np.asarray([item[0] for item in s])\n",
    "ys = np.asarray([item[1] for item in s])\n",
    "\n",
    "# plot\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1,1,1)\n",
    "ax.plot(xp, yp, color='tab:blue', label=\"P Wave\")\n",
    "ax.plot(xs, ys, color='tab:orange', label=\"S Wave\")\n",
    "\n",
    "# set minor ticks\n",
    "ax.xaxis.set_minor_locator(MultipleLocator(20)) #km\n",
    "ax.yaxis.set_minor_locator(MultipleLocator(5))  #seconds\n",
    "ax.grid(b=True,which='major',color='black', linestyle='-')\n",
    "ax.set_xlabel(\"Distance from epicenter (km)\")\n",
    "ax.set_ylabel(\"Arrival time (seconds)\")\n",
    "ax.set_title(\"Travel time plot for M5.0+ earthquakes \\nin the SCEDC catalog (1981-present)\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "plt.savefig(\"SCEDC-PvsS-M5-1981-2022.jpg\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}