Forensic Seismology and Acoustics¶
Notebook for class on February 17th.
Stephen Arrowsmith (sarrowsmith@smu.edu)
Data Exercise 1: Reading and analyzing DPRK data¶
This first class introduces you to analyzing seismic data with Python. We're going to dive in the deep end and try to tread water. First thing we'll do is import some 'modules'. These are different code packages that contain different things. We can import whole modules, parts of modules, and rename modules when we import them.
The module called 'forensic' is built for this class. The module called 'obspy' is an important one for analyzing seismic and acoustic data.
import forensic
from obspy import UTCDateTime
from pyproj import Geod
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import datetime, pickle
The next thing we'll do is define an event location and time. We're just defining the location and time of the DPRK nuclear test in September 2017, using the values from the Myers et al. (2018) paper that we read already. Defining times can be tricky, ObsPy (the main seismic data module we read in above) likes times in a type of time format called 'UTCDateTime'.
evla = 41.3; evlo = 129.08; evdp = 0.08
evt0 = UTCDateTime(datetime.datetime(2017,9,3,3,30,1))
evt0
Reading data¶
You can read data seismic directly into Python with an internet connection using the forensic module like below. We'll skip over the details of how to do this for now, because it's complex, but it might come in handy later. For now, we'll open the file called 'DPRK.p' that you downloaded already.
#st = forensic.getStreamFromIRIS(evt0-600, evt0+3000, '*', 'IU,II,IC', 'BHZ')
#pickle.dump(st, open('DPRK.p', 'wb'))
st = pickle.load(open('DPRK.p', 'rb'))
Exploring ObsPy¶
To explore data in ObsPy, you need to understand how the data are structured in an 'object' called an ObsPy Stream: https://docs.obspy.org/packages/obspy.core.html#waveform-data
print(st)
tr = st[0]
print(tr.stats)
tr.data
%pylab notebook
plt.plot(tr.data)
%pylab notebook
#fig, ax = plt.subplots()
fig = tr.plot(handle=True)
for tr in st:
print(tr.stats.station, tr.stats.network,
tr.stats.channel, tr.stats.sampling_rate,
round(forensic.get_dist(tr, evla, evlo)/111.1,2))
%pylab notebook
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
for tr in st:
ax.plot(tr.stats.sac.stlo, tr.stats.sac.stla, 'y^', transform=ccrs.Geodetic())
ax.plot(evlo, evla, 'r*', transform=ccrs.Geodetic())
ax.gridlines(draw_labels=True)
plt.show()
Exploring the DPRK data¶
We'll plot 'record sections' which are plots that show multiple seismograms, arranged according to the location and origin time of an event (in this case, the September 2017 DPRK test).
%pylab notebook
forensic.plot_recordsection_global(st, evlo, evla,
date2num(evt0.datetime),
scaleFactor=1.5)
plt.xlim([0., 160.])
plt.ylim([0., 50.])
%pylab notebook
st_f = st.copy()
st_f.taper(type='cosine', max_percentage=0.001)
st_f.filter('highpass', freq=1)
forensic.plot_recordsection_global(st_f, evlo, evla,
date2num(evt0.datetime),
scaleFactor=1.5)
plt.xlim([0., 160.])
plt.ylim([0., 50.])