This notebook goes with the blog post of 8 September 2017. We'll look at 3 non-earthquakes:
Plus one earthquake — the estimated M8.2 off Mexico this morning, 8 September 2017.
There was a landslide in Greenland in June. It triggered a fatal tsunami (video). When I saw Jascha Polet's tweet about it, I was curious. Ari Hartikainen's later tweet really made me want to look at the data.
Now let's get started!
First some preliminaries to make sure we get inline plots, and to turn off some matplotlib 2.0 warnings about a call ObsPy is making in its plotting functions somewhere.
In [1]:
import warnings
warnings.filterwarnings("ignore")
%matplotlib inline
In [2]:
from obspy.clients.fdsn import Client
client = Client("GFZ")
I'm just following the instructions in the ObsPy documentation.
The event was at about 2340 UTC, but I'm going to grab a 2 hours' worth of data, starting well before the landslide event. You'll see why in a minute.
In [3]:
from obspy import UTCDateTime
t = UTCDateTime("2017-06-17_22:00:00")
st = client.get_waveforms("DK", "NUUG", "*", "HHZ", t, t + 2*60*60)
st.plot()
Notice the landslide event, followed by the periodic signal of the tsunami.
What's HHZ? Seismic stations have lots of instruments and make their waveforms availables as various different channels, depending on the station. For example:
LHZ: long period ~1 Hz, with Z for vertical, or N or E for horizontal.BHZ: broadband 10–80 Hz, with Z for vertical, or N or E for horizontal.HHZ: broadband 80–250 Hz, with Z for vertical, or N or E for horizontal.
In [12]:
import matplotlib.pyplot as plt
inventory = client.get_stations(network="DK", station="NUUG")
inventory.plot(projection='local', color="0.05") # Can be global or ortho or local (at low res by default)
plt.show()
In [13]:
inventory
Out[13]:
In [14]:
st
Out[14]:
In [15]:
z = st.traces[0].data
NFFT = 2048
Fs = st.traces[0].meta.sampling_rate # Sample rate in Hz
label = "{network}.{station}.{location}.{channel}".format(**st.traces[0].meta)
plt.figure(figsize=(15, 5))
Pxx, freqs, bins, im = plt.specgram(z, NFFT=NFFT, Fs=Fs, noverlap=int(0.9*NFFT), cmap='viridis', vmin=-30, vmax=30)
plt.text(100, 45, label, size=15, color='white')
plt.show()
Notice the increasingly frequent signals preceeding the landslide event. You can hear them in the sonification of this event — and listen for the tsunami after the crash.
Let's move around the world to a different kind of event. North Korea made another nuclear test this morning just after 0330 UTC. Inspired by Steven Gibbons's tweets the morning of the test, I wanted to download the data for myself.
From his tweets, we want station MDJ in eastern China, which is in network code IC.
In [16]:
client = Client("IRIS")
inventory = client.get_stations(network="IC", station="MDJ")
inventory.plot(color="0.05")
plt.show()
In [17]:
t = UTCDateTime("2017-09-03_03:30:00")
st = client.get_waveforms("IC", "MDJ", "00", "BHZ", t, t + 5*60)
st.plot()
Again, following Steve's tweets, let's compare to the test last year, on 9 September 2016 at 0030 UTC:
In [18]:
t = UTCDateTime("2016-09-09_00:30:00")
st1 = client.get_waveforms("IC", "MDJ", "00", "BHZ", t, t + 5*60)
st1.plot()
In [19]:
label = "{network}.{station}.{location}.{channel}".format(**st.traces[0].meta)
plt.figure(figsize=(15,5))
plt.plot(st.traces[0], label="3 Sept 2017")
plt.plot(st1.traces[0], label="9 Sept 2016")
plt.text(200, 3e6, label, size=20)
plt.legend()
plt.grid()
plt.show()
I'll leave making the spectrogram for this event as an exercise for the reader. Hint: use the HHZ channel for better time resolution, and get a longer signal for better frequency resolution. See if you can see evidence for other events (chamber collapse? Landslides? about 10 minutes after the detonation.
This time we'll get inspired by @seismo_steve's tweets about Irma. He used stations on Barduda, Guadalupe, and St Maarten; we'll find the one in Barbuda, in the USGS Caribbean CU network.
In [20]:
client = Client("IRIS")
inventory = client.get_stations(network="CU", station="ANWB")
inventory.plot(projection='global', color=0.15) # Can be global or ortho or local (at low res by default_)
plt.show()
In [21]:
t0 = UTCDateTime("2017-09-04_18:00:00")
t1 = t0 + 36*60*60
st = client.get_waveforms("CU", "ANWB", "00", "LH2", t0, t1)
st.plot()
In [22]:
st.traces[0].meta.sampling_rate
Out[22]:
The LH2 channel is long period, with a 1 Hz sample rate. Let's clean it up a bit by eliminating the low frequency wiggle.
In [23]:
st.filter('highpass', freq=0.1)
Out[23]:
In [24]:
st.plot(endtime=t1)
This compares well with Steve's tweet.
You can clearly see the increase in noise as the storm approaches, and the failure of communications (or some other critical part of the system) and loss of signal.
Let's try listening to it!
In [25]:
from scipy.io import wavfile
import numpy as np
N = 24000 # Samples per second, a speed up of 24,000 x
wavfile.write('../data/irma.wav', N, st.traces[0].data.astype(np.int16))
In [26]:
import IPython
IPython.display.Audio("../data/irma.wav")
Out[26]:
OK, that's pretty cool.
OK then, let's look at an actual earthquake. This morning there was a very large earthquake off Mexico, killing at least 15 people.
In [27]:
client = Client("IRIS")
inventory = client.get_stations(network="MX")
inventory.plot(projection='local', color=0.15) # Can be global or ortho or local (at low res by default_)
plt.show()
According to the BBC article, the quake struck at 04:50 UTC.
In [28]:
t = UTCDateTime("2017-09-08_04:45:00")
st = client.get_waveforms("MX", "*", "*", "BHZ", t, t + 20*60)
st.plot(color='blue')
The Mars data is now open to the public so I thought I'd add a way to access the data from IRIS.
In [29]:
client = Client("IRIS")
t = UTCDateTime("2019-02-11_00:33:22")
st = client.get_waveforms("XB", "ELYSE", "00", "HHU", t, t + 10*60)
st.plot(color='blue')
© 2017 agilescientific.com and licensed CC-BY 4.0