This code is necessary on colab to install SeisBench. If SeisBench is already installed on your machine, you can skip this.

!pip install seisbench

This cell is required to circumvent an issue with colab and obspy. For details, check this issue in the obspy documentation: obspy/obspy#2547

try:
    import obspy
    obspy.read()
except TypeError:
    # Needs to restart the runtime once, because obspy only works properly after restart.
    print('Stopping RUNTIME. If you run this code for the first time, this is expected. Colaboratory will restart automatically. Please run again.')
    exit()

Quickly setting-up picking models

This tutorial is inspired by an INGV workshop, where a range of varying seismic examples were preseted and the models available through SeisBench were deployed to picks / detect events on the seismic streams. This notebook does not perform in-depth analysis of the resulting detection information, it just shows how you can quickly access and deploy a range of seismic models to seismic data.

Note: Some familiarity with obspy is helpful for this tutorial, but not required.

import seisbench
import seisbench.models as sbm

Creating a model

Here you can download your pretrained model for predicting on the example streams.

if you would like to see the full suite of available pre-trained weights for any model, you can call the list_pretrained method, which will return all the associated weights for the model, along with information on the training procedure for each set of weights. For the purposes of this plug-in-and-play example, we comment out a number of available models, for you to load and try out.

# Various pre-trained weights for PhaseNet
# pn_model = sbm.PhaseNet.from_pretrained("ethz")
# pn_model = sbm.PhaseNet.from_pretrained("instance")
# pn_model = sbm.PhaseNet.from_pretrained("scedc")
pn_model = sbm.PhaseNet.from_pretrained("stead")
# pn_model = sbm.PhaseNet.from_pretrained("geofon")
# pn_model = sbm.PhaseNet.from_pretrained("neic")

# Various pre-trained weights for EQT
eqt_model = sbm.EQTransformer.from_pretrained("original")
# eqt_model = sbm.EQTransformer.from_pretrained("ethz")
# eqt_model = sbm.EQTransformer.from_pretrained("instance")
# eqt_model = sbm.EQTransformer.from_pretrained("scedc")
# eqt_model = sbm.EQTransformer.from_pretrained("stead")
# eqt_model = sbm.EQTransformer.from_pretrained("geofon")

# Various pre-trained weights for GPD
# gpd_model = sbm.GPD.from_pretrained("original")
# gpd_model = sbm.GPD.from_pretrained("ethz")
# gpd_model = sbm.GPD.from_pretrained("scedc")
gpd_model = sbm.GPD.from_pretrained("stead")
# gpd_model = sbm.GPD.from_pretrained("geofon")
# gpd_model = sbm.GPD.from_pretrained("neic")

If GPU available, send models to GPU.

pn_model.cuda();
eqt_model.cuda();
gpd_model.cuda();

Read in waveforms

This waveform example contains part of the 2009 L’Acquila earthquake sequence occuring in the Abruzzo region of central Italy. The Mw 6.4 mainshock is located ~2 minutes into the seismic recording, with many afterhshocks occuring directly after. The seismic trace loaded in from obspy spans the ~1 hour period immediately following the mainshock.

from obspy.clients.fdsn import Client
from obspy import UTCDateTime
import matplotlib.pyplot as plt
import obspy as obs
import numpy as np
import matplotlib.pyplot as plt

client = Client("INGV")

t = UTCDateTime(2009, 4, 6, 1, 30)
stream = client.get_waveforms(network="MN", station="AQU", location="*", channel="HH?", starttime=t, endtime=t+3600)


fig = plt.figure(figsize=(20, 5))
ax = fig.add_subplot(111)
for i in range(3):
    ax.plot(stream[i].times(), stream[i].data, label=stream[i].stats.channel)
ax.legend();

Predict model on seismic stream

We can apply our varying models now to the seismic stream through the high level annotate method. This will output the charactersitic functions for the probability of a P-phase, S-phase, and noise class respectively.

pn_preds = pn_model.annotate(stream)
eqt_preds = eqt_model.annotate(stream)
gpd_preds = gpd_model.annotate(stream)

pn_preds

3 Trace(s) in Stream:
MN.AQU..PhaseNet_P | 2009-04-06T01:30:02.498597Z - 2009-04-06T02:29:57.498597Z | 100.0 Hz, 359501 samples
MN.AQU..PhaseNet_S | 2009-04-06T01:30:02.498597Z - 2009-04-06T02:29:57.498597Z | 100.0 Hz, 359501 samples
MN.AQU..PhaseNet_N | 2009-04-06T01:30:02.498597Z - 2009-04-06T02:29:57.498597Z | 100.0 Hz, 359501 samples

Step through to ‘zoom in’ on predictions

Enter window length in seconds, and define which model predictions you would like to plot

wlength = 5 * 60 #@param

color_dict = {"P": "C0", "S": "C1", "Detection": "C2"}

for s in range(0, int(stream[0].stats.endtime - stream[0].stats.starttime), wlength):
    t0 = stream[0].stats.starttime + s
    t1 = t0 + wlength
    subst = stream.slice(t0, t1)

    fig, ax = plt.subplots(4, 1, figsize=(15, 7), sharex=True, gridspec_kw={'hspace' : 0.05, 'height_ratios': [2, 1, 1, 1]})

    for i, preds in enumerate([eqt_preds, pn_preds, gpd_preds]):
        subpreds = preds.slice(t0, t1)
        offset = subpreds[0].stats.starttime - subst[0].stats.starttime
        for pred_trace in subpreds:
            model, pred_class = pred_trace.stats.channel.split("_")
            if pred_class == "N":
                # Skip noise traces
                continue
            c = color_dict[pred_class]
            ax[i + 1].plot(offset + pred_trace.times(), pred_trace.data, label=pred_class, c=c)
        ax[i + 1].set_ylabel(model)
        ax[i + 1].legend(loc=2)
        ax[i + 1].set_ylim(0, 1.1)

    ax[0].plot(subst[-1].times(), subst[-1].data / np.amax(subst[-1].data), 'k', label=subst[-1].stats.channel)
    ax[0].set_xlim(0, wlength)
    ax[0].set_ylabel('Normalised Amplitude')
    ax[3].set_xlabel('Time [s]')
    ax[0].legend(loc=2)
    plt.show()

If you would like to then make deterministic detections from the continuous characteristic functions, the classify method of each model will allow you to perform this task.

outputs = pn_model.classify(stream)

print(outputs)

ClassifyOutput(creator='PhaseNet', picks=PickList with 122 entries:

MN.AQU.	2009-04-06T01:32:42.308597Z	P
MN.AQU.	2009-04-06T01:36:31.508597Z	P
MN.AQU.	2009-04-06T01:36:32.868597Z	S
...
MN.AQU.	2009-04-06T02:27:50.218597Z	S
MN.AQU.	2009-04-06T02:28:53.088597Z	P
MN.AQU.	2009-04-06T02:28:55.708597Z	S)

With each pick storing the following information:

print(outputs.picks[0].__dict__)

{'trace_id': 'MN.AQU.', 'start_time': UTCDateTime(2009, 4, 6, 1, 32, 41, 938597), 'end_time': UTCDateTime(2009, 4, 6, 1, 32, 42, 828597), 'peak_time': UTCDateTime(2009, 4, 6, 1, 32, 42, 308597), 'peak_value': 0.5062704, 'phase': 'P'}