I’m trying to plot arrival time curves onto a section plot, and not having any success. I know other’s have done it, but I’ve exhausted all my ideas for now.
My idea was to make a crude program for finding small local quakes that aren’t on ShakeNet or other sources. I’ve observed a few interesting events like this now. The idea is to adjust the time and location of the quake until the arrival time curves match the waveform plots.
Here’s what the report looks like ATM with the plot_travel_times line commented out.
I’m struggling to get the plot_travel_times line to work.
This is my code:
# -*- coding: utf-8 -*-
"""
Created on Sun May 21 12:22:53 2023
@author: al72
"""
from obspy.clients.fdsn import Client
from obspy.core import UTCDateTime, Stream
#from obspy.signal import filter
import matplotlib.pyplot as plt
#from matplotlib.ticker import AutoMinorLocator
#import numpy as np
from obspy.taup import TauPyModel, plot_travel_times
#import math
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from matplotlib.transforms import blended_transform_factory
from obspy.geodetics import gps2dist_azimuth, kilometers2degrees, degrees2kilometers
rs = Client('https://data.raspberryshake.org/')
def stationList(sl):
sl += ['R21C0']
sl += ['R811A']
sl += ['R9AF3']
sl += ['R571C']
#sl += ['R6D2A']
#sl += ['RF35D']
#sl += ['R7AF5']
#sl += ['R3756']
#sl += ['R6707']
#sl += ['RB18D']
#print(sl)
return sl
def freqFilter(ft):
# bandpass filter - select to suit system noise and range of quake
#ft = [0.1, 0.1, 0.8, 0.9]
#ft = [0.3, 0.3, 0.8, 0.9]
#ft = [0.5, 0.5, 2, 2.1]
ft = [0.7, 0.7, 2, 2.1] #distant quake
#ft = [0.7, 0.7, 3, 3.1]
#ft = [0.7, 0.7, 4, 4.1]
#ft = [0.7, 0.7, 6, 6.1]
#ft = [0.7, 0.7, 8, 8.1]
#ft = [1, 1, 10, 10.1]
#ft = [1, 1, 20, 20.1]
#ft = [3, 3, 20, 20.1] #use for local quakes
#print(ft)
return ft
def buildStream(strm):
n = len(slist)
print(n)
for i in range(0, n):
inv = rs.get_stations(network='AM', station=slist[i], level='RESP') # get the instrument response
#print(inv[0])
#read each epoch to find the one that's active for the event
k=0
while True:
sta = inv[0][k] #station metadata
if sta.is_active(time=eventTime):
break
k += 1
latS = sta.latitude #active station latitude
lonS = sta.longitude #active station longitude
#eleS = sta.elevation #active station elevation
#print(slist[i], latS, lonS)
#print(sta)
trace = rs.get_waveforms('AM', station=slist[i], location = "00", channel = 'EHZ', starttime = start, endtime = end)
trace.merge(method=0, fill_value='latest') #fill in any gaps in the data to prevent a crash
trace.detrend(type='demean') #detrend the data
tr1 = trace.remove_response(inventory=inv,zero_mean=True,pre_filt=filt,output='VEL',water_level=.001, plot=False) # convert to Vel
tr1[0].stats.distance = gps2dist_azimuth(latS, lonS, latE, lonE)[0]
tr1[0].stats.latitude = latS
tr1[0].stats.longitude = lonS
tr1[0].stats.colour = colours[i]
strm += tr1
#strm.plot(method='full', equal_scale=False)
return strm
def k2d(x):
return kilometers2degrees(x)
def d2k(x):
return degrees2kilometers(x)
eventTime = UTCDateTime(2023, 5, 13, 5, 48, 00) # (YYYY, m, d, H, M, S) **** Enter data****
latE = -32.3 # quake latitude + N -S **** Enter data****
lonE = 149.8 # quake longitude + E - W **** Enter data****
depth = 1 # quake depth, km **** Enter data****
mag = 1 # quake magnitude **** Enter data****
eventID = 'unknown' # ID for the event **** Enter data****
locE = "Mudgeeish, NSW, Australia" # location name **** Enter data****
slist = []
filt = []
stationList(slist)
freqFilter(filt)
colours = ['r', 'b', 'g', 'k', 'c', 'm', 'purple', 'orange', 'gold', 'midnightblue']
#print(slist)
#print(filt)
#set up the plot
delay = 0
duration = 80
start = eventTime + delay
end = start + duration
st = Stream()
buildStream(st)
#print(st)
#set up the figure
fig = plt.figure(figsize=(20,14), dpi=100) # set to page size in inches
#build the section plot
ax1 = fig.add_subplot(1,2,1)
st.plot(type='section', plot_dx=50e3, recordlength=80, time_down=True, linewidth=.25, grid_linewidth=.25, show=False, fig=fig)
# Plot customization: Add station labels to offset axis
ax = ax1.axes
transform = blended_transform_factory(ax.transData, ax.transAxes)
for t in st:
ax.text(t.stats.distance / 1e3, 1.0, t.stats.station, rotation=90,
va="bottom", ha="center", color=t.stats.colour, transform=transform, zorder=10)
#setup secondary x axis
secax_x1 = ax1.secondary_xaxis('top', functions = (k2d, d2k))
#plot arrivals times
model = TauPyModel(model="iasp91")
#ax = plot_travel_times(source_depth=1, ax=ax, fig=fig, max_degrees=1.6, phase_list=['P', 'S'])
#plot the map
ax2 = fig.add_subplot(1,2,2, projection=ccrs.PlateCarree())
ax2.set_extent([148,152,-35,-30], crs=ccrs.PlateCarree())
#ax2.coastlines(resolution='110m')
#ax2.stock_img()
# Create a features
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='50m',
facecolor='none')
ax2.add_feature(cfeature.LAND)
ax2.add_feature(cfeature.OCEAN)
ax2.add_feature(cfeature.COASTLINE)
ax2.add_feature(states_provinces, edgecolor='gray')
ax2.add_feature(cfeature.LAKES, alpha=0.5)
ax2.add_feature(cfeature.RIVERS)
ax2.gridlines()
#plot event/earthquake position on map
ax2.plot(lonE, latE,
color='yellow', marker='*', markersize=20, markeredgecolor='black',
transform=ccrs.Geodetic(),
)
#plot station positions on map
for tr in st:
ax2.plot(tr.stats.longitude, tr.stats.latitude,
color=tr.stats.colour, marker='H', markersize=12, markeredgecolor='black',
transform=ccrs.Geodetic(),
)
ax2.plot
plt.show()
Line 134 is the problem. I’ve tried every permutation I can think of, so hoping some else has done this and can point me in the right direction.
Thanks,
Al.
