In this very short introduction to ObsPy we will see how we can start from scratch to acquire time series data (aka waveforms) of the ground shaking resulting from earthquakes, recorded at seismometer stations. Many global seismometer recordings are free to be downloaded by everybody. We will also see how these data can be handled in Python using the ObsPy framework.
Again, please execute the following cell, it contains a few adjustments to the notebook.
%matplotlib inline from __future__ import print_function import numpy as np import matplotlib.pyplot as plt plt.style.use('ggplot') plt.rcParams['figure.figsize'] = 12, 8
the most important protocol are FDSN web services (other protocols work very similar)
<img src="images/fdsn_short_urls.png" width=60%>
get_events()" method has many ways to restrict the search results (time, geographic, magnitude, ..)
from obspy.fdsn import Client client = Client("IRIS") catalog = client.get_events(minmagnitude=8.5, starttime="2010-01-01") print(catalog) catalog.plot();
Catalogobject, which is a collection (~list) of
Eventobjects are again collections of other resources (origins, magnitudes, picks, ...)
event = catalog print(event)
# you can use Tab-Completion to see what attributes/methods the event object has: event. origin.
from obspy import readEvents catalog2 = readEvents("./data/events_unterhaching.xml") print(catalog2) # these events contain picks, too print(catalog2.picks)
get_stations()" method has many ways to restrict the search results (time, geographic, station names, ..)
t = event.origins.time print(t) inventory = client.get_stations( network="TA", starttime=t, endtime=t+10) print(inventory) inventory.plot(projection="local");
Inventoryobject, which is a collection (~list) of
Networkobjects are again collections of
Stationobjects, which are collections of
<img src="images/Inventory.svg" width=100%> <img src="images/stationxml_light.png" width=100%>
# let's request full detail metadata for one particular station inventory = client.get_stations(network="TA", station="S36A", level="response") # when searching for many available stations response information is not included by default, # because of the huge size. we explicitly need to specify that we want response information included print(inventory)
network = inventory print(network)
station = network print(station)
station = station.select(channel="BH*") print(station)
channel = station print(channel) print(channel.response)
# you can use Tab-Completion to see what attributes/methods the inventory object has: inventory. channel.
format="STATIONXML"in that case!)
from obspy import read_inventory inventory2 = read_inventory("./data/station_uh1.xml", format="STATIONXML") print(inventory2)
get_waveforms()" method needs the unique identification of the station (network, station, location and channel codes) and a time range (start time, end time) of requested data
stream = client.get_waveforms(network="TA", station="S36A", location="*", channel="BH*", starttime=t+10*60, endtime=t+70*60) # for functions/methods that need a fixed set of parameters, # we usually omit the parameter names and specify them in the expected order # Note that new timestamp objects can be created by # adding/subtracting seconds to/from an existing timestamp object. # (for details search the ObsPy documentation pages for "UTCDateTime") stream = client.get_waveforms("TA", "S36A", "*", "BH*", t+10*60, t+70*60) print(stream) stream.plot()
the nested ObsPy Inventory class structure (Inventory/Station/Channel/Response/...) is closely modelled after FDSN StationXML (the international standard exchange format for station metadata)
Streamobject, which is a collection (~list) of
Traceobjects consist of one single, contiguous waveform data block (i.e. regularly spaced time series, no gaps)
Traceobjects consist of two attributes:
trace.data(the raw time series as an
trace.stats(metainformation needed to interpret the raw sample data)
<img src="images/Stream_Trace.svg" width=90%>
trace = stream # trace.data is the numpy array with the raw samples print(trace.data) print(trace.data[20:30]) # arithmetic operations on the data print(trace.data[20:30] / 223.5) # using attached numpy methods print(trace.data.mean()) # pointers to the array in memory can be passed to C/Fortran routines from inside Python!
# trace.stats contains the basic metadata identifying the trace print(trace) print(trace.stats) print(trace.stats.sampling_rate)
# many convenience routines are attahed to Stream/Trace, use Tab completion stream.
from obspy import read stream2 = read("./data/waveforms_uh1_eh?.mseed") print(stream2)
from obspy import UTCDateTime t1 = UTCDateTime("2015-08-04T10:30") print(t1.hour) print(t1.strftime("%c"))
+/-operator to add/subtract seconds
-operator to get diff of two timestamps in seconds
# calculate time since Tohoku earthquake in years now = UTCDateTime() print((now - t) / (3600 * 24 * 365)) # 2 hours from now print(now + 2 * 3600)
# first, let's have a look at this channel's instrument response channel.plot(min_freq=1e-4);
# we already have fetched the station metadata, # including every piece of information down to the instrument response # let's make on a copy of the raw data, to be able to compare afterwards stream_corrected = stream.copy() # attach the instrument response to the waveforms.. stream_corrected.attach_response(inventory) # and convert the data to ground velocity in m/s, # taking into account the specifics of the recording system stream_corrected.remove_response(water_level=10, output="VEL") stream.plot() stream_corrected.plot()
Inventory/... for easy manipulation
<img src="images/collage.png" width=100%>
This example prepares instrument corrected waveforms, 30 seconds around expected P arrival (buffer for processing!) for an epicentral range of 30-40 degrees for any TA station/earthquake combination with magnitude larger 7.
from obspy.taup import TauPyModel from obspy.core.util.geodetics import gps2DistAzimuth, kilometer2degrees model = TauPyModel(model="iasp91") catalog = client.get_events(starttime="2009-001", endtime="2012-001", minmagnitude=7) print(catalog) network = "TA" minradius = 30 maxradius = 40 phase = "P" for event in catalog: origin = event.origins try: inventory = client.get_stations(network=network, startbefore=origin.time, endafter=origin.time, longitude=origin.longitude, latitude=origin.latitude, minradius=minradius, maxradius=maxradius) except: continue print(event) for station in inventory[:3]: distance, _, _ = gps2DistAzimuth(origin.latitude, origin.longitude, station.latitude, station.longitude) distance = kilometer2degrees(distance / 1e3) arrivals = model.get_travel_times(origin.depth / 1e3, distance, phase_list=[phase]) traveltime = arrivals.time arrival_time = origin.time + traveltime st = client.get_waveforms(network, station.code, "*", "BH*", arrival_time - 200, arrival_time + 200, attach_response=True) st.remove_response(water_level=10, output="VEL") st.filter("lowpass", freq=1) st.trim(arrival_time - 10, arrival_time + 20) st.plot() break
b) Take one event out of the catalog. Print its summary. You can see, there is information included on picks set for the event. Print the information of at least one pick that is included in the event metadata. We now know one station that recorded the event.
c) Save the origin time of the event in a new variable. Download waveform data for the station around the time of the event (e.g. 10 seconds before and 20 seconds after), connecting to the FDSN server at the observatory at "
http://erde.geophysik.uni-muenchen.de" (or read from local file
./data/waveforms_uh1_eh.mseed). Put a "
*" wildcard as the third (and last) letter of the channel code to download all three components (vertical, north, east). Plot the waveform data.
d) Download station metadata for this station including the instrument response information (using the same client, or read it from file
./data/station_uh1.xml). Attach the response information to the waveforms and remove the instrument response from the waveforms. Plot the waveforms again (now in ground velocity in m/s).
e) The peak ground velocity (PGV) is an important measure to judge possible effects on buildings. Determine the maximum value (absolute, whether positive or negative!) of either of the two horizontal components' data arrays ("N", "E"). You can use Python's builtin "
max()" and "
abs()" functions. For your information, ground shaking (in the frequency ranges of these very close earthquakes) can become perceptible to humans at above 0.5-1 mm/s, damage to buildings can be safely excluded for PGV values not exceeding 3 mm/s.