Generalize plotting script for sigma tau.

array_processing/tools/plotting.py

@@ -3,64 +3,114 @@
 import numpy as np
 
 
-def array_plot(tvec, data, t, mdccm, vel, az, mcthresh):
+def array_plot(st, t, mdccm, vel, baz, ccmplot=False,
+ mcthresh=None, sigma_tau=None):
  """
- Module to run plot array processing results
- @ authors: David Fee
-
- example: array_plot(stf,tvec,t,mdccm,vel,az,mcthresh):
+ Creates plots for velocity - back-azimuth array processing.
+
+ Args:
+ st: Filtered obspy stream. Assumes response has been removed.
+ t: Array processing time vector.
+ mdccm: Array of median cross-correlation maxima.
+ vel: Array of trace velocity estimates.
+ baz: Array of back-azimuth estimates.
+ ccmplot: Boolean flag to plot the Mean/Median cross-correlation
+ maxima values on the y-axis in addition to the color scale.
+ mcthresh: Add a dashed line at the [float] level
+ in the ccmplot subplot.
+ sigma_tau: Array of sigma_tau values.
+ The flag to add the sigma_tau subplot.
+
+ Returns:
+ fig1: Output figure.
+ axs1: Output figure axes.
+
+ Usage:
+ fig1, axs1= array_plot(st, t, mdccm, vel, baz, ccmplot=False,
+ mcthresh=None, sigma_tau=None)
  """
 
- cm='RdYlBu_r' #colormap
- cax=0.2,1 #colorbar/y-axis for mccm
-
- fig1, axarr1 = plt.subplots(4, 1,sharex='col')
- fig1.set_size_inches(10,9)
- axs1=axarr1.ravel()
- axs1[0].plot(tvec,data[:,0],'k')
+ # Specify the colormap.
+ cm = 'RdYlBu_r'
+ # Colorbar/y-axis limits for MdCCM.
+ cax = (0.2, 1)
+ # Specify the time vector for plotting the trace.
+ tvec = st[0].times('matplotlib')
+
+ # Determine the number and order of subplots.
+ num_subplots = 3
+ vplot = 1
+ bplot = 2
+ splot = bplot
+ if ccmplot:
+ num_subplots += 1
+ vplot += 1
+ bplot += 1
+ splot = bplot
+ if sigma_tau is not None:
+ num_subplots += 1
+ splot = bplot + 1
+
+ # Start Plotting.
+ # Plot the trace.
+ fig1, axs1 = plt.subplots(num_subplots, 1, sharex='col')
+ fig1.set_size_inches(10, 9)
+ axs1[0].plot(tvec, st[0].data, 'k')
  axs1[0].axis('tight')
  axs1[0].set_ylabel('Pressure [Pa]')
- #set_xlim(st[0].stats.starttime,st[0].stats.endtime)
-
- sc=axs1[1].scatter(t,mdccm,c=mdccm,edgecolors='k',lw=.3,cmap=cm)
- axs1[1].plot([t[0],t[-1]],[mcthresh,mcthresh],'r--')
- axs1[1].axis('tight')
- axs1[1].set_xlim(t[0],t[-1])
- axs1[1].set_ylim(cax)
- sc.set_clim(cax)
- axs1[1].set_ylabel('MdCCM')
 
- sc=axs1[2].scatter(t,vel,c=mdccm,edgecolors='k',lw=.3,cmap=cm)
- axs1[2].set_ylim(.25,.45)
- axs1[2].set_xlim(t[0],t[-1])
+ # Plot MdCCM on its own plot.
+ if ccmplot:
+ sc = axs1[1].scatter(t, mdccm, c=mdccm,
+ edgecolors='k', lw=0.3, cmap=cm)
+ if mcthresh is not None:
+ axs1[1].plot([t[0], t[-1]], [mcthresh, mcthresh], 'k--')
+ axs1[1].axis('tight')
+ axs1[1].set_xlim(t[0], t[-1])
+ axs1[1].set_ylim(cax)
+ sc.set_clim(cax)
+ axs1[1].set_ylabel('MdCCM')
+
+ # Plot the trace/apparent velocity.
+ sc = axs1[vplot].scatter(t, vel, c=mdccm, edgecolors='k', lw=0.3, cmap=cm)
+ axs1[vplot].set_ylim(0.25, 0.45)
+ axs1[vplot].set_xlim(t[0], t[-1])
  sc.set_clim(cax)
- axs1[2].set_ylabel('Trace Velocity\n [km/s]')
+ axs1[vplot].set_ylabel('Trace Velocity\n [km/s]')
 
- sc=axs1[3].scatter(t,az,c=mdccm,edgecolors='k',lw=.3,cmap=cm)
- #axs1[3].plot([t[0],t[-1]],[azvolc,azvolc],'r--')
- axs1[3].set_ylim(0,360)
-
- axs1[3].set_xlim(t[0],t[-1])
+ # Plot the back-azimuth.
+ sc = axs1[bplot].scatter(t, baz, c=mdccm, edgecolors='k', lw=0.3, cmap=cm)
+ axs1[bplot].set_ylim(0, 360)
+ axs1[bplot].set_xlim(t[0], t[-1])
  sc.set_clim(cax)
- axs1[3].set_ylabel('Back-azimuth\n [deg]')
-
- axs1[3].xaxis_date()
- axs1[3].tick_params(axis='x',labelbottom='on')
- axs1[3].fmt_xdata = dates.DateFormatter('%HH:%MM')
- axs1[3].xaxis.set_major_formatter(dates.DateFormatter("%d-%H:%M"))
- axs1[3].set_xlabel('UTC Time')
-
- cbot=.1
- ctop=axs1[1].get_position().y1
- cbaxes=fig1.add_axes([.92,cbot,.02,ctop-cbot])
- hc=plt.colorbar(sc,cax=cbaxes)
+ axs1[bplot].set_ylabel('Back-azimuth\n [deg]')
+
+ # Plot sigma_tau if given.
+ if sigma_tau is not None:
+ sc = axs1[splot].scatter(t, sigma_tau, c=mdccm,
+ edgecolors='k', lw=0.3, cmap=cm)
+ axs1[splot].set_xlim(t[0], t[-1])
+ sc.set_clim(cax)
+ axs1[splot].set_ylabel(r'$\sigma_\tau$')
+
+ axs1[splot].xaxis_date()
+ axs1[splot].tick_params(axis='x', labelbottom='on')
+ axs1[splot].fmt_xdata = dates.DateFormatter('%HH:%MM')
+ axs1[splot].xaxis.set_major_formatter(dates.DateFormatter("%d-%H:%M"))
+ axs1[splot].set_xlabel('UTC Time')
+
+ # Add the MdCCM colorbar.
+ cbot = axs1[splot].get_position().y0
+ ctop = axs1[1].get_position().y1
+ cbaxes = fig1.add_axes([0.92, cbot, 0.02, ctop-cbot])
+ hc = plt.colorbar(sc, cax=cbaxes)
  hc.set_label('MdCCM')
 
- return fig1,axs1
+ return fig1,axs1
 
 
 def arraySigPlt(rij, sig, sigV, sigTh, impResp, vel, th, kvec, figName=None):
- """
+r"""
  Plots output of arraySig method
 
  Parameters
@@ -87,11 +137,11 @@ def arraySigPlt(rij, sig, sigV, sigTh, impResp, vel, th, kvec, figName=None):
  Name of output file, will be written as figName.png (optional)
  """
 
- # for plotting methods & scripts
+ # For plotting methods & scripts
  figFormat = 'png' # MUCH faster than pdf!!
  figDpi = 600 # good resolution
 
- # lower RHS is array geometry
+ # Lower RHS is array geometry
  fig = plt.figure()
  axRij = plt.subplot(2, 2, 4)
  for h in range(rij.shape[1]):
@@ -101,14 +151,14 @@ def arraySigPlt(rij, sig, sigV, sigTh, impResp, vel, th, kvec, figName=None):
  axRij.axis('square')
  axRij.grid()
 
- # upper RHS is impulse reponse
+ # Upper RHS is impulse reponse
  axImp = plt.subplot(2, 2, 2)
  plt.pcolormesh(kvec, kvec, impResp)
  plt.ylabel('k$_y$ (km$^{-1}$)')
  plt.xlabel('k$_x$ (km$^{-1}$)')
  axImp.axis('square')
 
- # upper RHS is th uncertainty
+ # Upper RHS is theta uncertainty
  plt.subplot(2, 2, 1)
  meshTh = plt.pcolormesh(th, vel, sigTh)
  plt.ylabel('vel. (km/s)')
@@ -117,22 +167,22 @@ def arraySigPlt(rij, sig, sigV, sigTh, impResp, vel, th, kvec, figName=None):
  sigStr = str(sig)
  cbrTh.set_label(r'$\delta\theta\;\;\sigma_\tau = $' + sigStr + ' s')
 
- # lower RHS is vel uncertainty
+ # Lower RHS is velocity uncertainty
  plt.subplot(2, 2, 3)
  meshV = plt.pcolormesh(th, vel, sigV)
  plt.ylabel('vel. (km/s)')
  plt.xlabel(r'$\theta (\circ)$')
  cbrV = plt.colorbar(meshV, )
- cbrV.set_label('$\delta v$')
+ cbrV.set_label(r'$\delta v$')
 
- # prepare output & display in iPython workspace
- plt.tight_layout() # IGNORE renderer warning from script! it does just fine
+ # Prepare output & display in iPython workspace
+ plt.tight_layout() # IGNORE renderer warning from script! It is fine.
  if figName is not None:
  plt.savefig(figName + '.' + figFormat, format=figFormat, dpi=figDpi)
 
 
 def arraySigContourPlt(sigV, sigTh, vel, th, trace_v):
- """
+r"""
  Plots output of arraySig method onto a polar plot for a specified trace
  velocity.
 
@@ -160,12 +210,12 @@ def arraySigContourPlt(sigV, sigTh, vel, th, trace_v):
  """
 
  tvel_ptr = np.abs(vel - trace_v).argmin()
- sigV_cont = sigV[tvel_ptr,:]
- sigTh_cont = sigTh[tvel_ptr,:]
+ sigV_cont = sigV[tvel_ptr,:]
+ sigTh_cont = sigTh[tvel_ptr,:]
  theta = np.linspace(0, 2 * np.pi, len(sigV_cont))
 
-
-fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, subplot_kw={'projection': 'polar'})
+ fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2,
+ subplot_kw={'projection': 'polar'})
 
  ax1.set_theta_direction(-1)
  ax1.set_theta_offset(np.pi/2.0)
@@ -174,7 +224,8 @@ def arraySigContourPlt(sigV, sigTh, vel, th, trace_v):
  ax1.yaxis.get_major_locator().base.set_params(nbins=6)
  ax1.set_rlabel_position(22.5)
  ax1.grid(True)
- ax1.set_title('Trace Velocity Uncertainty, V=%.2f' % trace_v, va='bottom', pad=20)
+ ax1.set_title('Trace Velocity Uncertainty, V=%.2f' % trace_v,
+ va='bottom', pad=20)
 
  ax2.set_theta_direction(-1)
  ax2.set_theta_offset(np.pi/2.0)
@@ -183,6 +234,7 @@ def arraySigContourPlt(sigV, sigTh, vel, th, trace_v):
  ax2.yaxis.get_major_locator().base.set_params(nbins=6)
  ax2.set_rlabel_position(22.5)
  ax2.grid(True)
- ax2.set_title('Back-Azimuth Uncertainty, V=%.2f' % trace_v, va='bottom', pad=20)
+ ax2.set_title('Back-Azimuth Uncertainty, V=%.2f' % trace_v,
+ va='bottom', pad=20)
 
  return fig

example.py

@@ -1,4 +1,4 @@
-#%% User-defined parameters
+#%% User-defined parameters
 
 from waveform_collection import gather_waveforms
 from obspy.core import UTCDateTime
@@ -26,7 +26,7 @@
 
 stf = st.copy()
 stf.filter('bandpass', freqmin=FMIN, freqmax=FMAX, corners=2, zerophase=True)
-stf.taper(max_percentage=.01)
+stf.taper(max_percentage=0.01)
 
 tvec = stf[0].times('matplotlib')
 
@@ -42,7 +42,7 @@
 
 vel, az, mdccm, t, data = wlsqva_proc(stf, rij, tvec, WINLEN, WINOVER)
 
-fig1, axs1 = array_plot(tvec, data, t, mdccm, vel, az, .6)
+fig1, axs1 = array_plot(stf, t, mdccm, vel, az, ccmplot=True, sigma_tau=False)
 
 #%% Array uncertainty