Adding a secondary axis to subplots with UTC time

Developers
1

I have a report with several subplots with the x axis in seconds after the event. I want to add a secondary axis, showing the UTC. I have seen plots that others have posted that does just this.

I have a successful secondary axis working on my PSD plot, so I suspect the problem lies in handling the date/time values and perhaps the format of UTCDateTime (though I have tried converting to datetime, etc). I’d appreciate any help anyone can offer about where I’m going wrong and/or what I’m missing.

Here is my code so far:


from obspy.clients.fdsn import Client
from obspy.core import UTCDateTime, Stream
from obspy.signal import filter
import matplotlib.pyplot as plt
import numpy as np
import datetime
import matplotlib.dates as mdates
from obspy.taup import TauPyModel
from matplotlib.ticker import AutoMinorLocator
import math
rs = Client('RASPISHAKE')

def plot_arrivals(ax):
    y1 = -1
    axb, axt = ax.get_ylim()               # calculate the y limits of the graph
    for q in range(0, no_arrs):            #plot each arrival in turn
        x1 = arrs[q].time              # extract the time to plot
        if (x1 >= delay):
            if x1 < delay+duration:
                ax.axvline(x=x1, linewidth=0.5, linestyle='--', color='black')      # draw a vertical line
                if y1 < 0:                      # alternate top and bottom for phase tags
                    y1 = axt*0.8
                else:
                    y1 = axb*0.95
                ax.text(x1,y1,arrs[q].name, alpha=0.5)                                 # print the phase name

def time2UTC(a):        #convert time (seconds) since event back to UTCDateTime
    tUTC = eventTime + a
    return tUTC.strftime('%H:%M:%S')

def uTC2time(a):        #convert UTCDateTime to seconds since the event
    return a - eventTime

def one_over(a):            # 1/x to convert frequency to period
    #Vectorized 1/a, treating a==0 manually
    a = np.array(a).astype(float)
    near_zero = np.isclose(a, 0)
    a[near_zero] = np.inf
    a[~near_zero] = 1 / a[~near_zero]
    return a

inverse = one_over          #function 1/x is its own inverse

# set the station name and download the response information
stn = 'R21C0'                            # your station name
inv = rs.get_stations(network='AM', station=stn, level='RESP')  # get the instrument response
latS = -33.7691     #station latitude
lonS = 149.8843     #station longitude
eleS = 1114         #station elevation

#enter event data
eventTime = UTCDateTime(2022, 7, 24, 19, 30, 51) # (YYYY, m, d, H, M, S)
latE = -14.9                            # quake latitude
lonE = 167.4                           # quake longitude
depth = 129                             # quake depth, km
mag = 5.0                              # quake magnitude
eventID = 'rs2022olyyms'               # ID for the event
locE = 'Vanuatu Islands'     # location name

#set start time
delay = 120                     # delay the start of the plot from the event
duration = 900                  # duration of plots
start = eventTime + delay       # calculate the plot start time from the event and delay
end = start + duration               # 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)
trace0.merge(method=0, fill_value='latest')         #fill in any gaps in the data to prevent a crash
trace0.detrend(type='demean')                       #demean the data
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)     #great circle angle between quake and station
distance = great_angle_rad*12742/2      #calculate distance between quake and station in km

model = TauPyModel(model='iasp91')
arrs = model.get_travel_times(depth, great_angle_deg)
print(arrs)
no_arrs = len(arrs)     # the number of arrivals

# bandpass filter
filt = [0.7, 0.7, 45, 50]      

# 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,12))       # 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,8), polar=True)       # 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")          # Filter details
fig.text(0.2, 0.92, 'Separation = '+str(round(great_angle_deg,3))+u"\N{DEGREE SIGN}"+' or '+str(int(distance))+'km.')   #distance between quake and station
fig.text(0.4, 0.92, 'Latitude: '+str(latE)+u"\N{DEGREE SIGN}"+' Longitude: '+str(lonE)+u"\N{DEGREE SIGN}"+' Depth: '+str(depth)+'km.')  #quake lat, lon and depth

#plot traces
ax1.plot(trace0[0].times(reftime=eventTime), outdisp[0].data, lw=1, color='b')      # displacement waveform
#ax1.set_ylim(-3e-8,3e-8)         # set manual y limits for displacement- comment this out for autoscaling
ax2.plot(trace0[0].times(reftime=eventTime), outvel[0].data, lw=1, color='g')       # velocity Waveform
#ax2.set_ylim(-1e-7,1e-7)         # set manual y limits for velocity - comment this out for autoscaling
ax3.plot(trace0[0].times(reftime=eventTime), outacc[0].data, lw=1, color='r')       # acceleration waveform
#ax3.set_ylim(-1e-6,1e-6)         # set manual y limits for acceleration - comment this out for auto scaling
ax4.specgram(x=trace1[0], NFFT=32, noverlap=2, Fs=100, cmap='viridis')         # velocity spectrogram
ax4.set_ylim(0,filt[3])
ax5.psd(x=trace1[0], NFFT=duration, noverlap=0, Fs=100, color='g', lw=1)             # velocity PSD
ax5.set_xscale('log')

#plot secondary axes
secax_x1 = ax1.secondary_xaxis('top', functions=(time2UTC, uTC2time))
secax_x1.set_xlabel('UTC')
secax_x5 = ax5.secondary_xaxis('top', functions=(one_over, inverse))
secax_x5.set_xlabel('Period, s', labelpad=0)

#test conversion functions for secondary time axes
print(time2UTC(120))
print(uTC2time(eventTime+360))
                               
# build array of arrival data
x2 = 0.49        #start near middle of page but maximise list space
dx = 0.0066      # linespacing
fig.text(x2, 0.03, 'Phase',size='small', rotation=90)         #print headings
fig.text(x2, 0.09, 'Time',size='small', rotation=90)
fig.text(x2, 0.15, 'UTC',size='small', rotation=90)
fig.text(x2, 0.2, 'Vertical Component', alpha = 0.5, size='small', rotation=90)
for i in range (0, no_arrs):                    #print data array
    x2 += dx
    fig.text(x2, 0.03, arrs[i].name, size='small', rotation=90)       #print phase name
    fig.text(x2, 0.09, str(round(arrs[i].time,3))+'s', size='small', rotation=90)     #print arrival time
    arrtime = eventTime + arrs[i].time
    fig.text(x2, 0.15, arrtime.strftime('%H:%M:%S'), size='small', rotation=90)
    if arrs[i].name.endswith('P') or arrs[i].name.endswith('p') or arrs[i].name.endswith('Pdiff') or arrs[i].name.endswith('Pn'):                    #calculate and print the vertical component of the signal
        fig.text(x2, 0.2, str(round(100*math.cos(math.radians(arrs[i].incident_angle)),1))+'%', alpha = 0.5, size='small', rotation=90)
    elif arrs[i].name.endswith('S') or arrs[i].name.endswith('s') or arrs[i].name.endswith('Sn'):
        fig.text(x2, 0.2, str(round(100*math.sin(math.radians(arrs[i].incident_angle)),1))+'%', alpha = 0.5, size='small', rotation=90)
x2 += 2*dx
fig.text(x2, 0.03, str(no_arrs)+' arrivals total.', size='small', rotation=90)     #print number of arrivals

# print phase key
x2 = 0.91
fig.text(x2, 0.03, 'Phase Key', size='small', rotation=90)
pkey = ['P:   compression wave', 'p:   strictly upward compression wave', 'S:   shear wave', 's:   strictly upward shear wave', 'K:   compression wave in outer core', 'I:   compression wave in inner core', 'c:   reflection off outer core', 'diff:   diffracted wave along core mantle boundary', 'i:   reflection off inner core', 'n:   wave follows the Moho (crust/mantle boundary)']
for i in range (0, 10):
    x2 +=dx
    fig.text(x2, 0.03, pkey[i], size='small', rotation=90)

#plot phase arrivals
plot_arrivals(ax1)          #plot arrivals on displacement plot
plot_arrivals(ax2)          #plot arrivals on velocity plot
plot_arrivals(ax3)          #plot arrivals on acceleration plot

# set up some plot details
ax1.set_ylabel("Vert. Disp, m")
ax2.set_ylabel("Vert. Vel, m/s")
ax3.set_ylabel("Vert. Acc, m/s/s")         
ax4.set_ylabel("Vel. Freq, Hz")
ax5.set_ylabel("Vel. PSD, dB")

# get the limits of the y axis so text can be consistently placed
ax4b, ax4t = ax4.get_ylim()
ax4.text(2, ax4t*0.85, 'Plot Start Time: '+start.strftime(' %d/%m/%Y %H:%M:%S.%f UTC (')+str(delay)+' seconds after event).', size='small')      # explain difference in x time scale

#adjust subplots for readability
plt.subplots_adjust(hspace=0.7, wspace=0.1, left=0.05, right=0.95, bottom=0.05)

#plot the ray paths
arrivals = model.get_ray_paths(depth, great_angle_deg, phase_list=['ttbasic'])      #calculate the ray paths
ax7 = arrivals.plot_rays(plot_type='spherical', ax=ax7, fig=fig, phase_list=['ttbasic'], show=False, legend=True)   #plot the ray paths

# Annotate regions
ax7.text(0, 0, 'Solid\ninner\ncore',
        horizontalalignment='center', verticalalignment='center',
        bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))
ocr = (model.model.radius_of_planet -
       (model.model.s_mod.v_mod.iocb_depth +
        model.model.s_mod.v_mod.cmb_depth) / 2)
ax7.text(math.radians(180), ocr, 'Fluid outer core',
        horizontalalignment='center',
        bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))
mr = model.model.radius_of_planet - model.model.s_mod.v_mod.cmb_depth / 2
ax7.text(math.radians(180), mr, 'Solid mantle',
        horizontalalignment='center',
        bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))

# save the final figure
#plt.savefig('E:\Pictures\Raspberry Shake and Boom\M'+str(mag)+'Quake '+locE+eventID+eventTime.strftime('%Y%m%d %H%M%S UTC.png'))  #comment this line out till figure is final

# show the final figure
plt.show()

Lines 131 and 132 are where I try to add the secondary axis to the first subplot ax1.

Lines 137 and 138 are currently to test that the functions for the time conversion work correctly.

FWIW here are the errors I am getting:

Exception in Tkinter callback
Traceback (most recent call last):
  File "F:\Program Files\Python310\lib\tkinter\__init__.py", line 1921, in __call__
    return self.func(*args)
  File "F:\Program Files\Python310\lib\tkinter\__init__.py", line 839, in callit
    func(*args)
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\backends\_backend_tk.py", line 251, in idle_draw
    self.draw()
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\backends\backend_tkagg.py", line 9, in draw
    super().draw()
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\backends\backend_agg.py", line 436, in draw
    self.figure.draw(self.renderer)
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\artist.py", line 73, in draw_wrapper
    result = draw(artist, renderer, *args, **kwargs)
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\artist.py", line 50, in draw_wrapper
    return draw(artist, renderer)
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\figure.py", line 2796, in draw
    artists = self._get_draw_artists(renderer)
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\figure.py", line 245, in _get_draw_artists
    child.apply_aspect(pos)
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\axes\_secondary_axes.py", line 124, in apply_aspect
    self._set_lims()
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\axes\_secondary_axes.py", line 226, in _set_lims
    lims = self._functions[0](np.array(lims))
  File "F:\Program Files\Python310\QReport6dev.py", line 28, in time2UTC
    tUTC = eventTime + a
  File "F:\Program Files\Python310\lib\site-packages\obspy\core\utcdatetime.py", line 996, in __add__
    return UTCDateTime(ns=self._ns + int(round(value * 1e9)))
TypeError: type numpy.ndarray doesn't define __round__ method
Exception in Tkinter callback
Traceback (most recent call last):
  File "F:\Program Files\Python310\lib\tkinter\__init__.py", line 1921, in __call__
    return self.func(*args)
  File "F:\Program Files\Python310\lib\tkinter\__init__.py", line 839, in callit
    func(*args)
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\backends\_backend_tk.py", line 251, in idle_draw
    self.draw()
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\backends\backend_tkagg.py", line 9, in draw
    super().draw()
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\backends\backend_agg.py", line 436, in draw
    self.figure.draw(self.renderer)
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\artist.py", line 73, in draw_wrapper
    result = draw(artist, renderer, *args, **kwargs)
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\artist.py", line 50, in draw_wrapper
    return draw(artist, renderer)
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\figure.py", line 2796, in draw
    artists = self._get_draw_artists(renderer)
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\figure.py", line 245, in _get_draw_artists
    child.apply_aspect(pos)
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\axes\_secondary_axes.py", line 124, in apply_aspect
    self._set_lims()
  File "F:\Program Files\Python310\lib\site-packages\matplotlib\axes\_secondary_axes.py", line 226, in _set_lims
    lims = self._functions[0](np.array(lims))
  File "F:\Program Files\Python310\QReport6dev.py", line 28, in time2UTC
    tUTC = eventTime + a
  File "F:\Program Files\Python310\lib\site-packages\obspy\core\utcdatetime.py", line 996, in __add__
    return UTCDateTime(ns=self._ns + int(round(value * 1e9)))
TypeError: type numpy.ndarray doesn't define __round__ method
1 Like
2

Hello sheeny,

From what I can see, the test lines work indeed well, but when the same function is passed for the plot, instead of a single time value, an array of numbers is fed in the function, causing the issue type numpy.ndarray doesn't define __round__ method that you can see.

One potential solution could be first acquiring all the datetime values that you want to convert/convert back in an array and then processing them one by one with a for cycle that uses the same function.

Another would instead be to acquire all the timestamp locations that you want to plot (likely, on the same x coordinate as the times in the lower axis) and then process them again one by one, but only those values and nothing else.

As said, the function works perfectly for single time values, but doesn’t know what to do when an array of multiple times is fed into it.

3

G’Day Stormchaser,

I’ve been playing around with this a bit (but also been away so not full time ;o) ) and I’m no further advanced - still getting the same errors.

I tried the array of values as suggested, but to no avail. I print them out before plotting the figure and they appear OK.

I thought also that the strftime format in the function may have been a problem, so I removed that back to basic UTCDateTime so the functions are truly the inverse of each other. That also made no difference.

I also noticed the 1e9 in the error is for some reason multiplying the nanoseconds (ns) part of the UTCDateTime, so perhaps I need to use another date/time format?

I feel I may not entirely understand how you intended I should handle the arrays, so would appreciate it if you could cast any eye over where I’ve got to now:


from obspy.clients.fdsn import Client
from obspy.core import UTCDateTime, Stream
from obspy.signal import filter
import matplotlib.pyplot as plt
import numpy as np
import datetime
import matplotlib.dates as mdates
from obspy.taup import TauPyModel
from matplotlib.ticker import AutoMinorLocator
import math
rs = Client('RASPISHAKE')

def plot_arrivals(ax):
    y1 = -1
    axb, axt = ax.get_ylim()               # calculate the y limits of the graph
    for q in range(0, no_arrs):            #plot each arrival in turn
        x1 = arrs[q].time              # extract the time to plot
        if (x1 >= delay):
            if x1 < delay+duration:
                ax.axvline(x=x1, linewidth=0.5, linestyle='--', color='black')      # draw a vertical line
                if y1 < 0:                      # alternate top and bottom for phase tags
                    y1 = axt*0.8
                else:
                    y1 = axb*0.95
                ax.text(x1,y1,arrs[q].name, alpha=0.5)                                 # print the phase name

def time2UTC(a):        #convert time (seconds) since event back to UTCDateTime
    return eventTime + a

def uTC2time(a):        #convert UTCDateTime to seconds since the event
    return a - eventTime

def one_over(a):            # 1/x to convert frequency to period
    #Vectorized 1/a, treating a==0 manually
    a = np.array(a).astype(float)
    near_zero = np.isclose(a, 0)
    a[near_zero] = np.inf
    a[~near_zero] = 1 / a[~near_zero]
    return a

inverse = one_over          #function 1/x is its own inverse

# set the station name and download the response information
stn = 'R21C0'                            # your station name
inv = rs.get_stations(network='AM', station=stn, level='RESP')  # get the instrument response
latS = -33.7691     #station latitude
lonS = 149.8843     #station longitude
eleS = 1114         #station elevation

#enter event data
eventTime = UTCDateTime(2022, 7, 24, 19, 30, 51) # (YYYY, m, d, H, M, S)
latE = -14.9                            # quake latitude
lonE = 167.4                           # quake longitude
depth = 129                             # quake depth, km
mag = 5.0                              # quake magnitude
eventID = 'rs2022olyyms'               # ID for the event
locE = 'Vanuatu Islands'     # location name

#set start time
delay = 120                     # delay the start of the plot from the event
duration = 900                  # duration of plots
start = eventTime + delay       # calculate the plot start time from the event and delay
end = start + duration               # 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)
trace0.merge(method=0, fill_value='latest')         #fill in any gaps in the data to prevent a crash
trace0.detrend(type='demean')                       #demean the data
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)     #great circle angle between quake and station
distance = great_angle_rad*12742/2      #calculate distance between quake and station in km

model = TauPyModel(model='iasp91')
arrs = model.get_travel_times(depth, great_angle_deg)
print(arrs)
no_arrs = len(arrs)     # the number of arrivals

# bandpass filter
filt = [0.7, 0.7, 45, 50]      

# 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,12))       # 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,8), polar=True)       # 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")          # Filter details
fig.text(0.2, 0.92, 'Separation = '+str(round(great_angle_deg,3))+u"\N{DEGREE SIGN}"+' or '+str(int(distance))+'km.')   #distance between quake and station
fig.text(0.4, 0.92, 'Latitude: '+str(latE)+u"\N{DEGREE SIGN}"+' Longitude: '+str(lonE)+u"\N{DEGREE SIGN}"+' Depth: '+str(depth)+'km.')  #quake lat, lon and depth

#plot traces
ax1.plot(trace0[0].times(reftime=eventTime), outdisp[0].data, lw=1, color='b')      # displacement waveform
#ax1.set_ylim(-3e-8,3e-8)         # set manual y limits for displacement- comment this out for autoscaling
ax2.plot(trace0[0].times(reftime=eventTime), outvel[0].data, lw=1, color='g')       # velocity Waveform
#ax2.set_ylim(-1e-7,1e-7)         # set manual y limits for velocity - comment this out for autoscaling
ax3.plot(trace0[0].times(reftime=eventTime), outacc[0].data, lw=1, color='r')       # acceleration waveform
#ax3.set_ylim(-1e-6,1e-6)         # set manual y limits for acceleration - comment this out for auto scaling
ax4.specgram(x=trace1[0], NFFT=32, noverlap=2, Fs=100, cmap='viridis')         # velocity spectrogram
ax4.set_ylim(0,filt[3])
ax5.psd(x=trace1[0], NFFT=duration, noverlap=0, Fs=100, color='g', lw=1)             # velocity PSD
ax5.set_xscale('log')

#plot secondary axes
if duration <= 90:
    dt=10
elif duration <= 200:
    dt=20
elif duration <= 300:
    dt=30
elif duration <= 600:
    dt=60
elif duration <= 1200:
    dt=120
elif duration <= 2400:
    dt=240
else:
    dt=300
tbase = start - start.second +(int(start.second/dt)+1)*dt
tlabels = []
tticks = []
nticks = int(duration/dt)-1
for k in range (0, nticks):
    tlabels.append((tbase+k*dt).strftime('%H:%M:%S'))
    tticks.append(uTC2time(tbase+k*dt))
print(tlabels)
print(tbase.strftime('%H:%M:%S')) 
print(tticks)
#secax_x1 = ax1.secondary_xaxis('top', functions=(time2UTC, uTC2time))
#secax_x1.set_xticks(ticks=tticks)
#secax_x1.set_xlabel('UTC')
secax_x2 = ax2.secondary_xaxis('top', functions=(time2UTC, uTC2time))
secax_x2.set_xticks(ticks=tticks)
secax_x2.set_xlabel('UTC')
secax_x5 = ax5.secondary_xaxis('top', functions=(one_over, inverse))
secax_x5.set_xlabel('Period, s', labelpad=0)

#test conversion functions for secondary time axes
print(time2UTC(120))
print(uTC2time(eventTime+360))
                               
# build array of arrival data
x2 = 0.49        #start near middle of page but maximise list space
dx = 0.0066      # linespacing
fig.text(x2, 0.03, 'Phase',size='small', rotation=90)         #print headings
fig.text(x2, 0.09, 'Time',size='small', rotation=90)
fig.text(x2, 0.15, 'UTC',size='small', rotation=90)
fig.text(x2, 0.2, 'Vertical Component', alpha = 0.5, size='small', rotation=90)
pphases=[]          #create an array of phases for the plot
allphases=True     # True if all phases to be plotted otherwise only those in the plotted time window are plotted
for i in range (0, no_arrs):                    #print data array
    x2 += dx
    fig.text(x2, 0.03, arrs[i].name, size='small', rotation=90)       #print phase name
    fig.text(x2, 0.09, str(round(arrs[i].time,3))+'s', size='small', rotation=90)     #print arrival time
    arrtime = eventTime + arrs[i].time
    fig.text(x2, 0.15, arrtime.strftime('%H:%M:%S'), size='small', rotation=90)
    if allphases or (arrs[i].time >= delay and arrs[i].time < delay+duration):      # build the list of phases to plot
        pphases.append(arrs[i].name)
    if arrs[i].name.endswith('P') or arrs[i].name.endswith('p') or arrs[i].name.endswith('Pdiff') or arrs[i].name.endswith('Pn'):                    #calculate and print the vertical component of the signal
        fig.text(x2, 0.2, str(round(100*math.cos(math.radians(arrs[i].incident_angle)),1))+'%', alpha = 0.5, size='small', rotation=90)
    elif arrs[i].name.endswith('S') or arrs[i].name.endswith('s') or arrs[i].name.endswith('Sn'):
        fig.text(x2, 0.2, str(round(100*math.sin(math.radians(arrs[i].incident_angle)),1))+'%', alpha = 0.5, size='small', rotation=90)
x2 += 2*dx
fig.text(x2, 0.03, str(no_arrs)+' arrivals total.', size='small', rotation=90)     #print number of arrivals

# print phase key
x2 = 0.91
fig.text(x2, 0.03, 'Phase Key', size='small', rotation=90)
pkey = ['P:   compression wave', 'p:   strictly upward compression wave', 'S:   shear wave', 's:   strictly upward shear wave', 'K:   compression wave in outer core', 'I:   compression wave in inner core', 'c:   reflection off outer core', 'diff:   diffracted wave along core mantle boundary', 'i:   reflection off inner core', 'n:   wave follows the Moho (crust/mantle boundary)']
for i in range (0, 10):
    x2 +=dx
    fig.text(x2, 0.03, pkey[i], size='small', rotation=90)

#plot phase arrivals
plot_arrivals(ax1)          #plot arrivals on displacement plot
plot_arrivals(ax2)          #plot arrivals on velocity plot
plot_arrivals(ax3)          #plot arrivals on acceleration plot

# set up some plot details
ax1.set_ylabel("Vert. Disp, m")
ax2.set_ylabel("Vert. Vel, m/s")
ax3.set_ylabel("Vert. Acc, m/s/s")         
ax4.set_ylabel("Vel. Freq, Hz")
ax5.set_ylabel("Vel. PSD, dB")

# get the limits of the y axis so text can be consistently placed
ax4b, ax4t = ax4.get_ylim()
ax4.text(2, ax4t*0.85, 'Plot Start Time: '+start.strftime(' %d/%m/%Y %H:%M:%S.%f UTC (')+str(delay)+' seconds after event).', size='small')      # explain difference in x time scale

#adjust subplots for readability
plt.subplots_adjust(hspace=0.7, wspace=0.1, left=0.05, right=0.95, bottom=0.05)

#plot the ray paths
arrivals = model.get_ray_paths(depth, great_angle_deg, phase_list=pphases)      #calculate the ray paths
ax7 = arrivals.plot_rays(plot_type='spherical', ax=ax7, fig=fig, phase_list=pphases, show=False, legend=True)   #plot the ray paths

# Annotate regions
ax7.text(0, 0, 'Solid\ninner\ncore',
        horizontalalignment='center', verticalalignment='center',
        bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))
ocr = (model.model.radius_of_planet -
       (model.model.s_mod.v_mod.iocb_depth +
        model.model.s_mod.v_mod.cmb_depth) / 2)
ax7.text(math.radians(180), ocr, 'Fluid outer core',
        horizontalalignment='center',
        bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))
mr = model.model.radius_of_planet - model.model.s_mod.v_mod.cmb_depth / 2
ax7.text(math.radians(180), mr, 'Solid mantle',
        horizontalalignment='center',
        bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))

# save the final figure
#plt.savefig('E:\Pictures\Raspberry Shake and Boom\M'+str(mag)+'Quake '+locE+eventID+eventTime.strftime('%Y%m%d %H%M%S UTC.png'))  #comment this line out till figure is final

# show the final figure
plt.show()
2 Likes
4

Hello sheeny, sorry for the delayed response,

Thank you for the new code, I think that, at first glance, it should work. However, since it is not, there must be something minor that is causing it to fail, as it usually goes.

I will review it and get back as soon as I find a possible solution that you can use.

5

Thanks for that Stormchaser. I do appreciate it, but if anyone has done this and is prepared to share some code, that may be a simpler way to see what I’m doing wrong.

Thanks,

Al.

1 Like
6

Hello sheeny, it’s no problem at all.

Thanks for the additional code as usual, I think you were almost there, especially with what you have written here in this section:

tbase = start - start.second +(int(start.second/dt)+1)*dt
tlabels = []
tticks = []
nticks = int(duration/dt)-1
for k in range (0, nticks):
    tlabels.append((tbase+k*dt).strftime('%H:%M:%S'))
    tticks.append(uTC2time(tbase+k*dt))
print(tlabels)
print(tbase.strftime('%H:%M:%S')) 
print(tticks)

Just after it, you start with the secax, where I added this line:

secax_x1.set_xticklabels(tlabels)

to this portion of the code:

secax_x1 = ax1.secondary_xaxis('top')
secax_x1.set_xticks(ticks=tticks)
secax_x1.set_xticklabels(tlabels)
secax_x1.set_xlabel('UTC')

Since you are already producing an array of times, there is no need to recall again the functions time2UTC and UTC2time on the main line of the secax. And this is the result, with no errors:

image
image1276×765 208 KB

From there, it will just be a matter of adapting the code to all plots and maybe display it in the way you have in mind. And yes, if someone else has some better code to provide, feel free to share it with all the community!

7

Good on you, Stormchaser! Thank you!

I actually tried that, but it was still failing because I didn’t realise I need to remove the function calls from the secondary axis command.

I had a feeling it would be something relatively minor - I just wasn’t seeing it.

Thanks again!

Al.

2 Likes
8

Happy to have been of help!

Yes, I think that the fact that both the time arrays and the functions call on the secax was the root of the problem. Eager to see the new plots on Twitter!