Just wanted to share a little script I made with Obspy, I got sick of having to manually input stuff like dates and lat/lon numbers etc so I made it so you only have to enter your Shake’s station name and a USGS Event ID and it does the rest for you.

Here are some examples although sadly I can’t view local earthquakes since the only Geo Science Australia fdsn catalogue updates once a month

Here is the code for anyone interested:

```
from obspy.clients.fdsn import Client
from obspy.core import UTCDateTime, Stream
from obspy.taup import TauPyModel
import matplotlib.pyplot as plt
rsClient = Client('RASPISHAKE')
usgsClient = Client('USGS')
#---------------------------------------------[Set Station and Event Data]---------------------------------------------#
# Set the station name and download the response information
stn = 'RD371'
inv = rsClient.get_stations(network='AM', station=stn, level='RESP')
latStn = inv[0][0].latitude
lonStn = inv[0][0].longitude
eleStn = inv[0][0].elevation
print('Station:', stn)
print('Latitude:', latStn)
print('Longitude:', lonStn)
print('Elevation:', eleStn)
# Get event data
eventID = "us7000jvl3"
event = usgsClient.get_events(eventid=eventID)
event_time = event[0].origins[0].time
latEq = event[0].origins[0].latitude
lonEq = event[0].origins[0].longitude
depth_km = event[0].origins[0].depth / 1000 # Convert depth to km
desired_magnitude_type = "Mww"
desired_magnitude = None
for magnitude in event[0].magnitudes:
if magnitude.magnitude_type == desired_magnitude_type:
desired_magnitude = magnitude.mag
break
locEq = event[0].event_descriptions[0].text
print('Event ID:', eventID)
print('Event Time:', event_time)
print('Latitude:', latEq)
print('Longitude:', lonEq)
print('Depth:', depth_km,"km")
print('Magnitude:', desired_magnitude)
print('Location:', locEq)
# Set data start/end time
start = event_time
end = event_time + 30 * 60 # 15 minutes = 15 seconds * 60 minutes
# Set channel names
channels = ['EHZ', 'HDF']
# Get waveforms
stream = Stream()
for ch in channels:
trace = rsClient.get_waveforms('AM', stn, '00', ch, start, end)
trace.merge(method=0, fill_value='latest')
trace.detrend(type='demean')
stream += trace
#---------------------------------------------[Calculate distance between station and event]---------------------------------------------#
from math import radians, sin, cos, sqrt, atan2
def calculate_distance(lat1, lon1, lat2, lon2):
# Convert degrees to radians
lat1_rad = radians(lat1)
lon1_rad = radians(lon1)
lat2_rad = radians(lat2)
lon2_rad = radians(lon2)
# Earth's radius in kilometers
radius = 6371
# Difference between latitudes and longitudes
dlat = lat2_rad - lat1_rad
dlon = lon2_rad - lon1_rad
# Haversine formula
a = sin(dlat / 2) ** 2 + cos(lat1_rad) * cos(lat2_rad) * sin(dlon / 2) ** 2
c = 2 * atan2(sqrt(a), sqrt(1 - a))
# Calculate distance
distance_km = radius * c
return distance_km
def convert_km_to_degrees(distance_km):
# Earth's circumference in kilometers
earth_circumference_km = 40075
# Convert distance to degrees
distance_deg = distance_km * (1 / earth_circumference_km) * 360
return distance_deg
# Coordinates
lat1 = latStn
lon1 = lonStn
lat2 = latEq
lon2 = lonEq
# Calculate distance between the coordinates
distance_km = calculate_distance(lat1, lon1, lat2, lon2)
# Convert distance to degrees
distance_deg = convert_km_to_degrees(distance_km)
# Print the distance in kilometers and degrees
print("Distance: {:.2f} km".format(distance_km))
print("Distance: {:.4f} degrees".format(distance_deg))
#---------------------------------------------[Plot]---------------------------------------------#
# Get arrivals
model = TauPyModel(model='iasp91')
arrivals = model.get_travel_times(source_depth_in_km=depth_km, distance_in_degree=distance_deg)
# Set up figure and subplots
fig = plt.figure(figsize=(20,5))
ax = fig.add_subplot(1,1,1)
# plot the trace
ax.plot(stream[0].times(reftime=start), stream[0].data, 'k')
# Plot P and S waves
p_arrival = None
s_arrival = None
for arr in arrivals:
if arr.name == 'P' and p_arrival is None:
p_arrival = arr
elif arr.name == 'S' and s_arrival is None:
s_arrival = arr
# Plot P and S waves
if p_arrival is not None:
p_arrival_time = p_arrival.time + start.timestamp - stream[0].stats.starttime.timestamp
ax.axvline(x=p_arrival_time, color='b', linestyle='-', linewidth=1)
ax.text(p_arrival_time, max(stream[0].data), p_arrival.name, va='bottom', ha='left')
if s_arrival is not None:
s_arrival_time = s_arrival.time + start.timestamp - stream[0].stats.starttime.timestamp
ax.axvline(x=s_arrival_time, color='r', linestyle='-', linewidth=1)
ax.text(s_arrival_time, max(stream[0].data), s_arrival.name, va='bottom', ha='left')
# Set up plot details
ax.set_ylabel('Counts')
ax.set_xlabel('Time (S)')
fig.suptitle("AM."+stn+" M"+str(desired_magnitude)+" "+locEq+" "+eventID+event_time.strftime(' %d/%m/%Y %H:%M:%S UTC')) #Title of the figure
plt.tight_layout()
plt.show()
```

Any feedback would be appreciated thanks!