This demo introduces some simple code that requests data using GeoNet's FDSN webservices and the obspy module in python. This notebook uses Python 3.
In [71]:
from obspy import UTCDateTime
from obspy.clients.fdsn import Client as FDSN_Client
from obspy import read_inventory
In [72]:
client = FDSN_Client("GEONET")
This example is to demonstrate how you can use FDSN webservices to request data for a particular event. Unlike our other notebooks that demonstrate the features of each individual FDSN service this will demonstrate using multiple services to make this common request.
First, get the event information for the earthquake. In this case, we have used the event PublicID. It can be found on the GeoNet website or through QuakeSearch.
This example will demonstrate how to get additional information about the Kaikoura Earthquake
In [73]:
cat = client.get_events(eventid="2016p858000")
print(cat)
To request the waveforms from the webservice we can't use the eventid, instead we need the start and end times of the period we are interested in. This can be defined in many ways but we are going to use a window around the origin time of the event.
In [74]:
event = cat[0]
origin = event.origins[0]
otime = origin.time
print(otime)
Or we can do the above call all in one line.
In [75]:
otime = cat[0].origins[0].time
print(otime)
We probably don't want to request all available GeoNet stations. We have a lot of data! There are a few different ways to do this. Here are a couple of options.
You can use the station service to get a list of stations for which you want to request data.
This example demonstrates a request for high rate (200Hz) weak motion stations within 1 degree of the earthquake.
In [76]:
inventory = client.get_stations(latitude=cat[0].origins[0].latitude,
longitude=cat[0].origins[0].longitude,
maxradius=1.0,
channel="HHZ",
level="channel",
starttime = otime-300,
endtime = otime+600)
print(inventory)
In [77]:
from obspy import Stream
st = Stream()
for network in inventory:
for station in network:
try:
st += client.get_waveforms(network.code, station.code, "*", "HHZ",
otime-300, otime + 600)
except:
pass
print(st)
st.plot()
In [78]:
st = Stream()
for p in range(len(event.picks)):
for a in range(len(origin.arrivals)):
if event.picks[p].resource_id == origin.arrivals[a].pick_id:
if origin.arrivals[a].phase == 'P' and origin.arrivals[a].distance <= 1.0:
try:
st += client.get_waveforms("NZ", event.picks[p].waveform_id['station_code'], "*", "HHZ",
otime-300, otime + 600)
except:
pass
print(st)
st.plot()