I’m working on a Python Report to show filtered Displacement, Velocity and Acceleration waveforms, the velocity spectrogram and velocity PSD as well as spherical rays and the phase arrival times.
I’ve got the basic parts of the graphs done for D, V and A and the spectrogram and PSD, but I need help with plotting the spherical rays and plotting the arrival times on the D,V and A plots. I don’t have problem producing a stand alone spherical ray plot, but I can’t work out how to get it into a subplot.
Below is the code so far, including commented out code that either I’ve used for development and testing or has simply failed and I don’t know how to fix ATM. I’ll also post the output so far FYI. The RH half of the page is where I was the spherical ray plot.
I would appreciate any tips and pointers in the right direction, or any example code anyone is prepared to share. I apologise in advance that the code may not be the tidiest, or adhering to any conventions as I’m pretty new to Python. ;o)
Al.
from obspy.clients.fdsn import Client
from obspy.core import UTCDateTime, Stream
from obspy.signal import filter
import matplotlib.pyplot as plt
from obspy.taup import TauPyModel
import math
rs = Client('RASPISHAKE')
# set the station name and download the response information
stn = 'R21C0' # your station name
inv = rs.get_stations(network='AM', station=stn, level='RESP')
latS = -33.7691 #station latitude
lonS = 149.8843 #station longitude
eleS = 1114 #station elevation
#enter event data
eventTime = UTCDateTime(2022, 6, 30, 18, 40, 37) # (YYYY, m, d, H, M, S)
latE = 19.0 # quake latitude
lonE = 121.3 # quake longitude
depth = 35 # quake depth, km
mag = 6.4 # quake magnitude
eventID = 'rs2022mtvway' # ID for the event
locE = 'Philippine Islands Region' # location name
#set start time
delay = 540 # delay the start of the plot from the event
start = eventTime + delay # calculate the plot start time from the event and delay
end = start + 300 # calculate the end time from the start and duration
# set the FDSN server location and channel names
ch = 'EHZ' # ENx = accelerometer channels; EHx or SHZ = geophone channels
# get waveform and copy it for independent removal of instrument response
trace0 = rs.get_waveforms('AM', stn, '00', ch, start, end)
trace1 = trace0.copy()
trace2 = trace0.copy()
# calculate great circle angle of separation
# convert angles to radians
latSrad = math.radians(latS)
lonSrad = math.radians(lonS)
latErad = math.radians(latE)
lonErad = math.radians(lonE)
if lonSrad > lonErad:
lon_diff = lonSrad - lonErad
else:
lon_diff = lonErad - lonSrad
great_angle_rad = math.acos(math.sin(latErad)*math.sin(latSrad)+math.cos(latErad)*math.cos(latSrad)*math.cos(lon_diff))
great_angle_deg = math.degrees(great_angle_rad)
model = TauPyModel(model='iasp91')
arrs = model.get_travel_times(depth, great_angle_deg)
print(arrs)
#arrs.plot_times(plot_all=True)
arrivals = model.get_ray_paths(depth, great_angle_deg, phase_list=['ttbasic'])
#arrivals.plot_rays(plot_type='spherical', phase_list=['ttbasic'], legend=True)
# bandpass filter
filt = [0.1, 0.3, 10, 12]
# Create output traces
outdisp = trace0.remove_response(inventory=inv,pre_filt=filt,output='DISP',water_level=60, plot=False) # convert to Disp
outvel = trace1.remove_response(inventory=inv,pre_filt=filt,output='VEL',water_level=60, plot=False) # convert to Vel
outacc = trace2.remove_response(inventory=inv,pre_filt=filt,output='ACC',water_level=60, plot=False) # convert to Acc
# set up plot
fig = plt.figure(figsize=(20,10)) # set to page size in inches
ax1 = fig.add_subplot(5,2,1) # displacement waveform
ax2 = fig.add_subplot(5,2,3) # velocity Waveform
ax3 = fig.add_subplot(5,2,5) # acceleration waveform
ax4 = fig.add_subplot(5,2,7) # velocity spectrogram
ax5 = fig.add_subplot(5,2,9) # velocity PSD
ax7 = fig.add_subplot(5,2,(2,10)) # TAUp plot
fig.suptitle("AM."+stn+" M"+str(mag)+" "+locE+" "+eventID+eventTime.strftime(' %d/%m/%Y %H:%M:%S.%f UTC')) #Title of the figure
fig.text(0.05, 0.92, "Filter: "+str(filt[1])+" to "+str(filt[2])+"Hz")
#plot traces
ax1.plot(trace0[0].times(reftime=eventTime), outdisp[0].data, lw=1, color='b') # displacement waveform
#ax1.plot(arrs.times())
ax2.plot(trace0[0].times(reftime=eventTime), outvel[0].data, lw=1, color='g') # velocity Waveform
ax3.plot(trace0[0].times(reftime=eventTime), outacc[0].data, lw=1, color='r') # acceleration waveform
ax4.specgram(x=trace1[0], NFFT=64, noverlap=16, Fs=100, cmap='viridis') # velocity spectrogram
ax4.set_ylim(filt[0],filt[3])
ax5.psd(x=trace1[0], NFFT=512, noverlap=2, Fs=100, color='g', lw=1) # velocity PSD
ax5.set_xscale('log')
#ax7.arrivals.plot_rays(plot_type='spherical', phase_list=['ttbasic'], legend=True)
# set up some plot details
ax1.set_ylabel("Disp, m")
ax2.set_ylabel("Vel, m/s")
ax3.set_ylabel("Acc, m/s/s")
ax4.set_ylabel("V Freq, Hz")
ax5.set_ylabel("V PSD")
# get the limits of the y axis so text can be consistently placed
ax4b, ax4t = ax4.get_ylim()
ax4.text(2, ax4t*0.8, 'Plot Start Time: '+start.strftime(' %d/%m/%Y %H:%M:%S.%f UTC'))
#adjust subplots for readability
plt.subplots_adjust(hspace=0.3, wspace=0.1, left=0.05, right=0.95, bottom=0.05)
# save the final figure
#plt.savefig(str(mag)+'Quake'+eventID+eventTime.strftime('%Y%m%d %H%M%S UTC')) #comment this line out till figure is final
# show the final figure
plt.show()
