2020_04_29_practical.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "2020-04-29_practical.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyPVCV+nBkQe8M58IM+G2+NK",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/slapin81/PhaseNet_practical/blob/master/2020_04_29_practical.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wqC4tjNbp4Ta",
        "colab_type": "text"
      },
      "source": [
        "# **Using PhaseNet**\n",
        "## 29/04/2020 ML reading group practical\n",
        "\n",
        "Original PhaseNet paper: https://doi.org/10.1093/gji/ggy423\n",
        "\n",
        "Google Colab is a great resource for small deep learning projects or creating sharable notebooks. It comes with many pre-installed Python modules (pandas, numpy, matplotlib, tensorflow, keras, ...) but you can also install others (e.g., we install obspy below). This makes it a really quick environment to just dive into some deep learning, or bash together any kind of Python work/tutorial really, and will ensure all users are working with the same setup and that everything should work for everyone.\n",
        "\n",
        "A couple useful bits to know about Google Colab notebooks:\n",
        "- In the top right corner, you can see RAM and Disk usage (if you hover the mouse over it, you will see exactly how much of each you have). At the time of writing this, the default was to give users 12 GB RAM (this should automatically go up to 25 GB if you exceed it, although you will have to restart the notebook) and differing amounts of disk space depending on whether you are using a GPU, TPU or nothing (see next point)\n",
        "- If you look at **Runtime** -> **Change runtime** in taskbar menu at the top, you will see available hardware acceleration options. This allows you to accelearate your deep learning code with a **GPU** or **TPU**. If you don't actually need either of these, you should set the hardware acclearation option to **None**. There are only a set no. of GPU's and TPU's for all users in the world so sometimes you might not see GPU or TPU in this menu, hence why it's important not to select one if you don't need it. At the time of writing, selecting **None** gives you ~ 108 GB disk space, selecting **GPU** gives you ~ 68 GB of disk space, and selecting **TPU** also gives you ~ 108 GB. Around 30 GB appears to be used by the notebook's backend. \n",
        "\n",
        "In this tutorial, we will need to use a GPU, but we will get to that in a minute.\n",
        "\n",
        "Anyway, let's get on with it..."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G-6nuiiWdtD6",
        "colab_type": "text"
      },
      "source": [
        "# 1. Install obspy using pip\n",
        "\n",
        "Google Colab does not come with obspy so we need to install it. This will allow us to retrieve waveforms, read/write mseed files, etc.\n",
        "\n",
        "You can install modules that aren't \"shipped\" with Google Colab by using `pip` in the shell (i.e., outside of Python). You can pass commands to the shell by putting a `!` at the start of any command."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mfl9xmIGKB3n",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!pip install obspy"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5xHgcLlOXnKj",
        "colab_type": "text"
      },
      "source": [
        "# 2. Restart runtime\n",
        "For some reason, you need to restart runtime after installing obspy or it doesn't work. So, once obspy has finished installing, click **Runtime** on the taskbar at the top, then **Restart runtime...** and then run the first cell of code above again before continuing.\n",
        "\n",
        "If you have done it correctly, you should see something like \n",
        "\n",
        "```\n",
        "Requirement already satisfied: obspy in /usr/local/lib/python3.6/dist-packages (1.2.1)\n",
        "Requirement already satisfied: requests in /usr/local/lib/python3.6/dist-packages (from obspy) (2.21.0)\n",
        "Requirement already satisfied: setuptools in /usr/local/lib/python3.6/dist-packages (from obspy) (46.1.3)\n",
        "Requirement already satisfied: lxml in /usr/local/lib/python3.6/dist-packages (from obspy) (4.2.6)\n",
        "Requirement already satisfied: future>=0.12.4 in /usr/local/lib/python3.6/dist-packages (from obspy) (0.16.0)\n",
        "Requirement already satisfied: numpy>=1.6.1 in /usr/local/lib/python3.6/dist-packages (from obspy) (1.18.3)\n",
        "Requirement already satisfied: scipy>=0.9.0 in /usr/local/lib/python3.6/dist-packages (from obspy) (1.4.1)\n",
        "Requirement already satisfied: decorator in /usr/local/lib/python3.6/dist-packages (from obspy) (4.4.2)\n",
        "Requirement already satisfied: matplotlib>=1.1.0 in /usr/local/lib/python3.6/dist-packages (from obspy) (3.2.1)\n",
        "Requirement already satisfied: sqlalchemy in /usr/local/lib/python3.6/dist-packages (from obspy) (1.3.16)\n",
        "Requirement already satisfied: urllib3<1.25,>=1.21.1 in /usr/local/lib/python3.6/dist-packages (from requests->obspy) (1.24.3)\n",
        "Requirement already satisfied: idna<2.9,>=2.5 in /usr/local/lib/python3.6/dist-packages (from requests->obspy) (2.8)\n",
        "Requirement already satisfied: chardet<3.1.0,>=3.0.2 in /usr/local/lib/python3.6/dist-packages (from requests->obspy) (3.0.4)\n",
        "Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.6/dist-packages (from requests->obspy) (2020.4.5.1)\n",
        "Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib>=1.1.0->obspy) (2.4.7)\n",
        "Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib>=1.1.0->obspy) (1.2.0)\n",
        "Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.6/dist-packages (from matplotlib>=1.1.0->obspy) (0.10.0)\n",
        "Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib>=1.1.0->obspy) (2.8.1)\n",
        "Requirement already satisfied: six in /usr/local/lib/python3.6/dist-packages (from cycler>=0.10->matplotlib>=1.1.0->obspy) (1.12.0)\n",
        "``` \n",
        "\n",
        "as output from that cell when you run it a second time."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Kc3G0ma4Ky9Q",
        "colab_type": "text"
      },
      "source": [
        "*******************************************"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-Zhvn4XWZiU4",
        "colab_type": "text"
      },
      "source": [
        "# 3. Select GPU\n",
        "\n",
        "For this practical we need to use a GPU. Using a GPU will speed up model training and prediction but, more importantly, the PhaseNet model only seems to work when GPU is selected.\n",
        "\n",
        "Go to **Runtime** then **Change runtime type** and select **GPU**.\n",
        "\n",
        "# 4. Set version of Tensorflow to 1.x\n",
        "\n",
        "Google Colab uses Tensorflow **v2.x** as default. However, PhaseNet was produced using Tensorflow **v1.10.1** and doesn't work with version 2.x. Run the following code cell to tell Google Colab to use Tensorflow **v1.x**."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "413E0cVKW0Xa",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# PhaseNet was produced using TensorFlow (TF) v1.10.1, so we need to tell Google Colab to use v1:\n",
        "%tensorflow_version 1.x"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IjOrwPEAiZE_",
        "colab_type": "text"
      },
      "source": [
        "# 5. Clone PhaseNet Github repo"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QT5x2HJXib6W",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!git clone https://github.com/wayneweiqiang/PhaseNet.git"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SWjwIGxiwF1L",
        "colab_type": "text"
      },
      "source": [
        "The code line above retrieves the Github repo for PhaseNet so we can use their scripts. The code is stored in a folder called PhaseNet. You can see the contents of this folder using the `ls` command in the shell:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7URl6JZPwade",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!ls ./PhaseNet"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X8Cc03E0dm3X",
        "colab_type": "text"
      },
      "source": [
        "Alternatively, you can look through your Google Colab file structure by clicking the folder tab ![](https://drive.google.com/uc?id=12m584puiNNQvXqbAgwZTWl0FrvtJ6ZNh)on the left hand panel.\n",
        "\n",
        "\n",
        "# 6. Get PhaseNet demo data using obspy\n",
        "\n",
        "In the PhaseNet directory structure above, there is a directory called `demo`. In that directory, there is a notebook that retrieves data from SCEDC (Southern California Earthquake Data Center) for running PhaseNet on. For ease, I have copied out the steps from that demo notebook below, changing them slightly to work better here (e.g., we request data from NCEDC, the Northern California Earthquake Data Center, as its servers were more reliable when I wrote this notebook):"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sPJ8sO9CXFVz",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Import required modules for getting data and saving to mseed files\n",
        "import os\n",
        "import obspy\n",
        "from obspy import UTCDateTime\n",
        "from obspy.clients.fdsn import Client"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "d7urrGT4dh9x",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Start client for requesting data from NCEDC\n",
        "client = Client(\"NCEDC\")\n",
        "\n",
        "# Create directory called mseed to store data\n",
        "data_dir = \"mseed\"\n",
        "if not os.path.exists(data_dir):\n",
        "  os.makedirs(data_dir)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6NmwcVfPdpKt",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Give start and end time to extract data\n",
        "starttime = UTCDateTime(\"2019-07-04T17:00:00\")\n",
        "endtime = UTCDateTime(\"2019-07-05T00:00:00\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "WsNCeAGOuZAv",
        "colab_type": "text"
      },
      "source": [
        "The start and end times given above are from the PhaseNet demo notebook and encompass the Ridgecrest 4th July foreshock (17:33:49 UTC)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "r6ZsngtXeFXq",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Get waveforms for stations CCC and GSC from CI network between start and end times above\n",
        "CCC = client.get_waveforms(\"CI\", \"CCC\", \"*\", \"HHE,HHN,HHZ\", starttime, endtime)\n",
        "GSC = client.get_waveforms(\"CI\", \"GSC\", \"*\", \"HHE,HHN,HHZ\", starttime, endtime)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-rWiru-3E2Fz",
        "colab_type": "text"
      },
      "source": [
        "Lets check the data was requested successfully by plotting it:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sy7wZChBE_PJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Plot data from station CCC\n",
        "CCC.plot()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bIf5_c23FCOk",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Plot data from station GSC\n",
        "GSC.plot()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "CAShrN9yFNoG",
        "colab_type": "text"
      },
      "source": [
        "Hopefully you should see some 3-component waveforms for each station. When I run it, the notebook seems to plot each set of 3-component data twice - don't worry if that happens, it's still fine."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NIsH-ywZzZAq",
        "colab_type": "text"
      },
      "source": [
        "# 7. Create input files for PhaseNet (mseed and filelist)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "n4c8-iwY1ekK",
        "colab_type": "text"
      },
      "source": [
        "The PhaseNet Github repo (which we cloned in step 5) has a script called `run.py`. This runs the PhaseNet model on either `.mseed` files or NumPy arrays (saved as `.npz` files). We'll start by saving the NCEDC data as `mseed` files and running the model on those."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9uBt_4mIeKOc",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "CCC.write(os.path.join(data_dir, \"CCC.mseed\")) # Write CCC station data to CCC.mseed\n",
        "GSC.write(os.path.join(data_dir, \"GSC.mseed\")) # Write GSC station data to GSC.mseed"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "IE_7QioV19cY",
        "colab_type": "text"
      },
      "source": [
        "Their `run.py` script also requires a file list for the model, in `csv` format, with four columns:\n",
        "- filename (mseed or npz file)\n",
        "- East component name\n",
        "- North component name\n",
        "- Vertical component name\n",
        "\n",
        "The following code cell makes that file list for you."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2shae_gDgOK2",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "with open(\"fname.csv\", 'w') as fp:\n",
        "  fp.write(\"fname,E,N,Z\\n\")\n",
        "  fp.write(\"CCC.mseed,HHE,HHN,HHZ\\n\")\n",
        "  fp.write(\"GSC.mseed,HHE,HHN,HHZ\\n\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "c1PnR42lGBB8",
        "colab_type": "text"
      },
      "source": [
        "We can see what's in that `csv` file by running `cat` in the shell:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ypu_y-z1jGtQ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Print what's in fname.csv:\n",
        "!cat fname.csv"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Hjd_1SgTGJP_",
        "colab_type": "text"
      },
      "source": [
        "We can see above that the file list has a header line with comma-separated values (`fname`, `E`, `N` and `Z`) followed by the filenames and component names for our mseed files on subsequent lines."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gsNvOVpujCJe",
        "colab_type": "text"
      },
      "source": [
        "# 8. Running PhaseNet (the easy way)\n",
        "\n",
        "PhaseNet was produced using TensorFlow (`tf`) version 1.x. This has a rather clunky process of having to initialise a 'session', `tf.Session()`, and creating lots of placeholders and variables to accomodate your batches, weights, biases, etc before you can create or run a model. Here is an example of code to train a simple neural network with **no** hidden layers in `tf` version 1.x:\n",
        "\n",
        "![alt text](https://drive.google.com/uc?id=1DkKr5lTLx56sayBXcOMlZPNjsKw80qop)\n",
        "\n",
        "\n",
        "Keras greatly streamlines this process (placeholders, matrix multiplication, reshaping, converting to float, telling optimizer to minimise loss function, etc are all done automatically), making it far more user friendly. The Keras equivalent of the TF code above would be something like:\n",
        "\n",
        "```\n",
        "import tensorflow as tf\n",
        "import keras\n",
        "from keras.models import *\n",
        "from keras.layers import *\n",
        "from keras.optimizers import *\n",
        "from keras.losses import *\n",
        "\n",
        "input_size = (28, 28, 1) # Input B/W image size (28 x 28 pixels x 1 channel)\n",
        "inputs = Input(input_size) # Input layer\n",
        "outputs = Dense(10, activation=\"softmax\")(inputs) # Output layer (10 nodes, 0 to 9 for handwritten digits)\n",
        "\n",
        "# Create model\n",
        "model = Model(inputs = inputs, outputs = outputs) \n",
        "\n",
        "# Compile model with optimizer and loss function\n",
        "model.compile(optimizer = SGD(0.003), loss = categorical_crossentropy, metrics = ['accuracy'])\n",
        "\n",
        "# Train model\n",
        "model.fit(x=mnist.train.images, y=mnist.train.labels, validation_data=(mnist.test.images, mnist.test.labels), batch_size=100, epochs=10000)\n",
        "\n",
        "```\n",
        "\n",
        "This is basically 6 lines of easy to follow code in Keras vs. ~ 20 in TF.\n",
        "\n",
        "Frustratingly, PhaseNet was not produced using Keras and the authors' Github code is rather convoluted. They have embedded all the extra code you need in Tensorflow v1.x within their own custom functions, accepting user inputted arguments to create all the required variables etc. \n",
        "\n",
        "Depending on how you look at it, this makes it either easier or harder to use: easier if you just use their scripts/functions, harder if you want to read the model in and play around with it (e.g., for fine-tuning / transfer learning etc).\n",
        "\n",
        "For now, we'll just use their pre-made Python script for running PhaseNet (`run.py`) as that will be far simpler. We could look at the Tensorflow code in more detail after, or try and recreate the model in Keras and do our own training for some practice at that.\n",
        "\n",
        "In the code cell below, we run their Python script (located at `PhaseNet/run.py`) with a number of arguments:\n",
        "- We choose `--mode=pred`, indicating that the model should be run in prediction mode, not training mode\n",
        "- `--model_dir` gives the location of the PhaseNet model within the directory structure from the authors' Github repo\n",
        "- `--data_dir` is where our data from stations CCC and GSC were saved as mseed files in step 7\n",
        "- `--data_list` is the file list we created in step 7\n",
        "- `--output_dir` is where the script should save the PhaseNet output (P and S pick times and probabilities)\n",
        "- `--batch_size` is how many waveforms at a time the model should be run on - this is not so important in prediction mode, but is a rather important hyperparameter during model training.\n",
        "- `--input_mseed` is to indicate we are inputting mseed files, rather than numpy arrays."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wGY52NqmJAOD",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!python PhaseNet/run.py --mode=pred --model_dir=PhaseNet/model/190703-214543 --data_dir=mseed --data_list=fname.csv --output_dir=output --batch_size=20 --input_mseed"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Rv9vWlRhbRGl",
        "colab_type": "text"
      },
      "source": [
        "The above step should have run the model on our two mseed files from stations CCC and GSC. The output of that step should be a `picks.csv` file in the directory called `output` containing pick times. Here are some notes from the PhaseNet repo describing what the above step does in more detail:\n",
        "\n",
        "- The detected P&S phases are stored to file `picks.csv` inside `--output_dir`. The `picks.csv` has three columns: file name with the beginning sample index (`fname_index`), P-phase index (`itp`), P-phase probability (`tp_prob`), S-phase index (`its`), S-phase probability (`ts_prob`). The absolute phase index = `fname_index` + `itp` or `its`.\n",
        "- The activation thresholds for P&S phases are set to 0.3 as default. Specify `--tp_prob` and `--ts_prob` to change the two thresholds.\n",
        "- On default, the mseed file is processed twice with 50% overlap to avoid phases being cut in the middle. The second pass are appended to the end of first pass. For example, if the index of the input data is from 0-60000, the index of second pass is from 60000-120000. If the processing window is 3000, the fist 1500 samples of the second pass are the padded zeros.\n",
        "\n",
        "So, as we set `--output_dir=output`, our picks should be stored in `output/picks.csv`, let's have a look:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "wekZvIfGS3Rb",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# First import modules for examining, manipulating and plotting picks\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "import math\n",
        "import matplotlib.pyplot as plt\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "z1slsFopbiox",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Read in output/picks.csv and print to screen\n",
        "picks_raw = pd.read_csv(\"output/picks.csv\")\n",
        "picks_raw"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "arsKUXWBcFZB",
        "colab_type": "text"
      },
      "source": [
        "Wow, this output format is really quite bad and involves a fair bit of wrangling to provide anything useful in terms of creating a catalogue of picks.\n",
        "\n",
        "The files have been broken into 30 second chunks and PhaseNet has been run on each one separately.\n",
        "\n",
        "The absolute index values with regards to our original mseed input files are given under `fname` (e.g., the first two files start at index 0, then index 3000, or 30 secs in, then index 6000, or 60 secs in, and so on...).\n",
        "\n",
        "The pick times are given as index values relative to the start of each 30 second chunk under headers `itp` or `its`, depending on whether its a P-phase pick or an S-phase pick. These are saved as strings surrounded by square brackets `[]`, with a space in between them if there are multiple picks. Pick probabilities between 0 and 1 are also stored as strings inside square brackets, but under headers `tp_prob` and `ts_prob` for P- and S-phase respectively. These correspond with the pick times given in `itp` and `its` if there are multiple picks.\n",
        "\n",
        "Just to make it a little more difficult to work with, the `run.py` script runs the model twice with 50% overlap but doesn't put these overlapping windows in order. Rather, they run the model over consecutive non-overlapping 30-second chunks, pad the original mseed file with 15 secs of zeros, then run it again with index values in the `fname` column following on from the index value at the end of the first pass. \n",
        "\n",
        "This means that the absolute index values under `fname` will run to approximately twice the number of actual index values (e.g., in this case, our input files only have 2,520,000 samples but the index values under `fname` run to 5,037,000). This also means that the same point in our input data is represented with two different index values (one < 2,520,000 from the first pass, one > 2,520,000 from the second pass). This creates a slight issue when the model has picked a single phase twice (once on each pass), which we will see below.\n",
        "\n",
        "# 9. Wrangling PhaseNet output\n",
        "\n",
        "Let's make a new dataframe with station name, UTC phase time, phase type and phase probability.\n",
        "\n",
        "To do this, we'll need the sample interval(s) (`delta`), the number of points (`npts`) and the UTC start time for our input files:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QMMHGNNoG2re",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# deltas, npts and start times for each station (in case different)\n",
        "ccc_delta = CCC.traces[0].stats['delta']\n",
        "ccc_npts = CCC.traces[0].stats['npts']\n",
        "ccc_start = CCC.traces[0].stats['starttime']\n",
        "\n",
        "gsc_delta = GSC.traces[0].stats['delta']\n",
        "gsc_npts = GSC.traces[0].stats['npts']\n",
        "gsc_start = GSC.traces[0].stats['starttime']\n",
        "\n",
        "# Check these values by printing to screen\n",
        "(ccc_delta, ccc_npts, ccc_start, gsc_delta, gsc_npts, gsc_start)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XiTrYnhoXrpC",
        "colab_type": "text"
      },
      "source": [
        "There's probably a more Pythonic way to extract the necessary info from the raw pick output and convert to UTC pick times, but here we'll just loop through each row of `picks_raw`, creating a new row in another dataframe, called `picks`, when there are P- or S-phase picks:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YO1t_PLldUfL",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Empty pick dataframe\n",
        "picks = pd.DataFrame(columns=['sta', 'time', 'pha', 'prob'])\n",
        "\n",
        "for i in range(picks_raw.shape[0]):\n",
        "  \n",
        "  # Get absolute start index value\n",
        "  abs_index = int(picks_raw.iloc[i]['fname'].split(\"_\")[-1])\n",
        "\n",
        "  # Indicate 0 if first pass (no overlap), 1 if second pass (50% overlap)\n",
        "  which_pass = 1 if abs_index >= ccc_npts else 0 # npts the same for both stations\n",
        "  \n",
        "  # If itp is not empty, i.e., if there is a P-phase pick...\n",
        "  if picks_raw.iloc[i]['itp'] != '[]':\n",
        "    \n",
        "    new_start_time = (abs_index % ccc_npts) * ccc_delta # Start time of 30-sec window in seconds since UTC starttime, delta is the same for both stations\n",
        "    # The % modulo operator above is because it's run over the data twice\n",
        "    \n",
        "    new_pick_times = [new_start_time + (t * ccc_delta) - (which_pass * 15) for t in list(map(float, picks_raw.iloc[i]['itp'][1:-1].split()))]\n",
        "    # The line above:\n",
        "    # i) removes the square brackets from itp by taking all but the first and last elements of string\n",
        "    # ii) splits the string by blank space, giving us indices for each pick in a row of picks_raw but still as strings\n",
        "    # iii) maps (converts) the strings to floats (creating a map object)\n",
        "    # iv) converts map object to list\n",
        "    # v) multiplies pick index by delta to get pick time in seconds, relative to start of 30-sec window\n",
        "    # vi) adds aboslute start time of 30 sec window so that you have a time for pick in seconds since UTC starttime\n",
        "    \n",
        "    new_pick_probs = list(map(float, picks_raw.iloc[i]['tp_prob'][1:-1].split())) # Do the same as above but converting phase probability strings to list of floats\n",
        "\n",
        "    # Create list where first element is station name, second element is UTC pick time, third element is phase type (p or s) and fourth element is phase probability\n",
        "    data = [[picks_raw.iloc[i]['fname'].split(\".\")[0], starttime + pick_time, \"p\", pick_prob] for pick_time, pick_prob in zip(new_pick_times, new_pick_probs)] \n",
        "\n",
        "    # Append any picks to picks dataframe\n",
        "    picks = picks.append(pd.DataFrame(data, columns=['sta', 'time', 'pha', 'prob']), ignore_index=True)\n",
        "  \n",
        "  # If its is not empty, i.e., if there is an S-phase pick:\n",
        "  if picks_raw.iloc[i]['its'] != '[]':\n",
        "\n",
        "    # Repeat above using its and ts_prob instead:\n",
        "\n",
        "    new_start_time = (abs_index % ccc_npts) * ccc_delta # Start time of 30-sec window in seconds since UTC starttime\n",
        "    new_pick_times = [new_start_time + (t * ccc_delta) - (which_pass * 15) for t in list(map(float, picks_raw.iloc[i]['its'][1:-1].split()))] # Convert pick index values to list of pick times in seconds since UTC starttime\n",
        "    new_pick_probs = list(map(float, picks_raw.iloc[i]['ts_prob'][1:-1].split())) # Convert phase probability strings to list of floats\n",
        "    data = [[picks_raw.iloc[i]['fname'].split(\".\")[0], starttime + pick_time, \"s\", pick_prob] for pick_time, pick_prob in zip(new_pick_times, new_pick_probs)] # Station, pick time (UTC), phase, probability\n",
        "    picks = picks.append(pd.DataFrame(data, columns=['sta', 'time', 'pha', 'prob']), ignore_index=True) # Append to picks dataframe\n",
        "  \n",
        "  # Track progress:\n",
        "  tenth = math.floor(picks_raw.shape[0] / 10)\n",
        "  if i % tenth == 0:\n",
        "    print(str(round(10 * i / tenth))  + \"% done\")\n",
        "  if i == (picks_raw.shape[0] - 1):\n",
        "    print(\"100% done\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "SOLFf1_qvfaU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Print new picks dataframe to screen\n",
        "picks"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZjCi1L3Atsub",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Sort this new dataframe by pick time and print to screen again\n",
        "picks.sort_values(by='time', inplace=True, ignore_index=True)\n",
        "picks"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ip_g7wYKyycU",
        "colab_type": "text"
      },
      "source": [
        "This format looks a little more useful in terms of producing a catalogue or creating a pick file for another program (e.g., Hypoinverse or NonLinLoc).\n",
        "\n",
        "We can see in the third and fourth rows above that we get two P-wave arrival picks really close together (0.07 seconds apart). This is because we have generated a pick from each of the two passes over the data. It's clear that sometimes the two passes will not give the exact same pick time, as above, and the probability of a phase arrival from each pass is also different (0.929 vs 0.950 above). Pretty annoying and makes it hard to choose just one pick!\n",
        "\n",
        "One strategy here could be to remove any picks that are within, say, 0.1 seconds of each other, taking the one with the highest prob (the more confident pick). Let's loop through our picks dataframe again and do this:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PcI9ElTbzPdb",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "picks_final = pd.DataFrame() # New empty dataframe to accommodate picks\n",
        "stations = picks.sta.unique() # List of station names\n",
        "phases = [\"p\", \"s\"] # List of phase names\n",
        "# Loop over stations and phases\n",
        "for sta in stations:\n",
        "  for pha in phases:\n",
        "    this_pha = picks.loc[(picks.pha == pha) & (picks.sta == sta)].reset_index(drop=True) # Extract dataframe for each phase and station\n",
        "    this_pha['diff_time'] = this_pha.time.diff() # Find difference in times between picks across two passes\n",
        "    for i in range(this_pha.shape[0]): # For each pick\n",
        "      if this_pha.diff_time[i] < 0.1: # If difference in time is less than 0.1 secs\n",
        "        this_pha.drop(i - np.argmax((this_pha.prob[(i-1)], this_pha.prob[i])), inplace=True) # Drop the one with the lower probability\n",
        "    this_pha = this_pha.reset_index(drop=True) # At the end, reset index column of dataframe (not during the loop or iteration won't work properly)\n",
        "    # Recombine to make final pick df\n",
        "    picks_final = picks_final.append(this_pha, ignore_index=True) # Append all df's from loop\n",
        "# Sort final picks df by time\n",
        "picks_final.sort_values(by=\"time\", inplace=True, ignore_index=True) # Sort by time\n",
        "picks_final.drop('diff_time', axis=1, inplace=True) # Drop diff_time column"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ET6O6MYN0uSh",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Print new dataframe to screen\n",
        "picks_final"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hPaPFjuz8AUV",
        "colab_type": "text"
      },
      "source": [
        "Cool, kind of a lengthy process but we now have some pick times! The default probability threshold for a pick is 0.3.... Now that we have the pick data in the format above, if we wanted to reduce the no. of potential false picks, at the cost of potentially losing some true picks, we could increase the threshold to, say, 0.8:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "C5LMQsRF078W",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "picks_final.loc[picks_final.prob >= 0.8]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ugd4T2tZ8cgO",
        "colab_type": "text"
      },
      "source": [
        "This reduces the number of picks from a total of 3689 across the two stations to 1441, which still sounds like a lot. Let's plot a section of data with picks using the original probability threshold of 0.3.\n",
        "\n",
        "# 10. Plot some picks over waveforms\n",
        "\n",
        "As far as I'm away, obspy doesn't have an option for adding pick times to plots. One crude way is to just plot the waveform data with `pyplot` and then add vertical lines for our picks at the right index values.\n",
        "\n",
        "Below we plot the first hour of data with P- and S-phase picks."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "mL2YqspL18f_",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "CCC_cut = CCC.copy() # Copy station CCC original mseed (7 hours data)\n",
        "CCC_cut.trim(starttime, starttime + (60 * 60)) # Trim just the first hour"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hwI50zb48lvU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet P-phase picks for first hour of data\")\n",
        "plt.plot(CCC_cut[2].data, 'k') # Vertical component data only\n",
        "plt.vlines(((picks_final.loc[(picks_final.sta == \"CCC\") & (picks_final.pha == \"p\") & (picks_final.prob >= 0.3)].time  - starttime) / ccc_delta).values, ymin=-1e7, ymax=1e7, colors='c', linestyles='dashed')\n",
        "plt.xlim(0, len(CCC_cut[2].data))\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "60VKFKUk95U5",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet S-phase picks for first hour of data\")\n",
        "plt.plot(CCC_cut[2].data, 'k')\n",
        "plt.vlines(((picks_final.loc[(picks_final.sta == \"CCC\") & (picks_final.pha == \"s\") & (picks_final.prob >= 0.3)].time  - starttime) / ccc_delta).values, ymin=-1e7, ymax=1e7, colors='m', linestyles='dashed')\n",
        "plt.xlim(0, len(CCC_cut[2].data))\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MzxlW_YzfF99",
        "colab_type": "text"
      },
      "source": [
        "We can see an increase in the number of phases picked after the main foreshock a little over halfway through. Let's zoom on the first few mins, before the foreshock hit, to see how well it picks up smaller events."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_h8LRTFR8mec",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "CCC_cut = CCC.copy() # Copy station CCC original mseed (7 hours data)\n",
        "CCC_cut.trim(starttime, starttime + (5 * 60)) # Trim just the first 5 mins\n",
        "CCC_cut.detrend()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YnCHyV6VOWKE",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet P-phase picks for first 5 mins of data\")\n",
        "plt.plot(CCC_cut[2].data, 'k')\n",
        "trace_max = max(abs(CCC_cut[2].data))\n",
        "plt.vlines(((picks_final.loc[(picks_final.sta == \"CCC\") & (picks_final.pha == \"p\") & (picks_final.prob >= 0.3)].time  - starttime) / ccc_delta).values, ymin=-trace_max, ymax=trace_max, colors='c', linestyles='dashed')\n",
        "plt.xlim(0, len(CCC_cut[2].data))\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BKfLPCF3OZ7a",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet S-phase picks for first 5 mins of data\")\n",
        "plt.plot(CCC_cut[2].data, 'k')\n",
        "plt.vlines(((picks_final.loc[(picks_final.sta == \"CCC\") & (picks_final.pha == \"s\") & (picks_final.prob >= 0.3)].time  - starttime) / ccc_delta).values, ymin=-trace_max, ymax=trace_max, colors='m', linestyles='dashed')\n",
        "plt.xlim(0, len(CCC_cut[2].data))\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Fxi3l3XviBv4",
        "colab_type": "text"
      },
      "source": [
        "Looks good on that one event. Let's just check that first P- and S-phase pick..."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8vpYR403iKbJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet P-phase picks between 2 mins 40 secs in and 3 mins in - total of 20 secs\")\n",
        "plt.plot(CCC_cut[2].data[0:18000], 'k')\n",
        "trace_max = max(abs(CCC_cut[2].data[16000:18000]))\n",
        "plt.vlines(((picks_final.loc[(picks_final.sta == \"CCC\") & (picks_final.pha == \"p\") & (picks_final.prob >= 0.3)].time  - starttime) / ccc_delta).values, ymin=-trace_max, ymax=trace_max, colors='c', linestyles='dashed')\n",
        "plt.xlim(16000, 18000)\n",
        "plt.ylim(-600, 0)\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "v5O5ZVZpj6Zw",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet S-phase picks between 2 mins 40 secs in and 3 mins in - total of 20 secs\")\n",
        "plt.plot(CCC_cut[1].data[0:18000], 'k')\n",
        "trace_max = max(abs(CCC_cut[1].data[16000:18000]))\n",
        "plt.vlines(((picks_final.loc[(picks_final.sta == \"CCC\") & (picks_final.pha == \"s\") & (picks_final.prob >= 0.3)].time  - starttime) / ccc_delta).values, ymin=-trace_max, ymax=trace_max, colors='m', linestyles='dashed')\n",
        "plt.xlim(16000, 18000)\n",
        "plt.ylim(-1500,1000)\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rmWfnsBbpZho",
        "colab_type": "text"
      },
      "source": [
        "Can't be 100% certain on that last example, but generally looks quite good!\n",
        "\n",
        "# 11. Magnitude Estimation - can't do!\n",
        "\n",
        "I had hoped we could run MagNet here, the local/duration magnitude estimation model produced by the same authors, but that is unfortunately not available on their Github at the mo, just their results (https://github.com/smousavi05/MagNet). We could produce the magnitude estimation model from scratch, following the architecture and methods in their paper, but the authors only provide the raw dataset (https://github.com/smousavi05/STEAD), which is ~ 90 GB in size and would involve investing a decent amount of time in pre-processing to produce examples and labels for model training. Could be a practical for the future though?\n",
        "\n",
        "Instead, lets try the above model on 'out-of-sample' data, i.e., not included in model training and not from California.\n",
        "\n",
        "# 12. Out-of-sample examples\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ki9Eek8xO9lq",
        "colab_type": "text"
      },
      "source": [
        "First, let's define a function to do all that conversion from PhaseNet output to a nice dataframe of UTC pick times and probabilities."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RZDpAYu4Od49",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Function to read in PhaseNet csv outputs and convert to a more useful pandas df\n",
        "\n",
        "def PNcsv2df(csv_output, starttime, npts, delta, prune = True, difftime_thresh = 0.1):\n",
        "  \"\"\"\n",
        "  Inputs:\n",
        "  - csv_output (csv file) = csv output file from PhaseNet run.py script\n",
        "  - starttime (obspy UTC time) = starttime for 3-component trace inputted to PhaseNet\n",
        "  - npts (int?) = npts for 3-component trace inputted to PhaseNet\n",
        "  - delta (float) = delta for 3-component trace inputted to PhaseNet\n",
        "  - difftime (float) = prune picks within difftime seconds of each other\n",
        "  \"\"\"\n",
        "  picks_raw = pd.read_csv(csv_output)\n",
        "  if picks_raw.shape[0] < 1:\n",
        "    return(\"Error: no picks in PhaseNet csv file\")\n",
        "  \n",
        "  # Empty pick dataframe\n",
        "  picks = pd.DataFrame(columns=['sta', 'time', 'pha', 'prob'])\n",
        "\n",
        "  for i in range(picks_raw.shape[0]):\n",
        "    abs_index = int(picks_raw.iloc[i]['fname'].split(\"_\")[-1]) # Get absolute start index value\n",
        "    which_pass = 1 if abs_index >= npts else 0 # 0 if first pass, 1 if second pass\n",
        "    if picks_raw.iloc[i]['itp'] != '[]':\n",
        "      new_start_time = (abs_index % npts) * delta # Start time of mseed file (secs)\n",
        "      new_pick_times = [new_start_time + (t * delta) - (which_pass * 15) for t in list(map(float, picks_raw.iloc[i]['itp'][1:-1].split()))] # P-phase pick time in secs\n",
        "      new_pick_probs = list(map(float, picks_raw.iloc[i]['tp_prob'][1:-1].split())) # P-phase pick probability\n",
        "      data = [[picks_raw.iloc[i]['fname'].split(\".\")[0], starttime + pick_time, \"p\", pick_prob] for pick_time, pick_prob in zip(new_pick_times, new_pick_probs)]  # List of station name, pick time (UTC), phase type (p/s), probability\n",
        "      picks = picks.append(pd.DataFrame(data, columns=['sta', 'time', 'pha', 'prob']), ignore_index=True) # Append to dataframe\n",
        "    if picks_raw.iloc[i]['its'] != '[]':\n",
        "      # Same as above but for S-phase\n",
        "      new_start_time = (abs_index % npts) * delta\n",
        "      new_pick_times = [new_start_time + (t * delta) - (which_pass * 15) for t in list(map(float, picks_raw.iloc[i]['its'][1:-1].split()))]\n",
        "      new_pick_probs = list(map(float, picks_raw.iloc[i]['ts_prob'][1:-1].split()))\n",
        "      data = [[picks_raw.iloc[i]['fname'].split(\".\")[0], starttime + pick_time, \"s\", pick_prob] for pick_time, pick_prob in zip(new_pick_times, new_pick_probs)] \n",
        "      picks = picks.append(pd.DataFrame(data, columns=['sta', 'time', 'pha', 'prob']), ignore_index=True)\n",
        "    \n",
        "    # Track progress:\n",
        "    tenth = math.floor(picks_raw.shape[0] / 10)\n",
        "    if i == 0:\n",
        "      print(\"Converting raw csv file to more useful pick table\")\n",
        "    if i % tenth == 0:\n",
        "      print(str(round(10 * i / tenth))  + \"% done\")\n",
        "    if i == (picks_raw.shape[0] - 1):\n",
        "      print(\"100% done\")\n",
        "  \n",
        "  # Sort by pick time (any phase or station)\n",
        "  picks.sort_values(by='time', inplace=True, ignore_index=True)\n",
        "\n",
        "  # Prune double picks from two pass approach and return pruned dataframe:\n",
        "  if prune == True:\n",
        "    picks_final = pd.DataFrame()\n",
        "    print(\"Pruning picks within \" + str(difftime_thresh) + \" seconds of each other\")\n",
        "    stations = picks.sta.unique()\n",
        "    phases = [\"p\", \"s\"]\n",
        "    for s in stations:\n",
        "      for p in phases:\n",
        "        this_pha = picks.loc[(picks.pha == p) & (picks.sta == s)].reset_index(drop=True) # Extract dataframe for each phase and station\n",
        "        this_pha['diff_time'] = this_pha.time.diff() # Find difference in times between picks across two passes\n",
        "        for i in range(this_pha.shape[0]): # For each pick\n",
        "          if this_pha.diff_time[i] < difftime_thresh: # If difference in time is less than difftime_thresh sec\n",
        "            this_pha.drop(i - np.argmax((this_pha.prob[(i-1)], this_pha.prob[i])), inplace=True) # Drop the one with the lower probability\n",
        "        this_pha = this_pha.reset_index(drop=True) # At the end, reset index column of dataframe (not during the loop or iteration won't work properly)\n",
        "        # Recombine to make final pick df\n",
        "        picks_final = picks_final.append(this_pha, ignore_index=True) # Append all df's from loop\n",
        "    # Sort final picks df by time\n",
        "    picks_final.sort_values(by=\"time\", inplace=True, ignore_index=True) # Sort by time\n",
        "    picks_final.drop('diff_time', axis=1, inplace=True) # Drop diff_time column\n",
        "    return(picks_final) # Return pruned df\n",
        "\n",
        "  # Else return unpruned picks df:\n",
        "  else:\n",
        "    return(picks)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NrwnQGVAvUrk",
        "colab_type": "text"
      },
      "source": [
        "Now let's get some data from somewhere else... \n",
        "\n",
        "## 12.1 Costa Rica 6.5 Mw event @ 2017-11-13 02:28:23"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fRTdnZ3reUkk",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Let's get some data from Costa Rica\n",
        "\n",
        "# Start client for requesting data from IRIS\n",
        "client = Client(\"IRIS\")\n",
        "\n",
        "# Create directory to store data if needed\n",
        "data_dir = \"mseed_cr\"\n",
        "if not os.path.exists(data_dir):\n",
        "  os.makedirs(data_dir)\n",
        "\n",
        "# Give start and end time to extract data - this was date of a mag 6.5 eq\n",
        "starttime = UTCDateTime(\"2017-11-13T00:00:00\")\n",
        "endtime = UTCDateTime(\"2017-11-13T08:00:00\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OpGD2CMkecU8",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Get waveforms\n",
        "LAFE = client.get_waveforms(\"OV\", \"LAFE\", \"*\", \"HHE,HHN,HHZ\", starttime, endtime)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WzgskZAujs38",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "LAFE.plot()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WXT8EVmGmWEt",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "LAFE.write(os.path.join(data_dir, \"LAFE.mseed\")) # Write LAFE station data to LAFE.mseed"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yJgHWly1nDL9",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Write file list\n",
        "with open(\"fname_cr.csv\", 'w') as fp:\n",
        "  fp.write(\"fname,E,N,Z\\n\")\n",
        "  fp.write(\"LAFE.mseed,HHE,HHN,HHZ\\n\")"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "-SlHfIDhnTyk",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Check file list csv\n",
        "!cat fname_cr.csv"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2uGl0hI1nL2Y",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Run PhaseNet on above data\n",
        "!python PhaseNet/run.py --mode=pred --model_dir=PhaseNet/model/190703-214543 --data_dir=mseed_cr --data_list=fname_cr.csv --output_dir=output_cr --batch_size=20 --input_mseed"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1rWoZFjQnvgO",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Convert output from PhaseNet to nicer-looking pick dataframe:\n",
        "LAFE_picks = PNcsv2df(csv_output = \"output_cr/picks.csv\", \n",
        "                      starttime = LAFE.traces[0].stats['starttime'], \n",
        "                      npts = LAFE.traces[0].stats['npts'], \n",
        "                      delta = LAFE.traces[0].stats['delta'], \n",
        "                      prune = True, \n",
        "                      difftime_thresh = 0.1)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EoCal_nsqm3l",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "LAFE_picks"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qaFBRv2ircc7",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Look at first two hours 25 mins:\n",
        "LAFE_cut = LAFE.copy()\n",
        "LAFE_cut.trim(starttime, starttime + (120 * 60)) # Look at first 2 hours 25 mins (couple mins before mag 6.5 event)\n",
        "LAFE_cut.filter('highpass', freq=2.0)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "fFGycJ5mthkU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.plot(LAFE_cut[2].data, 'k')\n",
        "plt.vlines(((LAFE_picks.loc[(LAFE_picks.sta == \"LAFE\") & (LAFE_picks.pha == \"p\") & (picks_final.prob >= 0.3)].time  - starttime) / LAFE_cut.traces[2].stats['delta']).values, ymin=-1e4, ymax=1e4, colors='c', linestyles='dashed')\n",
        "plt.xlim(0, len(LAFE_cut[2].data))\n",
        "plt.ylim(-8000,8000)\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "y_kDTYzZt0jg",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.plot(LAFE_cut[0].data, 'k')\n",
        "plt.vlines(((LAFE_picks.loc[(LAFE_picks.sta == \"LAFE\") & (LAFE_picks.pha == \"s\") & (picks_final.prob >= 0.3)].time  - starttime) / LAFE_cut.traces[2].stats['delta']).values, ymin=-1e4, ymax=1e4, colors='m', linestyles='dashed')\n",
        "plt.xlim(0, len(LAFE_cut[0].data))\n",
        "plt.ylim(-8000,8000)\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "sO3KuQNZxnzj",
        "colab_type": "text"
      },
      "source": [
        "Looks like it might have missed some potential events there, but maybe they're just bits of noise... let's look at the data when the mag 6.5 event hits."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HH_MDUWguAew",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Look at third hour:\n",
        "LAFE_cut = LAFE.copy()\n",
        "LAFE_cut.trim(starttime + (120 * 60), starttime + (180 * 60)) # Trim just the third hour\n",
        "LAFE_cut.filter('highpass', freq=2.0)\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EMqKgnFXu3GN",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.plot(LAFE_cut[2].data, 'k')\n",
        "trace_max = max(abs(LAFE_cut[2].data))\n",
        "plt.vlines(((LAFE_picks.loc[(LAFE_picks.sta == \"LAFE\") & (LAFE_picks.pha == \"p\") & (picks_final.prob >= 0.3)].time  - (starttime + (120*60))) / LAFE_cut.traces[2].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='c', linestyles='dashed')\n",
        "plt.xlim(0, len(LAFE_cut[2].data))\n",
        "plt.ylim(-1e5,1e5)\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6JWBo2VevZtR",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.plot(LAFE_cut[2].data, 'k')\n",
        "trace_max = max(abs(LAFE_cut[2].data))\n",
        "plt.vlines(((LAFE_picks.loc[(LAFE_picks.sta == \"LAFE\") & (LAFE_picks.pha == \"s\") & (picks_final.prob >= 0.3)].time  - (starttime + (120*60))) / LAFE_cut.traces[2].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='m', linestyles='dashed')\n",
        "plt.xlim(0, len(LAFE_cut[2].data))\n",
        "plt.ylim(-1e5,1e5)\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "0WsYbqYUxiPR",
        "colab_type": "text"
      },
      "source": [
        "\n",
        "## 12.2 LASSO experiment\n",
        "\n",
        "Many of you work with data that might be vertical component only or have different sample rates to 100 Hz. Let's look at some data from the LASSO experiment (https://doi.org/10.1785/0220190094), where > 1800 vertical component nodal seismometers were deployed in Oklahoma. This is available on IRIS PH5 under station 2A. The data is recorded at 500 samples/s (the data for training PhaseNet had a variety of sample rates but were all 3-component).\n",
        "\n",
        "\n",
        "First, we can do an IRIS event data search to see what might be a good section of data to look at. Here we run a search using the IRIS DMC FDSN event web service, where the URL specifies search arguments:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "atc37xwMObjf",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "!wget -O- http://service.iris.edu/fdsnws/event/1/query\\?starttime\\=2016-04-15\\&endtime\\=2016-05-10\\&minlatitude\\=36.5\\&maxlatitude\\=37\\&minlongitude\\=-98.2\\&maxlongitude\\=-97.7\\&minmagnitude\\=0\\&includeallmagnitudes\\=true\\&orderby\\=time\\&format\\=text"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "r-thp9MHQeMk",
        "colab_type": "text"
      },
      "source": [
        "The three events on 2016-05-04 look to occur really close together so let's get data for that day. This data is only vertical component, and PhaseNet requires three-component input, but lets try it anyway to see what happens (we will use vertical component data for all 3 channels)."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9wJMoY33yJzU",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Start client for requesting data from IRIS - slightly different process for requesting PH5 data\n",
        "DATASELECT = 'http://service.iris.edu/ph5ws/dataselect/1'\n",
        "client = Client('http://service.iris.edu',\n",
        "                       service_mappings={\n",
        "                           'dataselect': DATASELECT\n",
        "                       },\n",
        "                       debug=True\n",
        "                      )\n",
        "\n",
        "# Create directory to store data if needed\n",
        "data_dir = \"mseed_lasso\"\n",
        "if not os.path.exists(data_dir):\n",
        "  os.makedirs(data_dir)\n",
        "\n",
        "# Give start and end time to extract data\n",
        "starttime = UTCDateTime(\"2016-05-04T05:00:00\") # 5am\n",
        "endtime = UTCDateTime(\"2016-05-04T14:00:00\") # to 2pm\n",
        "\n",
        "# Get waveforms\n",
        "L398 = client.get_waveforms(\"2A\", \"398\", \"*\", \"DPZ\", starttime, endtime)\n",
        "\n",
        "# Write data to mseed\n",
        "L398.write(os.path.join(data_dir, \"L398.mseed\"))\n",
        "\n",
        "# Write file list\n",
        "with open(\"fname_lasso.csv\", 'w') as fp:\n",
        "  fp.write(\"fname,E,N,Z\\n\")\n",
        "  fp.write(\"L398.mseed,DPZ,DPZ,DPZ\\n\")\n",
        "\n",
        "# Check file list csv\n",
        "!cat fname_lasso.csv"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "QDPo1XxzXc-u",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "L398.plot()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Vb-Eps3rVfoB",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Run PhaseNet on above data\n",
        "!python PhaseNet/run.py --mode=pred --model_dir=PhaseNet/model/190703-214543 --data_dir=mseed_lasso --data_list=fname_lasso.csv --output_dir=output_lasso --batch_size=100 --input_mseed"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BQnHq8fIXYcM",
        "colab_type": "text"
      },
      "source": [
        "Cool, seems to have worked. Let's look at picks:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "P_HB3yheW-af",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Convert output from PhaseNet output to nicer-looking pick dataframe:\n",
        "L398_picks = PNcsv2df(csv_output = \"output_lasso/picks.csv\", \n",
        "                      starttime = L398.traces[0].stats['starttime'], \n",
        "                      npts = L398.traces[0].stats['npts'], \n",
        "                      delta = L398.traces[0].stats['delta'], \n",
        "                      prune = True, \n",
        "                      difftime_thresh = 0.05)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "TcrW1vFBYRN9",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "L398_picks"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Qk9AmISgYXvs",
        "colab_type": "text"
      },
      "source": [
        "So we have 387 P and S phase picks, but are they any good??\n",
        "\n",
        "We seem to get some straight away, at around 5am, so let's just look at the first few mins:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "p00Vu3azYUoF",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Look at first few mins:\n",
        "L398_cut = L398.copy()\n",
        "L398_cut.trim(starttime, starttime + (5 * 60)) # Trim just the first 5 mins"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "0rUzIu79Yyej",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "L398_cut.plot()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5tV_8Nm5Y0aX",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet P-phase picks for first 5 mins of 500 sample/s data\")\n",
        "plt.plot(L398_cut[0].data, 'k')\n",
        "trace_max = max(abs(L398_cut[0].data))\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"p\")].time  - starttime) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='c', linestyles='dashed')\n",
        "plt.xlim(0, len(L398_cut[0].data))\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ha9EQOXCZn9A",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet P-phase picks for first 5 mins of 500 sample/s data\")\n",
        "plt.plot(L398_cut[0].data, 'k')\n",
        "trace_max = max(abs(L398_cut[0].data))\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"s\")].time  - starttime) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='m', linestyles='dashed')\n",
        "plt.xlim(0, len(L398_cut[0].data))\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kollT3Lraklw",
        "colab_type": "text"
      },
      "source": [
        "Not great, but looks like it might have picked up that S-wave:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "WdYzQbkJZuU7",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet S-phase picks zoomed-in\")\n",
        "plt.plot(L398_cut[0].data, 'k')\n",
        "trace_max = max(abs(L398_cut[0].data))\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"s\")].time  - starttime) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='m', linestyles='dashed')\n",
        "plt.xlim(13000,18000)\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eFhyGYkabFNp",
        "colab_type": "text"
      },
      "source": [
        "Yeah.... maybe... What about those S-wave picks just slightly before?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jJT50ePQazPF",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet S-phase picks zoomed-in\")\n",
        "plt.plot(L398_cut[0].data, 'k')\n",
        "trace_max = max(abs(L398_cut[0].data[8000:13000]))\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"s\")].time  - starttime) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='m', linestyles='dashed')\n",
        "plt.xlim(8000,13000)\n",
        "plt.ylim(-trace_max,trace_max)\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4rEnfejTBGZu",
        "colab_type": "text"
      },
      "source": [
        "Less convinced by these picks. That third one has a probability of 0.768472 of being an S-phase arrival. So maybe this approach doesn't work... understandable as there should obviously be differences between horizontal and vertical component data. We could maybe rotate the vertical component to two orthogonal horizontal components.\n",
        "\n",
        "Let's look at the three events from the IRIS event data as these are larger magnitude:\n",
        "\n",
        "```\n",
        "5187310|2016-05-04T13:18:01|36.8708|-98.1407|7.033|us|NEIC PDE|us|us10005dli|ML|2.6|tul|OKLAHOMA\n",
        "5187308|2016-05-04T12:55:20|36.8675|-98.1341|5.342|us|NEIC PDE|us|us10005dl7|ML|2.7|tul|OKLAHOMA\n",
        "5181718|2016-05-04T12:20:25|36.8731|-98.1372|5.0|us|NEIC PDE|us|us10005dl1|ML|3.3|tul|OKLAHOMA\n",
        "```"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "pleRZOx2A1dn",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "L398_cut = L398.copy()\n",
        "L398_cut.trim(UTCDateTime(\"2016-05-04T12:15:00\"), UTCDateTime(\"2016-05-04T13:30:00\")) # 75 min window around three events above"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ult8lIKMCQ5W",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "L398_cut.plot()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "P8gP2eX0CXSX",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "events = [UTCDateTime(\"2016-05-04T12:20:25\"), UTCDateTime(\"2016-05-04T12:55:20\"), UTCDateTime(\"2016-05-04T13:18:01\")]\n",
        "#[(ev - UTCDateTime(\"2016-05-04T12:15:00\")) / L398.traces[0].stats['delta'] for ev in events]\n",
        "\n",
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet P-phase picks between 12:15 and 13:30 - three event times in red\")\n",
        "plt.plot(L398_cut[0].data, 'k')\n",
        "trace_max = max(abs(L398_cut[0].data))\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"p\")].time  - UTCDateTime(\"2016-05-04T12:15:00\")) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='c', linestyles='dashed')\n",
        "plt.vlines(([(ev - UTCDateTime(\"2016-05-04T12:15:00\")) / L398.traces[0].stats['delta'] for ev in events]), ymin=-trace_max, ymax=trace_max, colors='r', linewidth=4)\n",
        "plt.xlim(0, len(L398_cut[0].data))\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "3hpvOLLICklJ",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet S-phase picks between 12:15 and 13:30 - three event times in red\")\n",
        "plt.plot(L398_cut[0].data, 'k')\n",
        "trace_max = max(abs(L398_cut[0].data))\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"s\")].time  - UTCDateTime(\"2016-05-04T12:15:00\")) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='m', linestyles='dashed')\n",
        "plt.vlines(([(ev - UTCDateTime(\"2016-05-04T12:15:00\")) / L398.traces[0].stats['delta'] for ev in events]), ymin=-trace_max, ymax=trace_max, colors='r', linewidth=4)\n",
        "plt.xlim(0, len(L398_cut[0].data))\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "oDjlIPh1Rem_",
        "colab_type": "text"
      },
      "source": [
        "Let's just look at picks with probability >= 0.8."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "tKOm-OQsm_Nq",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet P-phase picks between 12:15 and 13:30 - three event times in red\")\n",
        "plt.plot(L398_cut[0].data, 'k')\n",
        "trace_max = max(abs(L398_cut[0].data))\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"p\") & (L398_picks.prob >= 0.8)].time  - UTCDateTime(\"2016-05-04T12:15:00\")) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='c', linestyles='dashed')\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"s\") & (L398_picks.prob >= 0.8)].time  - UTCDateTime(\"2016-05-04T12:15:00\")) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='m', linestyles='dashed')\n",
        "plt.vlines(([(ev - UTCDateTime(\"2016-05-04T12:15:00\")) / L398.traces[0].stats['delta'] for ev in events]), ymin=-trace_max, ymax=trace_max, colors='r', linewidth=4)\n",
        "plt.xlim(0, len(L398_cut[0].data))\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "k3YDvc4hf1eu",
        "colab_type": "text"
      },
      "source": [
        "Let's zoom in on each event individually..."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "8vAmhasGCvWu",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "L398_cut = L398.copy()\n",
        "L398_cut.trim(events[0]-30, events[0]+30) # 30 secs either side of event\n",
        "\n",
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet P- and S-phase picks around event 1\")\n",
        "plt.plot(L398_cut[0].data, 'k')\n",
        "trace_max = max(abs(L398_cut[0].data))\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"p\")].time  - (events[0] - 30)) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='c', linestyles='dashed')\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"s\")].time  - (events[0] - 30)) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='m', linestyles='dashed')\n",
        "plt.xlim(0, len(L398_cut[0].data))\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ht87Srn4gRzD",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "L398_cut = L398.copy()\n",
        "L398_cut.trim(events[1]-30, events[1]+30) # 30 secs either side of event\n",
        "\n",
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet P- and S-phase picks around event 2\")\n",
        "plt.plot(L398_cut[0].data, 'k')\n",
        "trace_max = max(abs(L398_cut[0].data))\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"p\")].time  - (events[1] - 30)) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='c', linestyles='dashed')\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"s\")].time  - (events[1] - 30)) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='m', linestyles='dashed')\n",
        "plt.xlim(0, len(L398_cut[0].data))\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "O12_7mbFgb39",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "L398_cut = L398.copy()\n",
        "L398_cut.trim(events[2]-30, events[2]+30) # 30 secs either side of event\n",
        "\n",
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet P- and S-phase picks around event 3\")\n",
        "plt.plot(L398_cut[0].data, 'k')\n",
        "trace_max = max(abs(L398_cut[0].data))\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"p\")].time  - (events[2] - 30)) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='c', linestyles='dashed')\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"s\")].time  - (events[2] - 30)) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='m', linestyles='dashed')\n",
        "plt.xlim(0, len(L398_cut[0].data))\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DI9BgY51g-CR",
        "colab_type": "text"
      },
      "source": [
        "Seems to pick P-arrival pretty well but no S-'s, I guess because the S-wave is not clear in that vertical component trace and that's what we've used as our horizontal components.\n",
        "\n",
        "I'm intrigued by those 4 P-wave picks just before the event... let's zoom in again."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YC3bxyXGg0uE",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "L398_cut = L398.copy()\n",
        "L398_cut.trim(events[2]-30, events[2]) # 30 secs either side of event\n",
        "\n",
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet P- and S-phase picks before event 3\")\n",
        "plt.plot(L398_cut[0].data, 'k')\n",
        "trace_max = max(abs(L398_cut[0].data))\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"p\")].time  - (events[2] - 30)) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='c', linestyles='dashed')\n",
        "plt.vlines(((L398_picks.loc[(L398_picks.pha == \"s\")].time  - (events[2] - 30)) / L398.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='m', linestyles='dashed')\n",
        "plt.xlim(0, len(L398_cut[0].data))\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "VqCqitUyhXym",
        "colab_type": "text"
      },
      "source": [
        "I think they're all just false picks.\n",
        "\n",
        "How about some actual 3-component data with higher sample rate?"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4P5HSfQqSIqz",
        "colab_type": "text"
      },
      "source": [
        "# 12.3 Peanut Lake\n",
        "\n",
        "I've picked a random dataset from IRIS that was recorded on 3-component seismometers with 500 Hz sample rate. The site is called \"Peanut Lake\", somewhere in Antarctica, and the study is called the *Impact of Supraglacial Lakes on Ice-Shelf Stability* (https://doi.org/10.7914/SN/9G_2016). Not sure what we'll find (hopefully some ice quakes!) but wanted some actual 3-component high-sample-rate data to try PhaseNet on. The data here will unlikely look like anything in the PhaseNet training data.\n",
        "\n",
        "We'll run the same steps again, requesting the data from IRIS, saving to mseed, creating a file list and running the PhaseNet `run.py` script. We'll then convert the picks to a nice looking dataframe and plot some of them."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "x6rLO2JRTcjM",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Start client for requesting data from IRIS\n",
        "client = Client(\"IRIS\")\n",
        "\n",
        "# Create directory to store data if needed\n",
        "data_dir = \"mseed_peanut\"\n",
        "if not os.path.exists(data_dir):\n",
        "  os.makedirs(data_dir)\n",
        "\n",
        "# Give start and end time to extract data - picked a random day during deployment (2016-11-20 to 2017-01-20)\n",
        "starttime = UTCDateTime(\"2016-12-25T06:00:00\")\n",
        "endtime = UTCDateTime(\"2016-12-25T12:00:00\")\n",
        "\n",
        "# Get waveforms\n",
        "SP01 = client.get_waveforms(\"9G\", \"SP01\", \"*\", \"DHE,DHN,DHZ\", starttime, endtime)\n",
        "\n",
        "# Write data to mseed\n",
        "SP01.write(os.path.join(data_dir, \"SP01.mseed\"))\n",
        "\n",
        "# Write file list\n",
        "with open(\"fname_peanut.csv\", 'w') as fp:\n",
        "  fp.write(\"fname,E,N,Z\\n\")\n",
        "  fp.write(\"SP01.mseed,DHE,DHN,DHZ\\n\")\n",
        "\n",
        "# # Check file list csv\n",
        "!cat fname_peanut.csv"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "2WaMkOqSUTGF",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "SP01.plot()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "uDHPolXxUnTu",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Run PhaseNet on above data\n",
        "!python PhaseNet/run.py --mode=pred --model_dir=PhaseNet/model/190703-214543 --data_dir=mseed_peanut --data_list=fname_peanut.csv --output_dir=output_peanut --batch_size=100 --input_mseed"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jYLDgWgXU7VL",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "# Convert output from PhaseNet output to nicer-looking pick dataframe:\n",
        "SP01_picks = PNcsv2df(csv_output = \"output_peanut/picks.csv\", \n",
        "                      starttime = SP01.traces[0].stats['starttime'], \n",
        "                      npts = SP01.traces[0].stats['npts'], \n",
        "                      delta = SP01.traces[0].stats['delta'], \n",
        "                      prune = True, \n",
        "                      difftime_thresh = 0.05)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LNu9coKJVWWI",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "SP01_picks"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tLJYocrIVcKN",
        "colab_type": "text"
      },
      "source": [
        "426 picks. Could be good! Let's just look at ones with a threshold above 0.9:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9-8HhK3bVac0",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "SP01_picks.loc[SP01_picks['prob'] >= 0.9]"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "344pDQlNVrsI",
        "colab_type": "text"
      },
      "source": [
        "We get 9 picks with pretty high probability. Let's look at the one with the highest probability (0.982):"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "jpQFyabdVnkn",
        "colab_type": "code",
        "colab": {}
      },
      "source": [
        "SP01_cut = SP01.copy()\n",
        "SP01_cut.trim(UTCDateTime(\"2016-12-25T11:06:12\"), UTCDateTime(\"2016-12-25T11:06:22\"))\n",
        "SP01_cut.detrend()\n",
        "\n",
        "plt.rcParams[\"figure.figsize\"] = (20,6)\n",
        "plt.title(\"PhaseNet P- and S-phase picks from Peanut Lake\")\n",
        "plt.plot(SP01_cut[0].data, 'k')\n",
        "trace_max = max(abs(SP01_cut[0].data))\n",
        "plt.vlines(((SP01_picks.loc[(SP01_picks.pha == \"p\") & (SP01_picks.prob >= 0.97)].time  - UTCDateTime(\"2016-12-25T11:06:12\")) / SP01.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='c', linestyles='dashed')\n",
        "plt.vlines(((SP01_picks.loc[(SP01_picks.pha == \"s\") & (SP01_picks.prob >= 0.97)].time  - UTCDateTime(\"2016-12-25T11:06:12\")) / SP01.traces[0].stats['delta']).values, ymin=-trace_max, ymax=trace_max, colors='m', linestyles='dashed')\n",
        "plt.xlim(0, len(SP01_cut[0].data))\n",
        "plt.show()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zJpQ4xeDYL-3",
        "colab_type": "text"
      },
      "source": [
        "Okay, so that's not a promising pick."
      ]
    }
  ]
}