Any ideas to improve on this code?

[Matt]

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!

[Stormchaser]

Hello Matt,

This is amazing; thank you for sharing your work and code with everyone in the community!

Personal observations are only cosmetic ones:

1. You can try offsetting the P and S texts a bit more so that they are not attached to the line
2. You can think about removing the white spaces at the start/end of the plot with plt.margins(x=0) or ax.margins(x=0)
3. You can also think about filtering the raw data to remove unwanted local spikes and show the events more prominently: Filtering Seismograms — ObsPy 1.4.0 documentation Usually, bandpass filters do the job, and you can take examples from the pre-sets that we have on DataView https://dataview.raspberryshake.org/.

[Matt]

Thanks for the feedback Stormchaser, I haven’t figured out the white spaces yet. Got to mess around with the filter but it currently looks like this
Filtered:

Raw Data:

[Stormchaser]

Yes, that appears to be a bandpass filter highlighting the initial P-wave arrival over everything else.

You can play with them and see how the various filtering values can be applied to earthquakes at different distances to make their waves more evident.

[sheeny]

G’Day Matt,

You might like to have a look at this post to ensure you have the latest epoch inventory:

Station Lat Long not matching Map Position

OK, I’ve come up with a solution, but if anyone knows a simpler or better way, I’d appreciate them sharing. This is the code I’ve come up with to get the latest epoch data in the inventory: inv = rs.get_stations(network='AM', station=slist[i], level='RESP') # get the instrument response #print(inv[0]) k=0 while True: sta = inv[0][k] #station metadata staOK = sta.is_active(time=eventTime) if staOK: break …

It may not be a problem to you now, but if you move your shake or want to use your code on another shake, it might help.

Now to update my old report…

Al.