<jupyter_start><jupyter_text>Running U-GPD Transfer Learning Model on Data from Nabro VolcanoLast updated 28th Jul 2021.___The code in this notebook accompanies the JGR: Solid Earth paper **_"Lapins et al. (2021) A Little Data goes a Long Way: Automating Seismic Phase Arrival Picking at Nabro Volcano with Transfer Learning"_**, https://doi.org/10.1029/2021JB021910 (open access).If you find the code, data or results in this notebook useful, please consider citing the above paper or the dataset directly (https://doi.org/10.5281/zenodo.4498549).This notebook has been developed and tested in Google Colab. A conda environment yml file (https://github.com/sachalapins/U-GPD/ugpd.yml) is provided for those who wish to download this notebook and run it locally but note that small changes may need to be made.If you have any suggestions for improvements, issues running the code or notice any errors then please feel free to contact me ([<EMAIL>](mailto:<EMAIL>)).Many thanks.<NAME>, University of Bristol---   In this notebook, we demonstrate how to run our U-GPD transfer learning model, trained after incorporating pre-trained weights from the extensively-trained GPD model (Ross et al., 2018, BSSA, https://doi.org/10.1785/0120180080) within a U-Net architecture. Please see the first accompanying notebook to this paper, [**1_UGPD_training.ipynb:** Training U-GPD Transfer Learning Model for Nabro Volcano](https://github.com/sachalapins/U-GPD/blob/main/1_UGPD_training.ipynb), for details on how to produce and train this model.  1. Install obspy(**IMPORTANT:** You must restart runtime after this step and then run notebook again from the beginning, otherwise obspy will not read or download data correctly. Click on 'Runtime' menu at top, then 'Restart runtime')<jupyter_code>!pip install obspy<jupyter_output>Requirement already satisfied: obspy in /usr/local/lib/python3.7/dist-packages (1.2.2)
Requirement already satisfied: numpy>=1.6.1 in /usr/local/lib/python3.7/dist-packages (from obspy) (1.19.5)
Requirement already satisfied: decorator in /usr/local/lib/python3.7/dist-packages (from obspy) (4.4.2)
Requirement already satisfied: scipy>=0.9.0 in /usr/local/lib/python3.7/dist-packages (from obspy) (1.4.1)
Requirement already satisfied: setuptools in /usr/local/lib/python3.7/dist-packages (from obspy) (57.2.0)
Requirement already satisfied: future>=0.12.4 in /usr/local/lib/python3.7/dist-packages (from obspy) (0.16.0)
Requirement already satisfied: sqlalchemy in /usr/local/lib/python3.7/dist-packages (from obspy) (1.4.20)
Requirement already satisfied: lxml in /usr/local/lib/python3.7/dist-packages (from obspy) (4.2.6)
Requirement already satisfied: matplotlib>=1.1.0 in /usr/local/lib/python3.7/dist-packages (from obspy) (3.2.2)
Requirement already satisfied: requests in /usr/local/lib/pyt[...]<jupyter_text>2. Load modules etcMake sure Hardware Accelerator is set to GPU (click on 'Runtime' menu at top, then 'Change runtime type', then make sure GPU is selected from drop down menu under 'Hardware accelerator').Then, let's load required modules and define some functions.<jupyter_code># Original GPD model was produced using TensorFlow (TF) version 1, so we need to tell Google Colab to use v1 (comment out this step if using conda environment above)
%tensorflow_version 1.x
# Version of h5py that comes pre-installed with Colab now doesn't like loading older keras models. Need to install earlier version before importing.
!pip install 'h5py<3.0.0'
# Import modules
import string
import time
import argparse as ap
import sys
import os
import glob

import numpy as np
import pandas as pd
import obspy.core as oc
from obspy.signal.trigger import trigger_onset
from obspy.core.utcdatetime import UTCDateTime

import matplotlib.pyplot as plt

import math # for pi, sqrt, etc

import tensorflow as tf
import keras
from keras.models import *
from keras.layers import *
from keras.optimizers import *
from keras import backend as K
from keras import losses
from keras.legacy import interfaces

import h5py
import json

import random

import scipy.stats as stats
import scipy.signal as signal

from datetime import datetime
# Adam optimiser with variable learning rates (applied through a learning rate multiplier)
# This was used during model training to set a learning rate 'multiplier' for different layers (so we can have lower learning rate for GPD layers)
# It needs to be defined so we can load our trained model

# Adapted from: https://erikbrorson.github.io/2018/04/30/Adam-with-learning-rate-multipliers/


class Adam_lr_mult(Optimizer):
    """Adam optimizer.
    Adam optimizer, with learning rate multipliers built on Keras implementation
    # Arguments
        lr: float >= 0. Learning rate.
        beta_1: float, 0 < beta < 1. Generally close to 1.
        beta_2: float, 0 < beta < 1. Generally close to 1.
        epsilon: float >= 0. Fuzz factor. If `None`, defaults to `K.epsilon()`.
        decay: float >= 0. Learning rate decay over each update.
        amsgrad: boolean. Whether to apply the AMSGrad variant of this
            algorithm from the paper "On the Convergence of Adam and
            Beyond".
    # References
        - [Adam - A Method for Stochastic Optimization](http://arxiv.org/abs/1412.6980v8)
        - [On the Convergence of Adam and Beyond](https://openreview.net/forum?id=ryQu7f-RZ)
        
    AUTHOR: <NAME>
    """

    def __init__(self, lr=0.001, beta_1=0.9, beta_2=0.999,
                 epsilon=None, decay=0., amsgrad=False,
                 multipliers=None, debug_verbose=False,**kwargs):
        super(Adam_lr_mult, self).__init__(**kwargs)
        with K.name_scope(self.__class__.__name__):
            self.iterations = K.variable(0, dtype='int64', name='iterations')
            self._lr = K.variable(lr, name='lr')
            self.beta_1 = K.variable(beta_1, name='beta_1')
            self.beta_2 = K.variable(beta_2, name='beta_2')
            self.decay = K.variable(decay, name='decay')
        if epsilon is None:
            epsilon = K.epsilon()
        self.epsilon = epsilon
        self.initial_decay = decay
        self.amsgrad = amsgrad
        self.multipliers = multipliers
        self.debug_verbose = debug_verbose

    @interfaces.legacy_get_updates_support
    def get_updates(self, loss, params):
        grads = self.get_gradients(loss, params)
        self.updates = [K.update_add(self.iterations, 1)]

        lr = self._lr
        if self.initial_decay > 0:
            lr *= (1. / (1. + self.decay * K.cast(self.iterations,
                                                  K.dtype(self.decay))))

        t = K.cast(self.iterations, K.floatx()) + 1
        lr_t = lr * (K.sqrt(1. - K.pow(self.beta_2, t)) /
                     (1. - K.pow(self.beta_1, t)))

        ms = [K.zeros(K.int_shape(p), dtype=K.dtype(p)) for p in params]
        vs = [K.zeros(K.int_shape(p), dtype=K.dtype(p)) for p in params]
        if self.amsgrad:
            vhats = [K.zeros(K.int_shape(p), dtype=K.dtype(p)) for p in params]
        else:
            vhats = [K.zeros(1) for _ in params]
        self.weights = [self.iterations] + ms + vs + vhats

        for p, g, m, v, vhat in zip(params, grads, ms, vs, vhats):

            # Learning rate multipliers
            if self.multipliers:
                multiplier = [mult for mult in self.multipliers if mult in p.name]
            else:
                multiplier = None
            if multiplier:
                new_lr_t = lr_t * self.multipliers[multiplier[0]]
                if self.debug_verbose:
                    print('Setting {} to learning rate {}'.format(multiplier[0], new_lr_t))
                    print(K.get_value(new_lr_t))
            else:
                new_lr_t = lr_t
                if self.debug_verbose:
                    print('No change in learning rate {}'.format(p.name))
                    print(K.get_value(new_lr_t))
            m_t = (self.beta_1 * m) + (1. - self.beta_1) * g
            v_t = (self.beta_2 * v) + (1. - self.beta_2) * K.square(g)
            if self.amsgrad:
                vhat_t = K.maximum(vhat, v_t)
                p_t = p - new_lr_t * m_t / (K.sqrt(vhat_t) + self.epsilon)
                self.updates.append(K.update(vhat, vhat_t))
            else:
                p_t = p - new_lr_t * m_t / (K.sqrt(v_t) + self.epsilon)

            self.updates.append(K.update(m, m_t))
            self.updates.append(K.update(v, v_t))
            new_p = p_t

            # Apply constraints.
            if getattr(p, 'constraint', None) is not None:
                new_p = p.constraint(new_p)

            self.updates.append(K.update(p, new_p))
        return self.updates

    def get_config(self):
        config = {'lr': float(K.get_value(self._lr)),
                  'beta_1': float(K.get_value(self.beta_1)),
                  'beta_2': float(K.get_value(self.beta_2)),
                  'decay': float(K.get_value(self.decay)),
                  'epsilon': self.epsilon,
                  'amsgrad': self.amsgrad,
                  'multipliers':self.multipliers}
        base_config = super(Adam_lr_mult, self).get_config()
        return dict(list(base_config.items()) + list(config.items()))
# Define function to produce sliding windows of signal (from original GPD GitHub repo https://github.com/interseismic/generalized-phase-detection/blob/master/gpd_predict.py)

def sliding_window(data, size, stepsize=1, padded=False, axis=-1, copy=True):
    """
    Calculate a sliding window over a signal
    Parameters
    ----------
    data : numpy array
        The array to be slided over.
    size : int
        The sliding window size
    stepsize : int
        The sliding window stepsize. Defaults to 1.
    axis : int
        The axis to slide over. Defaults to the last axis.
    copy : bool
        Return strided array as copy to avoid sideffects when manipulating the
        output array.
    Returns
    -------
    data : numpy array
        A matrix where row in last dimension consists of one instance
        of the sliding window.
    Notes
    -----
    - Be wary of setting `copy` to `False` as undesired sideffects with the
      output values may occurr.
    Examples
    --------
    >>> a = numpy.array([1, 2, 3, 4, 5])
    >>> sliding_window(a, size=3)
    array([[1, 2, 3],
           [2, 3, 4],
           [3, 4, 5]])
    >>> sliding_window(a, size=3, stepsize=2)
    array([[1, 2, 3],
           [3, 4, 5]])
    See Also
    --------
    pieces : Calculate number of pieces available by sliding
    """
    if axis >= data.ndim:
        raise ValueError(
            "Axis value out of range"
        )

    if stepsize < 1:
        raise ValueError(
            "Stepsize may not be zero or negative"
        )

    if size > data.shape[axis]:
        raise ValueError(
            "Sliding window size may not exceed size of selected axis"
        )

    shape = list(data.shape)
    shape[axis] = np.floor(data.shape[axis] / stepsize - size / stepsize + 1).astype(int)
    shape.append(size)

    strides = list(data.strides)
    strides[axis] *= stepsize
    strides.append(data.strides[axis])

    strided = np.lib.stride_tricks.as_strided(
        data, shape=shape, strides=strides
    )

    if copy:
        return strided.copy()
    else:
        return strided<jupyter_output><empty_output><jupyter_text>3. Download and load trained U-GPD Transfer Learning modelNext let's clone the repo for our trained U-GPD Transfer Learning model and load our trained model with weights from the training epoch that gave lowest validation loss.<jupyter_code>!git clone https://github.com/sachalapins/U-GPD.git
# Load model and weights
model = keras.models.load_model('./U-GPD/model/best/model.h5',
              custom_objects={'Adam_lr_mult': Adam_lr_mult})
print("*** 1 of 2: MODEL LOADED FROM DISK ***")

model.load_weights('./U-GPD/model/best/weights_lowest_loss.hdf5') # Load in weights from training epoch with lowest validation loss

print ("*** 2 of 2:BEST TRAINING WEIGHTS LOADED ***")<jupyter_output>WARNING:tensorflow:From /tensorflow-1.15.2/python3.7/tensorflow_core/python/ops/resource_variable_ops.py:1630: calling BaseResourceVariable.__init__ (from tensorflow.python.ops.resource_variable_ops) with constraint is deprecated and will be removed in a future version.
Instructions for updating:
If using Keras pass *_constraint arguments to layers.
WARNING:tensorflow:From /tensorflow-1.15.2/python3.7/keras/backend/tensorflow_backend.py:4070: The name tf.nn.max_pool is deprecated. Please use tf.nn.max_pool2d instead.

WARNING:tensorflow:From /tensorflow-1.15.2/python3.7/tensorflow_core/python/ops/math_grad.py:1424: where (from tensorflow.python.ops.array_ops) is deprecated and will be removed in a future version.
Instructions for updating:
Use tf.where in 2.0, which has the same broadcast rule as np.where
WARNING:tensorflow:From /tensorflow-1.15.2/python3.7/keras/backend/tensorflow_backend.py:422: The name tf.global_variables is deprecated. Please use tf.compat.v1.global_variables inst[...]<jupyter_text>4. Download 24 hours data from NabroGitHub doesn't allow uploading of large files so we need to download our Nabro data from the [IRIS dataselect Web Service](http://service.iris.edu/fdsnws/dataselect/1/) before we can run our model.Let's choose a day that wasn't included in any of the training, validation or test data (2011-09-26).**Note: This step can take a long time (over half an hour) so I have downloaded this data previously and saved it to my Google Drive. To run the model on this data yourself, simply uncomment the code block below (but be prepared to wait). You can also download it through the obspy client but this appears to take even longer.**<jupyter_code># # Give start and end time to extract data
# station = "NAB*"
# day = "2011-09-26"
# starttime = UTCDateTime(day + "T00:00:00")
# endtime = UTCDateTime(day + "T23:59:59.999999")

# # Get waveforms for network 4H on this day and save as mseed (seems to be a little quicker than using the obspy client to request data)
# before = datetime.now()
# os.system("wget -O nabro_2011-09-26.mseed http://service.iris.edu/fdsnws/dataselect/1/query?net=4H\&sta=" + station + "\&loc=*\&cha=HH?\&start=" + str(starttime)[0:-1] + "\&end=" + str(endtime)[0:-1])
# after = datetime.now()
# print(after-before)
# Load waveforms into obspy
st = oc.Stream()
# st += oc.read("nabro_2011-09-26.mseed") # Uncomment this line if downloaded using code cell above
st += oc.read("./drive/MyDrive/nabro_2011-09-26.mseed") # Using data downloaded previously to my Google Drive - comment out this line if data downloaded using code cell above
st.sort() # Sort obspy stream data by station and channel names<jupyter_output><empty_output><jupyter_text>The obspy Stream has 513 traces, even though there are only 7 operational stations on this day. This is because station NAB7 drops in and out every few minutes.As such, to avoid lots of spurious triggers caused by these data gaps, let's remove NAB7 from this Stream.<jupyter_code>for s in st:
    if s.stats.station == "NAB7":
        st.remove(s)

print(st)<jupyter_output>18 Trace(s) in Stream:
4H.NAB1..HHE | 2011-09-26T00:00:00.000000Z - 2011-09-26T23:59:59.990000Z | 100.0 Hz, 8640000 samples
4H.NAB1..HHN | 2011-09-26T00:00:00.000000Z - 2011-09-26T23:59:59.990000Z | 100.0 Hz, 8640000 samples
4H.NAB1..HHZ | 2011-09-26T00:00:00.000000Z - 2011-09-26T23:59:59.990000Z | 100.0 Hz, 8640000 samples
4H.NAB2..HHE | 2011-09-26T00:00:00.000000Z - 2011-09-26T23:59:59.990000Z | 100.0 Hz, 8640000 samples
4H.NAB2..HHN | 2011-09-26T00:00:00.000000Z - 2011-09-26T23:59:59.990000Z | 100.0 Hz, 8640000 samples
4H.NAB2..HHZ | 2011-09-26T00:00:00.000000Z - 2011-09-26T23:59:59.990000Z | 100.0 Hz, 8640000 samples
4H.NAB3..HHE | 2011-09-26T00:00:00.000000Z - 2011-09-26T23:59:59.990000Z | 100.0 Hz, 8640000 samples
4H.NAB3..HHN | 2011-09-26T00:00:00.000000Z - 2011-09-26T23:59:59.990000Z | 100.0 Hz, 8640000 samples
4H.NAB3..HHZ | 2011-09-26T00:00:00.000000Z - 2011-09-26T23:59:59.990000Z | 100.0 Hz, 8640000 samples
4H.NAB4..HHE | 2011-09-26T00:00:00.000000Z - 2011-09-26T23:59:59.990[...]<jupyter_text>5. Define some parameters for sliding window and phase arrival trigger<jupyter_code># Model prediction parameters
shift_size = 200 # Overlap (in samples) for sliding window
batch_size = 100 # No. of windows to process at a time

# Trigger parameters
p_thresh = 0.4 # Trigger on threshold for P-waves
s_thresh = 0.4 # Trigger on threshold for S-waves
trig_off = 0.2 # Trigerr off threshold for both phase arrivals<jupyter_output><empty_output><jupyter_text>6. Run model using sliding window on continuous dataIn this step, we loop through each available station, dividing the data from each instrument channel into overlapping 400-sample windows (with 200-sample overlap) and then run our model on these windows.We keep only the middle 200 samples from each window to avoid poor predictions at the window edges.We then apply a simple trigger algorithm in obspy to a concatenated array of all these window predictions (one long prediction trace with same sample rate as the input data). Any P- or S-wave predictions that exceed the trigger thresholds given above (`p_thresh` and `s_thresh`) will be added to the dataframe `df` of phase arrival pick times, along with their 'probabilities' (a measure of model confidence in that phase arrival classification).<jupyter_code># List available stations
stations = list(set([st[s].stats.station for s in np.arange(0, len(st))]))

# Dataframe to store picks from trigger
df = pd.DataFrame(columns=['time', 'sta', 'pha', 'prob'])

# Loop through list of available stations
before = datetime.now() # Get time before and after to see how long this takes
for s in np.arange(0, len(stations)):

    st2 = st.select(station=stations[s])
    
    # If there are data gaps (i.e., more than 3 traces per station) then merge the traces together
    no_of_traces_per_station = st2.count() / 3
    
    if no_of_traces_per_station > 1:
        st2.merge(method=0, fill_value=0)
    
    # Trim to same length:
    latest_start = np.max([x.stats.starttime for x in st2])
    earliest_stop = np.min([x.stats.endtime for x in st2])
    st2.trim(latest_start, earliest_stop)        
    
    # Expand overlapping windows using sliding_window function
    sliding_E = sliding_window(st2[0].data, 400, stepsize=shift_size) # Slide window across first component
    tr_win = np.zeros((sliding_E.shape[0], 400, 3)) # 3D array (of zeros)
    tr_win[:,:,0] = sliding_E # Add to 3D array
    sliding_E = None # clear some memory
    
    sliding_N = sliding_window(st2[1].data, 400, stepsize=shift_size) # Slide window across second component
    tr_win[:,:,1] = sliding_N # Add to 3D array
    sliding_N = None # clear some memory
    
    sliding_Z = sliding_window(st2[2].data, 400, stepsize=shift_size) # Slide window across third component
    tr_win[:,:,2] = sliding_Z # Add to 3D array
    sliding_Z = None # clear some memory
    
    tr_win = signal.detrend(tr_win, axis=1) # Detrend 
    tr_win = tr_win / np.max(np.abs(tr_win), axis=(1,2))[:,None,None] # Normalize between 0 and 1 (divide by max across 3 components for each window)
    
    # Model prediction step...
    ts = model.predict(tr_win, verbose=True, batch_size=batch_size)
    
    # Merge overlapping predictions:
    n_steps_per_win = 400/shift_size # No. of steps per window (i.e. how many overlapping windows)
    class_trace = np.zeros((len(st2[0].data), int(3))) # Create prediction traces to be same size as input data
    for i in np.arange(0, ts.shape[0]): # For each window
        class_trace[((i * shift_size) + 100):((i * shift_size) + 300), :] += ts[i, 100:300, :] # Keep middle 200 predictions

    # Use obspy trigger to find predictions above threshold (i.e. p_thresh and s_thresh)
    # P-wave picks (trigger):
    trigs = trigger_onset(class_trace[:,0], p_thresh, trig_off) # Look for triggers (prediction trace exceeding p_thresh)
    p_picks = []
    p_probs = []
    for trig in trigs:
        if trig[1] == trig[0]:
            continue # If trigger on and off times are the same then ignore
        pick = np.argmax(class_trace[trig[0]:trig[1], 0]) + trig[0] # Use argmax between trigger on and off times as pick time
        stamp_pick = st2[0].stats.starttime + (pick * st2[0].stats.delta) # Pick time as UTCDateTime
        p_picks.append(stamp_pick) # Append time to list of P pick times
        p_probs.append(class_trace[pick, 0]) # Append prob to list of P pick probabilities
    # Add picks to df for this day            
    if len(p_picks) > 0:
        df_add = pd.DataFrame(p_picks, columns=['time'])
        df_add['sta'] = stations[s]
        df_add['pha'] = 'p'
        df_add['prob'] = p_probs
        df = pd.concat([df, df_add])
        
    # Repeat above for S-wave picks and s_thresh:
    trigs = trigger_onset(class_trace[:,1], s_thresh, trig_off)
    s_picks = []
    s_probs = []
    for trig in trigs:
        if trig[1] == trig[0]:
            continue
        pick = np.argmax(class_trace[trig[0]:trig[1], 1]) + trig[0]
        stamp_pick = st2[0].stats.starttime + (pick * st2[0].stats.delta)
        s_picks.append(stamp_pick)
        s_probs.append(class_trace[pick, 1])
    # Add picks to df for this day
    if len(s_picks) > 0:
        df_add = pd.DataFrame(s_picks, columns=['time'])
        df_add['sta'] = stations[s]
        df_add['pha'] = 's'
        df_add['prob'] = s_probs
        df = pd.concat([df, df_add])

df = df.sort_values(by='time').reset_index(drop=True) # Sort df by pick times
#df.to_pickle("trigs.pkl") # Write df to pkl file - uncomment if you want to save your model picks

after = datetime.now()
time_taken = after-before
print("24 hours data from " + str(len(stations)) + " stations processed in " + str(time_taken.seconds) + "." + str(time_taken.microseconds) + " seconds.")<jupyter_output>43199/43199 [==============================] - 13s 295us/step
43199/43199 [==============================] - 5s 106us/step
43199/43199 [==============================] - 5s 107us/step
43199/43199 [==============================] - 5s 108us/step
43199/43199 [==============================] - 5s 108us/step
43199/43199 [==============================] - 5s 108us/step
24 hours data from 6 stations processed in 57.39385 seconds.<jupyter_text>7. Associate predicted phase arrivals for event locationIn our paper, we locate events using model picks using the following approach:- First group P-wave arrivals into 4-second 'bins'. This should encompass the largest arrival time difference between all our stations (i.e., 4 seconds is larger than the maximum theoretical travel time between any two stations).- We then locate events using just these P-wave groupings and use the origin times to determine which S-wave picks should be associated to each P-wave grouping. A Vp/Vs ratio of 1.76 and travel time error of 25% is used to do this.- Finally, events are located again using associated P- and S-wave picks.Before associating picks, let's look at our dataframe of model pick times:<jupyter_code># All picks
df
# By station
for s in list(set(df['sta'])):
    no_p_picks = df.loc[(df['sta'] == s) & (df['pha'] == "p")].shape[0]
    no_s_picks = df.loc[(df['sta'] == s) & (df['pha'] == "s")].shape[0]
    print("Station " + s + ": " + str(no_p_picks) + " P-waves and " + str(no_s_picks) + " S-waves.")<jupyter_output>Station NAB5: 212 P-waves and 170 S-waves.
Station NAB2: 584 P-waves and 315 S-waves.
Station NAB1: 617 P-waves and 269 S-waves.
Station NAB8: 370 P-waves and 292 S-waves.
Station NAB4: 498 P-waves and 824 S-waves.
Station NAB3: 289 P-waves and 165 S-waves.<jupyter_text>We have a few thousand picks and potentially hundreds of events in this day (this is an active volcano after all!).We can see at the start and end of this dataframe (`df`) that we have P- and S-wave picks across several stations that occur very close together in time. This is a good sign that we are detecting true phase arrivals and that we will be able to locate their events.Let's group P-wave arrivals into 4 second bins and see how many locatable events we have.<jupyter_code>df_p = df.loc[df['pha'] == "p"] # Extract just P-wave picks
df_p.reset_index(drop=True, inplace=True) # Reset dataframe indexing

pick_groups = [] # Create empty list to store P-wave groupings

# First loop through each pick finding all picks within 4 seconds of that pick - this takes around 20 secs
for i in range(df_p.shape[0]):
    this_time = df_p['time'][i]
    this_bin = df_p.loc[(df_p['time'] >= this_time) & (df_p['time'] <= (this_time + 4))]
    
    # If 4+ stations then add to pick_groups list
    if this_bin.shape[0] >= 4:
        if len(list(set(this_bin['sta']))) >= 4:
            pick_groups.append(this_bin)

# Remove groups that are fully contained in another
# First 'backwards' (remove groups that are fully contained in the previous group)
for i in reversed(range(len(pick_groups))):
    if i == 0:
        continue
    else:
        if pick_groups[i].isin(pick_groups[(i-1)]).all().all():
            del(pick_groups[i])
# Then 'forwards' (remove groups that are fully contained in the next group)
for i in reversed(range(len(pick_groups)-1)):
    if pick_groups[i].isin(pick_groups[(i+1)]).all().all():
        del(pick_groups[i])

print("Total no. of P-wave groups = " + str(len(pick_groups)))<jupyter_output>Total no. of P-wave groups = 171<jupyter_text>Let's look at our groupings so far...<jupyter_code>pick_groups<jupyter_output><empty_output><jupyter_text>These groupings look pretty good already, although some have more than one P-wave pick per station and some overlap with other groups, usually where one group incorporates an earlier pick and the next group incorporates a later pick (outside of the 4 second window). We'll see below that this is a relatively small proportion of groups.To deal with having more than one pick per station, let's increase the probability threshold in those groups (so we can have higher confidence in these picks) and then keep the picks with highest probability for each station.After that, let's remove any overlapping groups by keeping the group with highest mean pick probability. Another option could be to merge overlapping groups, but that is a matter of choice. We are broadly trying to evaluate our model's suitability for detecting events in continuous data rather than assess our phase association method.<jupyter_code># Check for any groups that have more than 1 pick - first increase threshold for all picks then use phase w/ highest prob
for i in reversed(range(len(pick_groups))):
    if pick_groups[i]['sta'].duplicated().any(): # If more than one pick for any station
        pick_groups[i] = pick_groups[i].loc[pick_groups[i]['prob'] >= 0.6] # First only keep picks above higher threshold (0.6)

        if len(list(set(pick_groups[i]['sta']))) < 4:
            del(pick_groups[i]) # Remove group if increased threshold reduces number of stations to less than 4
        else:
            # Otherwise check if any stations still have more than one pick and keep only highest probability pick
            if pick_groups[i]['sta'].duplicated().any():
                pick_groups[i] = pick_groups[i].sort_values('prob').drop_duplicates('sta', keep='last')

print("No. of P-wave groups after increasing threshold for groups with multiple picks per station and reducing no. of picks to one per station = " + str(len(pick_groups)))

# Then check for overlapping groups and choose group with highest mean probability
for i in reversed(range(len(pick_groups)-1)):
    if pick_groups[i].index.isin(pick_groups[i+1].index).any():
        if np.mean(pick_groups[i]['prob']) >= np.mean(pick_groups[i+1]['prob']):
            del(pick_groups[i+1])
        else:
            del(pick_groups[i])

print("No. of P-wave groups after removing overlapping groups = " + str(len(pick_groups)))<jupyter_output>No. of P-wave groups after increasing threshold for groups with multiple picks per station and reducing no. of picks to one per station = 156
No. of P-wave groups after removing overlapping groups = 146<jupyter_text>The above grouping algorithm leaves us with 146 locatable events.The next step is to put them in a format for locating events in [NonLinLoc](http://alomax.free.fr/nlloc/), a popular software package for probabilistic earthquake location.<jupyter_code># Make directories to store NLL pickfile, velocity model, travel time grid and locations
!mkdir ./NLL
!mkdir ./NLL/picks
!mkdir ./NLL/model
!mkdir ./NLL/time
!mkdir ./NLL/loc
# Loop through groups of picks, producing strings to go in NonLinLoc phase input file (e.g., http://alomax.free.fr/nlloc/soft7.00/formats.html#_phase_nlloc_)
# And then write strings to NLL phase input file

phase_file_str = [] # List to store strings for NLL phase input file

# Loop through groups and picks producing appropriate strings for phase input file
for i in range(len(pick_groups)):
    pick_groups[i].reset_index(drop=True, inplace=True) # Reset indices for each dataframe for looping through
    for j in range(pick_groups[i].shape[0]):
        if (i > 0) and (j == 0):
            phase_file_str.append(" ") # Add a blank line after each event
        this_line = pick_groups[i]['sta'][j].ljust(7) # Station
        this_line += '?'.ljust(5) # Instrument (leave as ?)
        this_line += "?".ljust(5) # Component (leave as ?)
        this_line += "? P".ljust(9) # Onset (leave as ? as don't know if impulsive or emergent) and phase descriptor (P)
        this_line += "? " # First motion (i.e. compression/dilation, leave as ? as don't know)
        this_line += str(pick_groups[i]['time'][j].year) # Year
        this_line += str(pick_groups[i]['time'][j].month).zfill(2) # Month (2 digits)
        this_line += str(pick_groups[i]['time'][j].day).zfill(2) # Day (2 digits)
        this_line += " " + str(pick_groups[i]['time'][j].hour).zfill(2) # Hour (2 digits)
        this_line += str(pick_groups[i]['time'][j].minute).zfill(2) # Minute (2 digits)
        this_line += " " + str(pick_groups[i]['time'][j].second).zfill(2) # Second (2 digits)
        this_line += "." + str(pick_groups[i]['time'][j].microsecond).zfill(4)[:4] # Microseconds (4 digits)
        this_line += " GAU " # Error type = Gaussian
        if pick_groups[i]['prob'][j] >= 0.85:
            this_line += "5.00e-02 " # 0.05 sec err for prob >= 0.85
        elif pick_groups[i]['prob'][j] >= 0.7:
            this_line += "1.00e-01 " # 0.10 sec err for prob >= 0.7
        elif pick_groups[i]['prob'][j] >= 0.55:
            this_line += "2.00e-01 " # 0.20 sec err for prob >= 0.55
        else:
            this_line += "3.33e-01 " # 0.333 sec err for prob >= 0.4
        this_line += "-1.00e+00 " # Coda duration (negative = don't use)
        this_line += "-1.00e+00 " # Max peak-to-peak amplitude (negative = don't use)
        this_line += "-1.00e+00" # Period of amplitude reading (negative = don't use)

        phase_file_str.append(this_line) # Append line to list

# Write phase input file
with open('./NLL/picks/pickfile.hpf', 'w') as f:
    for item in phase_file_str:
        f.write("%s\n" % item)<jupyter_output><empty_output><jupyter_text>Now we need to install NonLinLoc (seems to work but lots of warnings during `make distrib`).<jupyter_code>!wget http://alomax.free.fr/nlloc/soft7.00/tar/NLL7.00_src.tgz; tar -zxf NLL7.00_src.tgz; cd src; make distrib<jupyter_output>--2021-07-28 16:20:11--  http://alomax.free.fr/nlloc/soft7.00/tar/NLL7.00_src.tgz
Resolving alomax.free.fr (alomax.free.fr)... 192.168.3.11
Connecting to alomax.free.fr (alomax.free.fr)|192.168.3.11|:80... connected.
HTTP request sent, awaiting response... 200 OK
Length: 407847 (398K) [application/x-gzip]
Saving to: ‘NLL7.00_src.tgz’

NLL7.00_src.tgz     100%[===================>] 398.29K  1.01MB/s    in 0.4s    

2021-07-28 16:20:12 (1.01 MB/s) - ‘NLL7.00_src.tgz’ saved [407847/407847]

gcc -c -O3 -Wall -std=gnu99  NLLoc_main.c  
gcc -c -O3 -Wall -std=gnu99  NLLoc1.c  
[01m[KNLLoc1.c:[m[K In function ‘[01m[KNLLoc[m[K’:
[01m[KNLLoc1.c:277:9:[m[K [01;35m[Kwarning: [m[Kignoring return value of ‘[01m[Ksystem[m[K’, declared with attribute warn_unused_result [[01;35m[K-Wunused-result[m[K]
         [01;35m[Ksystem(sys_command)[m[K;
         [01;35m[K^~~~~~~~~~~~~~~~~~~[m[K
[01m[KNLLoc1.c:280:9:[m[K [01;35m[Kwarning: [m[Kignoring return value of ‘[01[...]<jupyter_text>Then we need to create velocity model and travel time grids for each station for locating events (using pre-made Control File in U-GPD GitHub repo).<jupyter_code>!./src/Vel2Grid ./U-GPD/NLL/NLL_master.in
!./src/Grid2Time ./U-GPD/NLL/NLL_master.in<jupyter_output>Recursive initialization: level 1
Rediscretized timefield dimensions: 31 31 16
Homogeneous region: x[14->15] y[15->16] z[0->2]

Starting F.D. computation...
Starting F.D. computation...
x_side 81->80: ff fb bb bf 
y_side 79->78: ff fb bb bf 
x_side 95->96: ff fb bb bf 
y_side 93->94: ff fb bb bf 
z_side 7->8: ff fb bb bf 
x_side 80->79: ff fb bb bf 
y_side 78->77: ff fb bb bf 
x_side 96->97: ff fb bb bf 
y_side 94->95: ff fb bb bf 
z_side 8->9: ff fb bb bf 
x_side 79->78: ff fb bb bf 
y_side 77->76: ff fb bb bf 
x_side 97->98: ff fb bb bf 
y_side 95->96: ff fb bb bf 
z_side 9->10: ff fb bb bf 
x_side 78->77: ff fb bb bf 
y_side 76->75: ff fb bb bf 
x_side 98->99: ff fb bb bf 
y_side 96->97: ff fb bb bf 
z_side 10->11: ff fb bb bf 
x_side 77->76: ff fb bb bf 
y_side 75->74: ff fb bb bf 
x_side 99->100: ff fb bb bf 
y_side 97->98: ff fb bb bf 
z_side 11->12: ff fb bb bf 
x_side 76->75: ff fb bb bf 
y_side 74->73: ff fb bb bf 
x_side 100->101: ff fb bb bf 
y_side 98->99: ff fb bb bf 
z_s[...]<jupyter_text>And then finally locate events in NonLinLoc<jupyter_code>!./src/NLLoc ./U-GPD/NLL/NLL_master.in<jupyter_output>Finished event location, output files: ./NLL/loc/ALL.20110926.000137.* <./NLL/loc/ALL.20110926.000137.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc/ALL.20110926.002946.* <./NLL/loc/ALL.20110926.002946.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc/ALL.20110926.003202.* <./NLL/loc/ALL.20110926.003202.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc/ALL.20110926.003231.* <./NLL/loc/ALL.20110926.003231.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc/ALL.20110926.004338.* <./NLL/loc/ALL.20110926.004338.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc/ALL.20110926.004549.* <./NLL/loc/ALL.20110926.004549.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc/ALL.20110926.005500.* <./NLL/loc/ALL.20110926.005500.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc/ALL.20110926.005737.* <./NLL/loc/ALL.20110926.005737.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc/ALL.20110926.01[...]<jupyter_text>Now we have event locations from our P-wave groupings, we can use these to determine which S-wave picks should be included for each event.The code cell below extracts event origin time info from NonLinLoc output files, calculates the difference between our P-wave arrival times and event origin times, and then looks through our model's S-wave picks to see if any should be associated to our events (using a Vp/Vs ratio of 1.76 and 25% travel time error).<jupyter_code># Read in origin times from NLL location summary file:
with open("./NLL/loc/ALL.sum.grid0.loc.hyp","r") as fi:
    origins = []
    for ln in fi:
        if ln.startswith("GEOGRAPHIC"):
            origins.append(ln[15:42]) # Take origin time info from GEOGRAPHIC line in summary hyp file

# Loop through and convert origin times to UTCDateTime
for i in range(len(origins)):
    origins[i] = UTCDateTime(origins[i][:4] + origins[i][5:7] + origins[i][8:10] + "T" + origins[i][12:14] + origins[i][15:17] + ("%.4f" % float(origins[i][18:27])).zfill(7))

# Get dataframe of model S-wave picks
df_s = df.loc[df['pha'] == "s"] # Extract just S-wave picks
df_s.reset_index(drop=True, inplace=True) # Reset dataframe indexing

# Then loop through P-wave pick groups finding any S-waves to add:
for i in range(len(pick_groups)):
    # Add columns for new S-wave times and probabilities:
    pick_groups[i]['s_time'] = pick_groups[i]['time'] # Just put S-wave time = P-wave time initially
    pick_groups[i]['s_prob'] = 0. # Put zero probability initially
    # Loop through picks in each group
    for j in range(pick_groups[i].shape[0]):
        p_origin_diff = pick_groups[i]['time'][j] - origins[i] # P arrival time minus origin time
        s_window = df_s.loc[(df_s['time'] >= origins[i] + (p_origin_diff * 1.76 * 0.75)) & (df_s['time'] <= origins[i] + (p_origin_diff * 1.76 * 1.25))] # Any picks within predicted S arrival window?
        s_window = s_window.loc[s_window['sta'] == pick_groups[i]['sta'][j]] # Any S picks within window and from this station?
        # If so...
        if s_window.shape[0] > 0:
            if s_window.shape[0] > 1:
                # If more than one pick for this station, use pick with highest probability
                s_window = s_window.sort_values('prob').drop_duplicates('sta', keep='last')
            s_window.reset_index(inplace=True, drop=True) # Reset index
            # Add pick to pick_groups dataframe
            pick_groups[i]['s_time'][j] = s_window['time'][0]
            pick_groups[i]['s_prob'][j] = s_window['prob'][0]
# A quick look at our new groupings
pick_groups<jupyter_output><empty_output><jupyter_text>Then finally re-write our NonLinLoc phase input file and run it again to locate events with P's and S's.<jupyter_code># Rewrite NLL phase input file:

phase_file_str = [] # Empty list to store strings for NLL phase input file

# Loop through groups and picks producing appropriate strings for phase input file
for i in range(len(pick_groups)):
    pick_groups[i].reset_index(drop=True, inplace=True) # Reset indices for each dataframe for looping through
    for j in range(pick_groups[i].shape[0]):
        if (i > 0) and (j == 0):
            phase_file_str.append(" ") # Add a blank line after each event
        this_line = pick_groups[i]['sta'][j].ljust(7) # Station
        this_line += '?'.ljust(5) # Instrument (leave as ?)
        this_line += "?".ljust(5) # Component (leave as ?)
        this_line += "? P".ljust(9) # Onset (leave as ? as don't know if impulsive or emergent) and phase descriptor (P)
        this_line += "? " # First motion (i.e. compression/dilation, leave as ? as don't know)
        this_line += str(pick_groups[i]['time'][j].year) # Year
        this_line += str(pick_groups[i]['time'][j].month).zfill(2) # Month (2 digits)
        this_line += str(pick_groups[i]['time'][j].day).zfill(2) # Day (2 digits)
        this_line += " " + str(pick_groups[i]['time'][j].hour).zfill(2) # Hour (2 digits)
        this_line += str(pick_groups[i]['time'][j].minute).zfill(2) # Minute (2 digits)
        this_line += " " + str(pick_groups[i]['time'][j].second).zfill(2) # Second (2 digits)
        this_line += "." + str(pick_groups[i]['time'][j].microsecond).zfill(4)[:4] # Microseconds (4 digits)
        this_line += " GAU " # Error type = Gaussian
        if pick_groups[i]['prob'][j] >= 0.85:
            this_line += "5.00e-02 " # 0.05 sec err for prob >= 0.85
        elif pick_groups[i]['prob'][j] >= 0.7:
            this_line += "1.00e-01 " # 0.10 sec err for prob >= 0.7
        elif pick_groups[i]['prob'][j] >= 0.55:
            this_line += "2.00e-01 " # 0.20 sec err for prob >= 0.55
        else:
            this_line += "3.33e-01 " # 0.333 sec err for prob >= 0.4
        this_line += "-1.00e+00 " # Coda duration (negative = don't use)
        this_line += "-1.00e+00 " # Max peak-to-peak amplitude (negative = don't use)
        this_line += "-1.00e+00" # Period of amplitude reading (negative = don't use)

        phase_file_str.append(this_line) # Append line to list

        # Check if there is an associated S-wave pick for this station and append to the phase file:
        if pick_groups[i]['s_prob'][j] >= 0.4:
            this_line = pick_groups[i]['sta'][j].ljust(7) # Station
            this_line += '?'.ljust(5) # Instrument (leave as ?)
            this_line += "?".ljust(5) # Component (leave as ?)
            this_line += "? S".ljust(9) # Onset (leave as ? as now S-wave) and phase descriptor (S)
            this_line += "? " # First motion (i.e. compression/dilation, leave as ? as now S-wave)
            this_line += str(pick_groups[i]['s_time'][j].year) # Year
            this_line += str(pick_groups[i]['s_time'][j].month).zfill(2) # Month (2 digits)
            this_line += str(pick_groups[i]['s_time'][j].day).zfill(2) # Day (2 digits)
            this_line += " " + str(pick_groups[i]['s_time'][j].hour).zfill(2) # Hour (2 digits)
            this_line += str(pick_groups[i]['s_time'][j].minute).zfill(2) # Minute (2 digits)
            this_line += " " + str(pick_groups[i]['s_time'][j].second).zfill(2) # Second (2 digits)
            this_line += "." + str(pick_groups[i]['s_time'][j].microsecond).zfill(4)[:4] # Microseconds (4 digits)
            this_line += " GAU " # Error type = Gaussian
            # Keep same error magnitudes as P-wave arrivals
            if pick_groups[i]['prob'][j] >= 0.85:
                this_line += "5.00e-02 " # 0.05 sec err for prob >= 0.85
            elif pick_groups[i]['prob'][j] >= 0.7:
                this_line += "1.00e-01 " # 0.10 sec err for prob >= 0.7
            elif pick_groups[i]['prob'][j] >= 0.55:
                this_line += "2.00e-01 " # 0.20 sec err for prob >= 0.55
            else:
                this_line += "3.33e-01 " # 0.333 sec err for prob >= 0.4
            this_line += "-1.00e+00 " # Coda duration (negative = don't use)
            this_line += "-1.00e+00 " # Max peak-to-peak amplitude (negative = don't use)
            this_line += "-1.00e+00" # Period of amplitude reading (negative = don't use)

            phase_file_str.append(this_line) # Append line to list

# Write new phase file
with open('./NLL/picks/pickfile_s.hpf', 'w') as f:
    for item in phase_file_str:
        f.write("%s\n" % item)
# Make new directory for locating events with S-waves included:
!mkdir ./NLL/loc_s

# Run NonLinLoc again (using pre-made Control File for including S-waves)
!./src/NLLoc ./U-GPD/NLL/NLL_master_s.in<jupyter_output>Finished event location, output files: ./NLL/loc_s/ALL.20110926.000137.* <./NLL/loc_s/ALL.20110926.000137.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc_s/ALL.20110926.002946.* <./NLL/loc_s/ALL.20110926.002946.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc_s/ALL.20110926.003202.* <./NLL/loc_s/ALL.20110926.003202.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc_s/ALL.20110926.003231.* <./NLL/loc_s/ALL.20110926.003231.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc_s/ALL.20110926.004338.* <./NLL/loc_s/ALL.20110926.004338.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc_s/ALL.20110926.004549.* <./NLL/loc_s/ALL.20110926.004549.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc_s/ALL.20110926.005500.* <./NLL/loc_s/ALL.20110926.005500.grid0.loc.hyp>
Finished event location, output files: ./NLL/loc_s/ALL.20110926.005737.* <./NLL/loc_s/ALL.20110926.005737.grid0.loc.hyp>
Finished event location, output [...]<jupyter_text>8. Final event locationsThe very last step is to read in our final event locations for plotting or further analysis. You can obviously extract other information from the NonLinLoc output files in `./NLL/loc_s`, such as horizontal and vertical Gaussian errors (as presented in our paper) or importance sampling outputs (e.g., pdf scatter and likelihood values for each event).Phase information from our model picks can be accessed from the NonLinLoc phase input file we produced earlier (`./NLL/picks/pickfile_s.hpf`) and used to isolate phase arrivals for calulating moment magnitudes, etc.The plots in our original paper were all produced in R so for now we'll just read the event locations into a pandas dataframe and leave you to experiment further.<jupyter_code>origins = []
lats = []
lons = []
depths = []

with open("./NLL/loc_s/ALL.sum.grid0.loc.hyp","r") as fi:
    for ln in fi:
        if ln.startswith("GEOGRAPHIC"):
            origins.append(ln[15:42])
            lats.append(float(ln[47:57]))
            lons.append(float(ln[62:72]))
            depths.append(float(ln[78:]))

# # Loop through and convert origin times to UTCDateTime
for i in range(len(origins)):
    origins[i] = UTCDateTime(origins[i][:4] + origins[i][5:7] + origins[i][8:10] + "T" + origins[i][12:14] + origins[i][15:17] + ("%.4f" % float(origins[i][18:27])).zfill(7))

# Create dataframe of event locations
df_events = pd.DataFrame(list(zip(origins, lats, lons, depths)), 
               columns =['time', 'lat', 'lon', 'depth'])
# Quick look at events dataframe
df_events<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12731795065402296, 'fraction_numerical': 0.05105598614226387, 'mean_word_length': 3.3331408775981526, 'pattern_counts': {'":': 2, '<': 91, '<?xml version=': 0, '>': 169, 'https://': 9, 'lorem ipsum': 0, 'www.': 0, 'xml': 1}, 'pii_count': 2, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '73218', 'license': 'MIT', 'n_tokens_mistral': 16541, 'n_tokens_neox': 14831, 'n_words': 4608, 'path': '2_UGPD_run.ipynb', 'repo_name': 'sachalapins/U-GPD', 'size': '45081'}