cleaned up tilt figure; saving tilt to file

Author: paudetseis

obstools/atacr/plotting.py

@@ -27,6 +27,7 @@
 
 import numpy as np
 from matplotlib import pyplot as plt
+from matplotlib.gridspec import GridSpec
 from obstools.atacr import utils
 from obspy import Trace
 
@@ -237,7 +238,7 @@ def fig_av_cross(f, field, gooddays, ftype, ncomp, key='',
     return plt
 
 
-def fig_coh_ph(coh, ph, direc):
+def fig_coh_ph(coh, ph, direc, tilt, date):
     """
     Function to plot the coherence and phase between the rotated H and Z
     components, used to characterize the tilt direction.
@@ -250,21 +251,46 @@ def fig_coh_ph(coh, ph, direc):
         Phase between rotated H and Z components
     direc : :mod:`~numpy.ndarray`
         Directions considered in maximizing coherence between H and Z
+    date : list
+        List of datetime dates for plotting as function of time
+    tilt : list
+        List of tilt azimuths determined for each date
 
     """
 
     colors = plt.cm.cividis(np.linspace(0, 1, coh.shape[0]))
 
     if coh.ndim > 1:
-        f, (ax1, ax2) = plt.subplots(1, 2)
-        for i, (co, p) in enumerate(zip(coh, ph)):
+
+        meantilt = np.mean(tilt)
+        stdetilt = np.std(tilt, ddof=1)/np.sqrt(coh.shape[0])
+        mediantilt = np.median(tilt)
+
+        plt.rcParams['date.converter'] = 'concise'
+        f = plt.figure()
+        gs = GridSpec(3, 2, figure=f)
+        ax1 = f.add_subplot(gs[0:-1, 0])
+        ax2 = f.add_subplot(gs[0:-1, 1])
+        ax3 = f.add_subplot(gs[2, :])
+        ax3.axhline(
+            np.mean(tilt),
+            ls='--',
+            label='Mean: {0:.0f} $\pm$ {1:.1f}'.format(meantilt, stdetilt))
+        ax3.axhline(
+            np.median(tilt),
+            ls=':',
+            label='Median: {0:.0f}'.format(np.median(tilt)))
+        for i, (co, p, d, t) in enumerate(zip(coh, ph, date, tilt)):
             ax1.plot(direc, co, c=colors[i])
             ax2.plot(direc, p*180./np.pi, c=colors[i])
-        ax1.set_ylabel('Coherence')
+            ax3.plot(d, t, 'o', c=colors[i])
+        ax3.legend(loc='best', fontsize=8)
+        ax1.set_ylabel('Coherence Z-H1')
         ax1.set_ylim((0, 1.))
-        ax2.set_ylabel('Phase')
-        ax1.set_xlabel('Angle from H1')
-        ax2.set_xlabel('Angle from H1')
+        ax2.set_ylabel('Phase Z-H1')
+        ax3.set_ylabel('Tilt direction (deg)')
+        ax1.set_xlabel('Angle clockwise from H1 (deg)')
+        ax2.set_xlabel('Angle clockwise from H1 (deg)')
         plt.tight_layout()
 
     else:

obstools/atacr/utils.py

@@ -408,9 +408,8 @@ def calculate_tilt(ft1, ft2, ftZ, ftP, f, goodwins, tiltfreqs):
     phi = np.arange(0., 360., 10.)
     coh = np.zeros(len(phi))
     ph = np.zeros(len(phi))
-    cZZ = np.abs(
-        np.mean(ftZ[goodwins, :] * \
-            np.conj(ftZ[goodwins, :]), axis=0))[0:len(f)]
+    cZZ = np.abs(np.mean(ftZ[goodwins, :] *
+                         np.conj(ftZ[goodwins, :]), axis=0))[0:len(f)]
 
     for i, d in enumerate(phi):
 
@@ -486,6 +485,10 @@ def calculate_tilt(ft1, ft2, ftZ, ftP, f, goodwins, tiltfreqs):
     phase_value = rph[ind[0]][0]
     coh_value = rcoh[ind[0]][0]
     tilt = rphi[ind[0]][0]
+    if tilt > 360.:
+        tilt -= 360.
+    elif tilt < 0:
+        tilt += 360.
 
     # Now calculate spectra at tilt direction
     ftH = rotate_dir(ft1, ft2, tilt)

obstools/scripts/atacr_clean_spectra.py

@@ -162,11 +162,11 @@ def get_cleanspec_arguments(argv=None):
         help="Plot daily average figure. " +
         "[Default does not plot figure]")
     FigureGroup.add_argument(
-        "--figCoh",
+        "--figTilt",
         action="store_true",
-        dest="fig_coh_ph",
+        dest="fig_tilt",
         default=False,
-        help="Plot Coherence and Phase figure. " +
+        help="Plot Coherence, phase and tilt direction figure. " +
         "[Default does not plot figure]")
     FigureGroup.add_argument(
         "--figCross",
@@ -349,6 +349,7 @@ def main(args=None):
         dstart = str(tstart.year).zfill(4)+'.'+str(tstart.julday).zfill(3)+'-'
         dend = str(tend.year).zfill(4)+'.'+str(tend.julday).zfill(3)+'.'
         fileavst = avstpath / (dstart+dend+'avg_sta.pkl')
+        filetilt = avstpath / (dstart+dend+'tilt.txt')
 
         if fileavst.exists():
             if not args.ovr:
@@ -383,14 +384,18 @@ def main(args=None):
         # Initialize StaNoise object
         stanoise = StaNoise()
 
+        # Date + tilt list
+        date_list = []
+        tilt_list = []
+
         # Loop through each day withing time range
         while t1 < tend:
 
             year = str(t1.year).zfill(4)
             jday = str(t1.julday).zfill(3)
 
-            tstamp = year+'.'+jday+'.'
-            filespec = specpath / (tstamp + 'spectra.pkl')
+            tstamp = year+'.'+jday
+            filespec = specpath / (tstamp + '.spectra.pkl')
 
             # Load file if it exists
             if filespec.exists():
@@ -401,6 +406,8 @@ def main(args=None):
                 file = open(filespec, 'rb')
                 daynoise = pickle.load(file)
                 file.close()
+                tilt_list.append(daynoise.rotation.tilt)
+                date_list.append(t1.date)
                 stanoise += daynoise
             else:
                 t1 += 3600.*24.
@@ -526,21 +533,31 @@ def main(args=None):
                    ad_2Z_all, ad_2P_all, ad_ZP_all)
 
         # Quality control to identify outliers
-        stanoise.QC_sta_spectra(pd=args.pd, tol=args.tol, alpha=args.alpha,
-                                fig_QC=args.fig_QC, debug=args.debug,
-                                save=plotpath, form=args.form)
+        stanoise.QC_sta_spectra(
+            pd=args.pd,
+            tol=args.tol,
+            alpha=args.alpha,
+            fig_QC=args.fig_QC,
+            debug=args.debug,
+            save=plotpath,
+            form=args.form)
 
         # Average spectra for good days
         stanoise.average_sta_spectra(
             fig_average=args.fig_average,
-            save=plotpath, form=args.form)
+            save=plotpath,
+            form=args.form)
 
         if args.fig_av_cross:
             fname = stkey + '.' + 'av_coherence'
             plot = plotting.fig_av_cross(
-                stanoise.f, coh, stanoise.gooddays,
-                'Coherence', stanoise.ncomp, key=stkey, lw=0.5)
-            # if plotpath.is_dir():
+                stanoise.f,
+                coh, stanoise.gooddays,
+                'Coherence',
+                stanoise.ncomp,
+                key=stkey,
+                lw=0.5)
+
             if plotpath:
                 plot.savefig(
                     str(plotpath / (fname + '.' + args.form)),
@@ -550,8 +567,13 @@ def main(args=None):
 
             fname = stkey + '.' + 'av_admittance'
             plot = plotting.fig_av_cross(
-                stanoise.f, ad, stanoise.gooddays,
-                'Admittance', stanoise.ncomp, key=stkey, lw=0.5)
+                stanoise.f,
+                ad,
+                stanoise.gooddays,
+                'Admittance',
+                stanoise.ncomp,
+                key=stkey,
+                lw=0.5)
 
             if plotpath:
                 plot.savefig(
@@ -562,8 +584,14 @@ def main(args=None):
 
             fname = stkey + '.' + 'av_phase'
             plot = plotting.fig_av_cross(
-                stanoise.f, ph, stanoise.gooddays,
-                'Phase', stanoise.ncomp, key=stkey, marker=',', lw=0)
+                stanoise.f,
+                ph,
+                stanoise.gooddays,
+                'Phase',
+                stanoise.ncomp,
+                key=stkey,
+                marker=',',
+                lw=0)
 
             if plotpath:
                 plot.savefig(
@@ -572,9 +600,15 @@ def main(args=None):
             else:
                 plot.show()
 
-        if args.fig_coh_ph and stanoise.phi is not None:
+        if args.fig_tilt and stanoise.phi is not None:
             fname = stkey + '.' + 'coh_ph'
-            plot = plotting.fig_coh_ph(coh_all, ph_all, stanoise.phi)
+            plot = plotting.fig_coh_ph(
+                coh_all,
+                ph_all,
+                stanoise.phi,
+                tilt_list,
+                date_list)
+
             if plotpath:
                 plot.savefig(
                     str(plotpath / (fname + '.' + args.form)),
@@ -583,8 +617,23 @@ def main(args=None):
                 plot.show()
 
         # Save to file
+        print()
+        print("* Clean station spectra saved to: ")
+        print("*   "+str(fileavst))
         stanoise.save(fileavst)
 
+        # Write out events
+        print()
+        print("* Tilt direction as function of time saved to: ")
+        print("*   "+str(filetilt))
+        print()
+        fid = open(filetilt, 'w')
+        fid.writelines("Date, Tilt azimuth (deg from H1)\n")
+        for i in range(len(tilt_list)):
+            line1 = "{0},{1:.0f}\n".format(date_list[i], tilt_list[i])
+            fid.writelines(line1.replace(" ", ""))
+        fid.close()
+
 
 if __name__ == "__main__":