lab_response.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Applied Seismology, GEOS 626, University of Alaska Fairbanks\n",
    "\n",
    "- script lab_response.ipynb"
   ]
  },
  {
   "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 warnings\n",
    "\n",
    "from obspy.clients.fdsn import Client\n",
    "from obspy.core import UTCDateTime"
   ]
  },
  {
   "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'] = 10, 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# getting the inventory\n",
    "\n",
    "client = Client(\"IRIS\")\n",
    "\n",
    "starttime = UTCDateTime(2004, 12, 26, 0, 58, 50)\n",
    "endtime   = starttime + ( 5 * 24 * 60 * 60 )        # 5 days after the start time\n",
    "\n",
    "# check the page - https://docs.obspy.org/packages/obspy.core.inventory.html, \n",
    "# for detailed information on the inventory object\n",
    "\n",
    "waveforms = ['G.CAN.'    ,\n",
    "             'II.ALE.00' ,\n",
    "             'IU.ANMO.00']\n",
    "\n",
    "# choose a waveform based on the network whose response you want to look at\n",
    "picked_waveform = waveforms[0]\n",
    "\n",
    "net, sta, loc = picked_waveform.split('.')\n",
    "\n",
    "inv_LHZ = client.get_stations(network=net, station=sta, location=loc, channel=\"LHZ\", starttime=starttime,\n",
    "                              endtime=endtime, level=\"response\")\n",
    "inv_BHZ = client.get_stations(network=net, station=sta, location=loc, channel=\"BHZ\", starttime=starttime,\n",
    "                              endtime=endtime, level=\"response\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# see details about the instrument response\n",
    "\n",
    "# set inventory for the channel of interest\n",
    "#inv = inv_LHZ.copy()\n",
    "inv = inv_BHZ.copy()\n",
    "\n",
    "print(inv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(inv[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(inv[0][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(inv[0][0][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# collect key values from the inventory\n",
    "\n",
    "samples_per_second = inv[0][0][0].sample_rate\n",
    "dt = 1/samples_per_second\n",
    "Nyquist_frequency  = 1 / (2*dt)\n",
    "print(f'sample rate = {samples_per_second} 1/s, dt = {dt} s, Nyquist frequency = {Nyquist_frequency} Hz \\n')\n",
    "\n",
    "print(inv[0][0][0].response)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the page - https://docs.obspy.org/packages/autogen/obspy.core.inventory.inventory.Inventory.plot_response.html,\n",
    "# for detailed information on plotting response using this method\n",
    "\n",
    "#inv.plot_response?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#inv = inv_LHZ.copy()\n",
    "\n",
    "# plot the instrument response (note: semi-colon is needed to block a repeated plot)\n",
    "inv.plot_response(min_freq=0.0001);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# optional: plot the stages separately\n",
    "#inv[0][0][0].plot(min_freq=0.0001, start_stage=1, end_stage=1, output=\"VEL\");\n",
    "#inv[0][0][0].plot(min_freq=0.0001, start_stage=2, end_stage=2, output=\"VEL\");\n",
    "#inv[0][0][0].plot(min_freq=0.0001, start_stage=3, end_stage=3, output=\"VEL\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from here on out, we will specify the specific frequencies for which we will evaluate the instrument response\n",
    "f = np.logspace(-4.0, 2.0, num=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# retrieving the combined instrument response for all stages\n",
    "\n",
    "start_stage = None\n",
    "end_stage = None\n",
    "\n",
    "Id = inv[0][0][0].response.get_evalresp_response_for_frequencies(f, output='DISP', start_stage= start_stage,\n",
    "                                                                 end_stage=end_stage)\n",
    "Iv = inv[0][0][0].response.get_evalresp_response_for_frequencies(f, output='VEL', start_stage= start_stage,\n",
    "                                                                 end_stage=end_stage)\n",
    "Ia = inv[0][0][0].response.get_evalresp_response_for_frequencies(f, output='ACC', start_stage= start_stage,\n",
    "                                                                 end_stage=end_stage)\n",
    "\n",
    "max_Id   = max(Id)\n",
    "f_Id_max = f[np.nonzero(Id == max_Id)]\n",
    "max_Iv   = max(Iv)\n",
    "f_Iv_max = f[np.nonzero(Iv == max_Iv)]\n",
    "max_Ia   = max(Ia)\n",
    "f_Ia_max = f[np.nonzero(Ia == max_Ia)]\n",
    "\n",
    "fig, ax = plt.subplots(3, 2, figsize=(14,16))\n",
    "\n",
    "ax[0,0].semilogx(f, np.angle(Id)*180/np.pi, color='b')\n",
    "ax[0,0].set(xlabel='Frequency (Hz)', ylabel='Phase (deg)')\n",
    "ax[0,0].axvline(Nyquist_frequency, c='r', lw=0.5, ls='--', dashes=[12,12])\n",
    "\n",
    "ax[0,1].loglog(f, abs(Id), color='b')\n",
    "ax[0,1].set(xlabel='Frequency (Hz)', ylabel='Amplitude (m Hz$^-$$^1$)')\n",
    "ax[0,1].axvline(Nyquist_frequency, c='r', lw=0.5, ls='--', dashes=[12,12])\n",
    "ax[0,1].set_title('(m to counts) max = %.2e at %.2e Hz'% (np.real(max_Id), f_Id_max), fontsize=11)\n",
    "\n",
    "ax[1,0].semilogx(f, np.angle(Iv)*180/np.pi, color='b')\n",
    "ax[1,0].set(xlabel='Frequency (Hz)', ylabel='Phase (deg)')\n",
    "ax[1,0].axvline(Nyquist_frequency, c='r', lw=0.5, ls='--', dashes=[12,12])\n",
    "\n",
    "ax[1,1].loglog(f, abs(Iv), color='b')\n",
    "ax[1,1].set(xlabel='Frequency (Hz)', ylabel='Amplitude (ms$^-$$^1$Hz$^-$$^1$)')\n",
    "ax[1,1].axvline(Nyquist_frequency, c='r', lw=0.5, ls='--', dashes=[12,12])\n",
    "ax[1,1].set_title('(m/s to counts) max = %.2e at %.2e Hz'% (np.real(max_Iv), f_Iv_max), fontsize=11)\n",
    "\n",
    "ax[2,0].semilogx(f, np.angle(Ia)*180/np.pi, color='b')\n",
    "ax[2,0].set(xlabel='Frequency (Hz)', ylabel='Phase (deg)')\n",
    "ax[2,0].axvline(Nyquist_frequency, c='r', lw=0.5, ls='--', dashes=[12,12])\n",
    "\n",
    "ax[2,1].loglog(f, abs(Ia), color='b')\n",
    "ax[2,1].set(xlabel='Frequency (Hz)', ylabel='Amplitude (ms$^-$$^2$Hz$^-$$^1$)')\n",
    "ax[2,1].set_title('(m/s$^2$ to counts) max = %.2e at %.2e Hz'% (np.real(max_Ia), f_Ia_max), fontsize=11)\n",
    "ax[2,1].axvline(Nyquist_frequency, c='r', lw=0.5, ls='--', dashes=[12,12])\n",
    "\n",
    "fig.suptitle(f'Network - {inv[0].code}, Station - {inv[0][0].code}, Location - {inv[0][0][0].location_code}, Channel - {inv[0][0][0].code}', fontsize=12)\n",
    "\n",
    "plt.subplots_adjust(hspace=0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# superimposing BHZ and LHZ instrument responses for velocity\n",
    "# NOTE - the plots below represent the combined instrument response for all stages\n",
    "\n",
    "Iv_BHZ = inv_BHZ[0][0][0].response.get_evalresp_response_for_frequencies(f, output='VEL')\n",
    "Iv_LHZ = inv_LHZ[0][0][0].response.get_evalresp_response_for_frequencies(f, output='VEL')\n",
    "\n",
    "fNy_BHZ = (inv_BHZ[0][0][0].sample_rate) / 2\n",
    "fNy_LHZ = (inv_LHZ[0][0][0].sample_rate) / 2\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 8))\n",
    "\n",
    "ax1.semilogx(f, np.angle(Iv_BHZ)*180/np.pi, color='b', label='BHZ')\n",
    "ax1.semilogx(f, np.angle(Iv_LHZ)*180/np.pi, color='r', label='LHZ')\n",
    "ax1.axvline(fNy_BHZ, c='b', lw=0.5, ls='--', dashes=[12,12])\n",
    "ax1.axvline(fNy_LHZ, c='r', lw=0.5, ls='--', dashes=[12,12])\n",
    "ax1.set(xlabel='Frequency (Hz)', ylabel='Phase (degrees)')\n",
    "ax1.legend(loc='lower left', shadow=True)\n",
    "\n",
    "ax2.loglog(f, abs(Iv_BHZ), color='b', label='BHZ')\n",
    "ax2.loglog(f, abs(Iv_LHZ), color='r', label='LHZ')\n",
    "ax2.axvline(fNy_BHZ, c='b', lw=0.5, ls='--', dashes=[12,12])\n",
    "ax2.axvline(fNy_LHZ, c='r', lw=0.5, ls='--', dashes=[12,12])\n",
    "ax2.set(xlabel='Frequency (Hz)', ylabel='Amplitude (ms$^-$$^1$Hz$^-$$^1$)')\n",
    "ax2.legend(loc='lower left', shadow=True)\n",
    "\n",
    "fig.suptitle(f'Network - {inv[0].code}, Station - {inv[0][0].code}, Location - {inv[0][0][0].location_code}', fontsize=12)"
   ]
  }
 ],
 "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
}