{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "## Get some 40m data using NDS\n",
    "### do some analysis of the seismometer data\n",
    "#### import packages and set up the matplotlib environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false",
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Library Imports and Python parameter settings\n",
    "\n",
    "# make your jupyter lab awesome (lab >> notebook)\n",
    "# https://github.com/mauhai/awesome-jupyterlab\n",
    "\n",
    "# HowTo install nds2 client:\n",
    "# from within your Anaconda virtualenv:\n",
    "# conda install -c conda-forge python-nds2-client\n",
    "import nds2\n",
    "\n",
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.signal as sig\n",
    "import scipy.constants as const\n",
    "from astropy.time import Time\n",
    "\n",
    "from timeit import default_timer as timer\n",
    "\n",
    "# Update the matplotlib configuration parameters:\n",
    "plt.rcParams.update({'text.usetex': False,\n",
    "                     'lines.linewidth': 4,\n",
    "                     'font.family': 'serif',\n",
    "                     'font.serif': 'Georgia',\n",
    "                     'font.sans-serif': 'Helvetica',\n",
    "                     'font.size': 20,\n",
    "                     'xtick.direction': 'in',\n",
    "                     'ytick.direction': 'in',\n",
    "                     'xtick.labelsize': 'medium',\n",
    "                     'ytick.labelsize': 'medium',\n",
    "                     'axes.labelsize': 'medium',\n",
    "                     'axes.titlesize': 'medium',\n",
    "                     'axes.grid.axis': 'both',\n",
    "                     'axes.grid.which': 'both',\n",
    "                     'axes.grid': True,\n",
    "                     'grid.color': 'xkcd:grey',\n",
    "                     'grid.alpha': 0.53,\n",
    "                     'lines.markersize': 12,\n",
    "                     'legend.borderpad': 0.2,\n",
    "                     'legend.fancybox': True,\n",
    "                     'legend.fontsize': 'small',\n",
    "                     'legend.framealpha': 0.8,\n",
    "                     'legend.handletextpad': 0.5,\n",
    "                     'legend.labelspacing': 0.33,\n",
    "                     'legend.loc': 'best',\n",
    "                     'figure.figsize': ((12, 8)),\n",
    "                     'savefig.dpi': 140,\n",
    "                     'savefig.bbox': 'tight',\n",
    "                     'pdf.compression': 9})\n",
    "\n",
    "debug = True\n",
    "\n",
    "\n",
    "\n",
    "# this is a better performing downsampling function\n",
    "# cf. E Quintero's study: some repo\n",
    "def downsample(x, down_factor):\n",
    "    # Using fir_aa[1:-1] cuts off a leading and trailing zero\n",
    "    fir_aa = sig.firwin(20 * down_factor + 1, 0.8 / down_factor,\n",
    "                        window='blackmanharris')[1:-1]\n",
    "    out = sig.decimate(x, down_factor, ftype=sig.dlti(fir_aa, 1.0),\n",
    "                       zero_phase=True, axis=-1)\n",
    "\n",
    "    return out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false"
   },
   "source": [
    "### Get the data\n",
    "* https://www.lsc-group.phys.uwm.edu/daswg/projects/nds-client/doc/manual/ch04.html\n",
    "* http://wiki.scipy.org/NumPy_for_Matlab_Users"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false",
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "# Setup connection to the NDS (Network Data Server)\n",
    "conn = nds2.connection('nds40.ligo.caltech.edu', 31200)\n",
    "\n",
    "# Setup start and stop times\n",
    "times = '2020-04-22 07:03:00'   # UTC time of data start\n",
    "t = Time(times, format='iso', scale='utc')\n",
    "t_start = int(t.gps)\n",
    "print('Start time = ' + str(t_start))\n",
    "dur = 60 * 3  # data duration [s]\n",
    "\n",
    "# channel names\n",
    "channels = []\n",
    "chead = 'C1:PEM-SEIS_'\n",
    "\n",
    "seiss = ['BS', 'EX', 'EY']\n",
    "dofs  = ['X', 'Y', 'Z']\n",
    "ctail = '_OUT_DQ'\n",
    "for ii in seiss:\n",
    "    for jj in dofs:\n",
    "        channels.append(chead + ii + '_' + jj + ctail)\n",
    "\n",
    "#channels = ['C1:LSC-TRX_OUT_DQ', 'C1:ALS-BEATX_FINE_PHASE_OUT_DQ']\n",
    "channels = ['C1:IOO-MC_F_DQ']\n",
    "#channels = ['H1:OMC-DCPD_B_OUT_DQ']\n",
    "#channels = ['C1:SUS-BS_OL_SUM_OUT16.mean,m-trend']\n",
    "\n",
    "if debug:\n",
    "    print('Channel list:')\n",
    "    print(*channels,  sep='\\n')\n",
    "\n",
    "# time the data getting function\n",
    "t0 = timer()\n",
    "data = conn.fetch(t_start, t_start + dur, channels)\n",
    "t1 = timer()\n",
    "print('Data download time = ' + str(round(t1-t0,2)) + ' s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Downsample the Data\n",
    "The seismometer data is stored on disk at 256 Hz, but the Guralp and Trillium seismometers have no useful information above ~20 Hz.\n",
    "\n",
    "So we use a modified 'decimation' function to reduce the sampling rate:\n",
    "1. This makes the data more lightweight and speeds up all computations\n",
    "1. It may also be helpful to eliminate the high frequency data in the construction of the seismic 'Heat' map\n",
    "1. Decimation involves low-pass filtering and then downsampling\n",
    "1. we use a custom FIR window function to do the filtering because its better(?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false"
   },
   "outputs": [],
   "source": [
    "# downsample the data\n",
    "new_sample_rate = 32\n",
    "old_sample_rate = data[0].channel.sample_rate\n",
    "down_factor = int(old_sample_rate / new_sample_rate)\n",
    "\n",
    "vdata = np.empty([int(len(data[0].data)/down_factor), len(channels)])\n",
    "j = 0\n",
    "for x in data:\n",
    "    \n",
    "    vdata[:,j] = downsample(x.data, down_factor)\n",
    "    j +=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Make some plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "Collapsed": "false",
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "outputs": [],
   "source": [
    "nrows = 3\n",
    "ncols = 3\n",
    "fig1, ax = plt.subplots(nrows=nrows, ncols=ncols,sharex=True, sharey=True,\n",
    "                        gridspec_kw={'hspace': 0, 'wspace': 0}, figsize=(23,13))\n",
    "#plt.loglog(aligo[:,0], sqrt(aligo[:,1]), color='Indigo', ls='--', alpha=0.65, lw=4)\n",
    "fs = new_sample_rate\n",
    "print('Sample rate = ' + str(fs))\n",
    "#y  = data[0].data\n",
    "tt = np.arange(0, 1/fs*np.size(vdata[:,0]), 1/fs)\n",
    "\n",
    "# loop over the channels and plot each one\n",
    "for ii in range(nrows):\n",
    "    for jj in range(ncols):\n",
    "        n = jj + ii*nrows\n",
    "        ax[ii,jj].plot(tt, vdata[:,n], linestyle='-', lw=2, alpha=0.79,\n",
    "             label=channels[n][12:-7], rasterized=True, c = np.random.rand(3))\n",
    "        ax[ii,jj].legend()\n",
    "        \n",
    "\n",
    "#ax.set_xlabel(r'Time [s]', fontsize=18)\n",
    "#ax.set_ylabel(r'Velocity [um/s]', fontsize=18)\n",
    "\n",
    "#ax.legend()\n",
    "\n",
    "#plt.savefig(\"TRY.pdf\", bbox_inches='tight')\n",
    "for x in ax.flat:\n",
    "    x.label_outer()\n",
    "ax[2,1].set_xlabel('Time [s] after ' + times + 'Z')\n",
    "ax[1,0].set_ylabel('Velocity [$\\mu m/s$]')\n",
    "\n",
    "plt.savefig('April22-EQ.pdf')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "Collapsed": "false",
    "jupyter": {
     "outputs_hidden": false
    }
   },
   "source": [
    "## Make some spectrograms of the seismo signals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tfft = 8\n",
    "nfft = fs * tfft\n",
    "\n",
    "f, t, Sxx = sig.spectrogram(vdata[:,4], fs, window=sig.get_window(('tukey', 0.25), nfft))\n",
    "S = np.zeros((nrows, ncols, np.shape(Sxx)[0], np.shape(Sxx)[1]))\n",
    "\n",
    "for i in range(nrows):\n",
    "    for j in range(ncols):\n",
    "        n = j + i*nrows\n",
    "        f, t, Sxx = sig.spectrogram(vdata[:,n], fs, window=sig.get_window(('tukey', 0.25), nfft))\n",
    "        S[i,j,:,:] = Sxx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(nrows=nrows, ncols=ncols,sharex=True, sharey=True,\n",
    "                        gridspec_kw={'hspace': 0, 'wspace': 0}, figsize=(23,13))\n",
    "\n",
    "for i in range(nrows):\n",
    "    for j in range(ncols):\n",
    "        # remove the first index of f, since its the zero Hz bin\n",
    "        ax[i,j].pcolormesh(t, f[1:], np.log10(S[i,j,1:,:]), shading='gouraud', cmap=plt.get_cmap('inferno'))\n",
    "        ax[i,j].set_yscale('log')\n",
    "        \n",
    "for x in ax.flat:\n",
    "    x.label_outer()\n",
    "ax[2,1].set_xlabel('Time [s] after ' + times + 'Z')\n",
    "ax[1,0].set_ylabel('Frequency [Hz]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f[0]"
   ]
  },
  {
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}