hw_sumatraA_CAN_response.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Applied Seismology, GEOS 626, University of Alaska Fairbanks\n",
    "\n",
    "- script hw_sumatraA_CAN_response.ipynb\n",
    "- template script for analyzing the amplitude spectra of the Sumatra earthquake\n",
    "- the example waveform is from CAN (Canberra, Australia) for the 2004 Mw 9.1 Sumatra-Andaman earthquake"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "\n",
    "import matplotlib.pylab as plt\n",
    "import numpy as np\n",
    "import warnings\n",
    "\n",
    "from lib_seis import markp_minutes\n",
    "from lib_seis import sumatra_event\n",
    "from lib_seis import wf_fft\n",
    "from matplotlib import dates\n",
    "from obspy.core import UTCDateTime\n",
    "from obspy.clients.fdsn import Client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# script settings\n",
    "\n",
    "warnings.filterwarnings('ignore')\n",
    "plt.style.use('ggplot')\n",
    "plt.rcParams['figure.figsize'] = 8, 6\n",
    "plt.rcParams['lines.linewidth'] = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Event details\n",
    "Mw 9.1 2004 Sumatra - Andaman Islands Earthquake \\\n",
    "https://earthquake.usgs.gov/earthquakes/eventpage/official20041226005853450_30"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# fetch and display event details\n",
    "\n",
    "event = sumatra_event()\n",
    "\n",
    "print()\n",
    "for item in event:\n",
    "    print(f'{item} = {event[item]}')\n",
    "print()    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# getting the waveform\n",
    "\n",
    "client = Client(\"IRIS\")\n",
    "\n",
    "starttime  = event['origin_time'] - ( 0.5 * 24 * 60 * 60 )           # half a day before origin time\n",
    "endtime    = event['origin_time'] + ( 9.5 * 24 * 60 * 60)            # 9 and a half days after origin time\n",
    "\n",
    "st = client.get_waveforms(\"G\", \"CAN\", \"*\", \"LHZ\", starttime, endtime, attach_response=True)\n",
    "print(st)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plotting the seismogram\n",
    "\n",
    "can = st.copy()                           \n",
    "\n",
    "t   = can[0].times()                                           # time axis for plotting in seconds\n",
    "tr  = can[0]\n",
    "\n",
    "title = f'{tr.id} − starting {tr.stats[\"starttime\"]}'\n",
    "\n",
    "# plot the seismogram with time in seconds\n",
    "plt.figure()\n",
    "plt.plot(t, tr.data, 'b')\n",
    "plt.title(title)\n",
    "#plt.ylim([-3E4, 3E4])\n",
    "plt.xlabel('time [s]')\n",
    "plt.ylabel('counts');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# label the x-axis using absolute time\n",
    "\n",
    "can_dates = st.copy().resample(0.1)\n",
    "\n",
    "date = can_dates[0].times(type='utcdatetime')                      # time axis for plotting absolute time\n",
    "date = dates.date2num(date)\n",
    "\n",
    "# plot the seismogram with respect to absolute time\n",
    "plt.figure()\n",
    "plt.plot_date(date, can_dates[0].data, 'b')\n",
    "plt.title(title)\n",
    "plt.xlabel('dates')\n",
    "plt.ylabel('counts')\n",
    "ax = plt.gca()\n",
    "ax.xaxis.set_tick_params(rotation=30, labelsize=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# spectral analysis\n",
    "\n",
    "can.detrend('constant')\n",
    "can.detrend('linear')\n",
    "can.taper(max_percentage=0.2, type=\"cosine\")\n",
    "\n",
    "tr = can[0]\n",
    "\n",
    "sampling_frequency = tr.stats.sampling_rate\n",
    "Nyquist_frequency  = sampling_frequency / 2\n",
    "\n",
    "fft_amp, fft_phase, fft_freq = wf_fft(tr.data, Nyquist_frequency)\n",
    "\n",
    "C = fft_amp * np.exp( 1j * fft_phase )                     # Fourier transform reconstructed from its amplitude and phase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# template code for using markp_minutes to identify the periods (in minutes) of spectral peaks\n",
    "\n",
    "#fig=plt.figure();\n",
    "#YOUR PLOT COMMAND HERE\n",
    "#cid = fig.canvas.mpl_connect('button_press_event', markp_minutes) #displays the value 1/x/60 for a selected point (x,y)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "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
}