hw_sumatraB_P_hf.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Applied Seismology, GEOS 626, University of Alaska Fairbanks\n",
    "\n",
    "- script hw_sumatraB_P_hf.ipynb\n",
    "- template script for analyzing the direct arrival waveforms from Sumatra:\n",
    "- channel BHZ, duration of 5 hours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "import warnings\n",
    "\n",
    "from lib_seis import smooth\n",
    "from lib_seis import get_dist_az\n",
    "from lib_seis import station_map_and_table\n",
    "from lib_seis import sumatra_event\n",
    "from obspy import read\n",
    "from obspy import Stream\n",
    "from pysep import Pysep\n",
    "from pysep.recsec import plotw_rs\n",
    "from scipy.signal import hilbert"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# script settings\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "plt.rcParams['figure.figsize'] = 10, 5\n",
    "plt.rcParams['lines.linewidth'] = 1\n",
    "\n",
    "class SEED_ID_Error(Exception):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Event details\n",
    "Mw 9.1 2004 Sumatra - Andaman Islands Earthquake \\\n",
    "https://earthquake.usgs.gov/earthquakes/eventpage/official20041226005853450_30"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# fetch and display event details\n",
    "\n",
    "event = sumatra_event()\n",
    "\n",
    "print()\n",
    "for i, item in enumerate(event):\n",
    "    print(f'{item} = {event[item]}')\n",
    "print()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# download BHZ data - response removed to displacement\n",
    "\n",
    "stations   = dict( networks                     = 'G,II,IU',\n",
    "                   stations                     = '*',\n",
    "                   locations                    = '*',\n",
    "                   channels                     = 'BHZ' )\n",
    "\n",
    "screening  = dict( mindistance_km               = 0  / 180 * np.pi * 6371,\n",
    "                   maxdistance_km               = 40 / 180 * np.pi * 6371,\n",
    "                   seconds_before_ref           = 1 * 60 * 60,\n",
    "                   seconds_after_ref            = 4 * 60 * 60,\n",
    "                   remove_clipped               = True,\n",
    "                   remove_insufficient_length   = True,\n",
    "                   fill_data_gaps               = 0,\n",
    "                   gap_fraction                 = 0.03,\n",
    "                   resample_freq                = 20     )\n",
    "\n",
    "processing = dict( remove_response              = True,\n",
    "                   output_unit                  = 'DISP',\n",
    "                   water_level                  = None,\n",
    "                   pre_filt                     = 'default',\n",
    "                   detrend                      = True,\n",
    "                   demean                       = True,\n",
    "                   taper_percentage             = 0.2 )\n",
    "\n",
    "output     = dict( taup_model                   = 'ak135',\n",
    "                   log_level                    = 'INFO',\n",
    "                   plot_files                   = 'map',\n",
    "                   output_dir                   = 'datawf',\n",
    "                   sac_subdir                   = '',\n",
    "                   overwrite_event_tag          = 'sumatra_BHZ_disp',\n",
    "                   overwrite                    = False )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# download data\n",
    "\n",
    "data_dir  = f'{output[\"output_dir\"]}/{output[\"overwrite_event_tag\"]}'\n",
    "overwrite = f'{output[\"overwrite\"]}'\n",
    "\n",
    "if (not os.path.isdir(data_dir)) or (overwrite == 'True'):\n",
    "    ps = Pysep(**event,**stations,**screening,**processing,**output)\n",
    "    ps.run()\n",
    "else:\n",
    "    print('\\nData directory already exists, no data will be downloaded\\n')"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# read downloaded raw BHZ data\n",
    "\n",
    "st = read(f'{data_dir}/*.sac', header=None)\n",
    "\n",
    "print('Details of downloaded waveforms - \\n')\n",
    "print(st.__str__(extended=True))\n",
    "print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event_path = f'{data_dir}/event.xml'\n",
    "station_map_and_table(st, event_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# OPTION A: PLOT ALL WAVEFORMS DOWNLOADED BY PYSEP\n",
    "# plot a record section using PySEP's record section plotting tool\n",
    "# change the following parameters as needed\n",
    "\n",
    "plot_info = dict( st                           = st,\n",
    "                  sort_by                      = 'distance',\n",
    "                  scale_by                     = 'normalize',\n",
    "                  time_shift_s                 = None,\n",
    "                  min_period_s                 = 1,\n",
    "                  max_period_s                 = 10,\n",
    "                  xlim_s                       = [-50, 7200],\n",
    "                  preprocess                   = 'st',\n",
    "                  max_traces_per_rs            = 40,\n",
    "                  distance_units               = 'deg',\n",
    "                  tmarks                       = [0],\n",
    "                  save                         = '',\n",
    "                  log_level                    = 'DEBUG')\n",
    "\n",
    "plotw_rs(**plot_info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# selecting waveforms\n",
    "\n",
    "# modify the following list of picked waveforms -\n",
    "picked_waveforms = ['IU.TATO.00.BHZ',\n",
    "                    'II.PALK.00.BHZ',\n",
    "                    'II.COCO.00.BHZ']\n",
    "\n",
    "st_select = Stream()\n",
    "\n",
    "for seed_id in picked_waveforms:\n",
    "    st_pick = st.select(id=seed_id)\n",
    "    if not st_pick:\n",
    "        raise SEED_ID_Error(f'{seed_id} not in inventory')\n",
    "    st_select += st_pick"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "event_path = f'{data_dir}/event.xml'\n",
    "station_map_and_table(st_select, event_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# OPTION B: PLOT SUBSET OF WAVEFORMS\n",
    "# plot a record section using PySEP's record section plotting tool\n",
    "# change the following parameters as needed\n",
    "\n",
    "plot_info = dict( st                           = st_select,\n",
    "                  sort_by                      = 'distance',\n",
    "                  scale_by                     = 'normalize',\n",
    "                  time_shift_s                 = None,\n",
    "                  min_period_s                 = 1,\n",
    "                  max_period_s                 = 10,\n",
    "                  xlim_s                       = [-50, 7200],\n",
    "                  preprocess                   = 'st',\n",
    "                  max_traces_per_rs            = 40,\n",
    "                  distance_units               = 'deg',\n",
    "                  tmarks                       = [0],\n",
    "                  save                         = '',\n",
    "                  log_level                    = 'DEBUG')\n",
    "\n",
    "plotw_rs(**plot_info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "language": "python",
   "display_name": "Python 3 (ipykernel)"
  },
  "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}