The following code snippets go over some of the fundamental features of FDSN-via-ObsPy

Requirements :

-A working internet connection

-Python3 (2.7 may also work but not recommended) & and interactive interpreter (e.g. ipython3)

-ObsPy installed. Use a package manager (such as conda or pip) to install all dependencies (numpy, SQLAlchemy, lxml, Cython, future, proj4, matplotlib and cartopy or basemap). For detailed instructions to install ObsPy, go to https://www.obspy.org/.

At its simplest, obspy can be used in a manner equivalent to curl or wget (or a web browser!) to grab miniseed data via a single target, similar to pointing to a file

Element	Signification
?	marks the end of the website adress serving the data
&	separates the arguments to select AusPass data
http://auspass.edu.au/	Client name (AusPass website serving FDSN data)
fdsnws/	FSDN Web Services on that client
dataselect/1/	Data retrieval services (distinct for station metadata and event metadata)
query	Do a download
net=S1	Seismic network "S1" - AUstralian Seismometers In Schools - AUSIS
sta=AUAYR	Station "AUAYR" located at Ayr, Queensland, Australia
cha=BH?	All broadband channels available starting with "BH" ie. BHE, BHN and BHZ
start=2018-03-29T20:00:00.000	Data start time using ISO standard date and time format
end=2018-03-29T20:10:00.000	Data end time using ISO standard date and time format

# a very simple version which grabs data using the above elements pasted together

import obspy
from matplotlib import pyplot as plt
st = obspy.read('http://auspass.edu.au/fdsnws/dataselect/1/query?net=S1&sta=AUAYR&cha=BH?&start=2018-03-29T20:00:00.000&end=2018-03-29T20:10:00.000')
print(st)
st.plot(); plt.close() #note that "plt.close()" may not be necessary for you! 

3 Trace(s) in Stream:
S1.AUAYR..BHZ | 2018-03-29T19:59:41.000000Z - 2018-03-29T20:10:01.500000Z | 10.0 Hz, 6206 samples
S1.AUAYR..BHE | 2018-03-29T19:59:48.000000Z - 2018-03-29T20:10:13.800000Z | 10.0 Hz, 6259 samples
S1.AUAYR..BHN | 2018-03-29T19:59:57.800000Z - 2018-03-29T20:10:00.900000Z | 10.0 Hz, 6032 samples

However the true utility of this protocol is its use in a modular or object-oriented sense. Let's try again.

Establish an FDSN connection

• User_agent field is not required, but it is polite to include and allows us to contact you in the unlikely event your requests are causing an issue. Your address is not stored or added to any lists.

• User and password should ONLY be given here if you are requesting restricted data. Restricted data is typically limited to the PI for up to 2 years following its availability here.

import obspy
from obspy.clients.fdsn import Client as FDSN_Client
from matplotlib import pyplot as plt
auspass = FDSN_Client('http://auspass.edu.au',user_agent='email@address.com') #,user='username',password='password')

Set request parameters

network = "S1"
station = "AUAYR"
location = "*" #n.b. AusPass data typically has a univeral blank (e.g. '') value for ALL station locations. "*" works just as well. 
channel = "BH*" #note that BH* channels here are 10 Hz, HH* are 100 
time0 = obspy.UTCDateTime(2018,3,29,1,19,5) #1 hour, 19 minutes, and 5 seconds after start of March 29 2018 (UTC) 
time1 = time0 + 180 #180 seconds after time0

Pull data from AusPass, and plot

• The commonly used "st" variable is short for "stream", comprising a group of ~similar waveforms clustered together
• A "trace" (usually "tr") is one element of a stream array e.g. tr = st[0]

st = auspass.get_waveforms(network=network,station=station,location=location,channel=channel,starttime=time0,endtime=time1)
print(st)
st.plot(); plt.close()

3 Trace(s) in Stream:
S1.AUAYR..BHZ | 2018-03-29T01:19:05.000000Z - 2018-03-29T01:22:05.000000Z | 10.0 Hz, 1801 samples
S1.AUAYR..BHE | 2018-03-29T01:19:05.000000Z - 2018-03-29T01:22:05.000000Z | 10.0 Hz, 1801 samples
S1.AUAYR..BHN | 2018-03-29T01:19:05.000000Z - 2018-03-29T01:22:05.000000Z | 10.0 Hz, 1801 samples

Network and staton metadata (and response information) can be downloaded from AusPass in a similar manner.

• This is commonly referred to as "inventory" data (with variable "inv")

inv = auspass.get_stations(network=network,level='channel') #this downloads all station information for a network.
# NOTE--- if you ALSO need response information, set level='response'. By default this is assumed to be 'channel'

#inventory/response can also be loaded via local file, or alternate formats such as RESP or DATALESS SEED (among others)
# inv = obspy.read_inventory("filename",format='SEED') #for DATALESS SEED, or format='RESP'

# # You can also search for stations/networks within a certain area, or filter by time with starttime= and endtime=

#bounding box
# inv = auspass.get_stations(network="*",minlatitude=-23.0, maxlatitude=-22.5, minlongitude=133.0, maxlongitude=134.0)

#or even search for stations/networks within a min/max radial distance (in degrees)
# inv = auspass.get_stations(network="*",latitude=-23.0, latitude=-22.5, minradius=0.1,maxradius=4)

# And then can plot the stations on a map
inv.plot(projection='local'); plt.close() #can quickly show where these stations are

• Like most python codes, references for these functions are accessible just by typing "auspass.get_stations?" or "auspass.get_waveforms?"

Finally let's work with this data!

Filtering is easy, and there are many types and methods to choose from in ObsPy.

# Here is a simple Butterworth Bandpass between 2-3 Hz
st.filter('bandpass',freqmin=0.2, freqmax=4.0)
st.plot(); plt.close()

NOTE that this operation is irreversible! So if you want to try a different filter you'll have to re-load your waveform data.

st = auspass.get_waveforms(network=network,station=station,location=location,channel=channel,starttime=time0,endtime=time1)
st_bak = st.copy()
st.filter('bandpass',freqmin=1.0, freqmax=4.0)
st = st_bak.copy() #revert to original 

Removing instrument response is also very very easy!

#first, be sure that response information exists within our inventory. this is done by setting level='response'
inv = auspass.get_stations(network=network,level='response')
st.remove_response(inventory=inv) #this automatically parses the inventory for the correct network/station/channel response!
st.plot(); plt.show()

Response information can also be examined and even plotted

#response is at the channel level (three subsets down) e.g. inv[network][station][channel]
inv0 = inv.select(network=network,station=station)
response = inv0[0][0][0].response
response.plot(min_freq=.001); plt.close()

#you can also directly compare the effect of removing the response on a waveform
st = st_bak.copy()
tr = st[0]
tr.remove_response(inventory=inv, output="DISP",pre_filt=[0.01,.1,3,4],water_level=60, plot=True); plt.show(); plt.close()

No handles with labels found to put in legend.

Downloading EVENTS

• AusPass does not host independent earthquake information, but you can download it from many different agencies. See https://docs.obspy.org/packages/obspy.clients.fdsn.html

#let's download earthquakes from USGS...

usgs = FDSN_Client("USGS") #establish connection in the same way as auspass. n.b."USGS" has its own shortcut, no address needed

#let's just download ALL earthquakes in the southern hemisphere during January 2016 greater than M6
events = usgs.get_events(starttime=obspy.UTCDateTime(2016,1,1), endtime=obspy.UTCDateTime(2016,2,1), maxlatitude=0, minmagnitude=6)

#show results & plot
print(events)
events.plot(); plt.close()

5 Event(s) in Catalog:
2016-01-31T17:38:59.950000Z | -63.263, +169.149 | 6.1 mww | manual
2016-01-26T03:10:20.750000Z |  -5.295, +153.245 | 6.1 mww | manual
2016-01-14T03:25:28.270000Z | -19.760,  -63.329 | 6.1 mww | manual
2016-01-05T09:34:14.620000Z | -54.291, -136.260 | 6.0 mww | manual
2016-01-01T02:00:39.950000Z | -50.557, +139.449 | 6.3 mww | manual

Some other common operations...

st.trim(starttime=time0,endtime=time1-5) #trims the last 5 seconds off all data in the stream

sub_st = st.slice(starttime=time0,endtime=time1-5) #returns a NEW stream "sub_st" instead, without altering original

st.decimate(factor=2) #performs a literal decimation bt "factor" amount (e.g. 100 Hz > 50 Hz)

##interpolate data to a set sampling rate (or via number or samples, or a time range, or...? 
#Performs better than straight decimation when the factor is large)
st.interpolate(method='quadratic',sampling_rate=10) 

3 Trace(s) in Stream:
S1.AUAYR..BHZ | 2018-03-29T01:19:05.000000Z - 2018-03-29T01:22:00.000000Z | 10.0 Hz, 1751 samples
S1.AUAYR..BHE | 2018-03-29T01:19:05.000000Z - 2018-03-29T01:22:00.000000Z | 10.0 Hz, 1751 samples
S1.AUAYR..BHN | 2018-03-29T01:19:05.000000Z - 2018-03-29T01:22:00.000000Z | 10.0 Hz, 1751 samples

• The waveform data contained in each trace is equivalent to a numpy array and all operations in numpy can be performed (e.g. calculating a mean or standard deviation)

import numpy as np
tr = st[0]
std = np.std(tr.data)
mean = np.mean(tr.data)

• You may also want to change the data type (in numpy, "dtype") to a higher precision float, or to an integer. For encoding in Steim2, an int32 is required.

#this loops through all traces in a stream and converts the data type
for tr in st: tr.data = tr.data.astype('int32') #or float64, or int8, or float16, or,..

Metadata can also be called/printed/altered

• It is kept at the "trace" level under "stats"

tr.stats #quickly prints what's available
sr = tr.stats.sampling_rate
delta = tr.stats.delta #the time in seconds per sample (e.g. 1/sr)
start = tr.stats.starttime + 10 #shifts starttime (as well as all subsequent time data) 10 seconds into the future
tr.stats.station = "NEWSTATIONNAME" #changes station name to whatever you like!
tr.stats.location = '00'
print("station name is now %s, and it now starts at %s" % (tr.stats.station,start))

station name is now NEWSTATIONNAME, and it now starts at 2018-03-29T01:19:15.000000Z

Saving to disk

#saves a local copy of inv as a 'STATIONXML', and validates it is structured correctly.
inv.write("S1.xml",format='STATIONXML',validate=True)

#lets save the events we downloaded as well
events.write("events.xml",format="QUAKEML")

#Saving miniseed data can be a little trickier, and you may have to organise filenames and data depending on what you're trying to do.
#To simply save the entire stream of data into a singular file, 
st.write("filename.ms",format='MSEED',encoding='STEIM2',reclen=4096) #note that STEIM2 requires int32 dtype