Merge pull request #49 from nfsi-canada/v021

V021

Author: paudetseis

docs/atacr.rst

@@ -487,7 +487,7 @@ Usage
       --debug          Plot intermediate steps for debugging. [Default does not
                        plot figure]
       --figAverage     Plot daily average figure. [Default does not plot figure]
-      --figCoh         Plot Coherence and Phase figure. [Default does not plot
+      --figTilt        Plot coherence, phase and tilt direction figure. [Default does not plot
                        figure]
       --figCross       Plot cross-spectra figure. [Default does not plot figure]
       --save-fig       Set this option if you wish to save the figure(s). [Default
@@ -1034,7 +1034,7 @@ therefore using all available data) and plot the results, we can type in a termi
 
 .. code-block::
 
-    $ atacr_clean_spectra --figQC --figAverage --figCoh --figCross M08A.pkl
+    $ atacr_clean_spectra --figQC --figAverage --figTilt --figCross M08A.pkl
 
     ###################################################################
     #       _                                          _              #
@@ -1071,8 +1071,13 @@ therefore using all available data) and plot the results, we can type in a termi
 
     ...
 
-And so on until all ``DayNoise`` objects are averaged into a ``StaNoise`` 
-object, which is saved to a newly created folder called ``AVG_STA/7D.M08A/``.
+    * Clean station spectra saved to: 
+    *   AVG_STA/7D.M08A/2011.293-2012.200.avg_sta.pkl
+
+    * Tilt direction as function of time saved to: 
+    *   AVG_STA/7D.M08A/2011.293-2012.200.tilt.csv
+
+The new ``StaNoise`` object and the daily tilt directions are saved to a newly created folder called ``AVG_STA/7D.M08A/``.
 
 .. note ::
 
@@ -1095,9 +1100,13 @@ Several figures are also produced, including Figures 4, 6-9.
 
 Figure 4: The orientation of maximum coherence between the vertical and the two
 horizontal components for M08A during March 2012. (Left) Coherence as a function
-of angle from the H1 component. (Right) Phase as a function of the angle. In this
-example, the coherence is low indicating the absence of dominant, uni-directional tilt
-noise.
+of angle from the H1 component. (Right) Phase as a function of the angle. 
+(Bottom) Tilt direction measured from the H1 component as a function of time. 
+In this example, the coherence is low indicating low tilt noise.
+
+.. note ::
+
+    In the MATLAB version of ``ATaCR``, the "Angle from H1" (Figure 4 in the manual) corresponds to a counter-clockwise rotation of H1 and H2, and the orientation is chosen at the lowest phase. Here, the angle corresponds to a clockwise rotation and the estimate is chosen to lie between 0 and 180 degrees. There is therefore a 225-degree difference between those two estimates, and the plots are mirrored.
 
 .. figure:: ../obstools/examples/figures/Figure_6.png
    :align: center

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,22 +251,53 @@ 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(layout='constrained')
+        gs = GridSpec(4, 2, figure=f)
+        ax1 = f.add_subplot(gs[0:-2, 0])
+        ax2 = f.add_subplot(gs[0:-2, 1])
+        ax3 = f.add_subplot(gs[2, :])
+        ax4 = f.add_subplot(gs[3, :])
+        ax3.axhline(
+            meantilt,
+            ls='--',
+            label=r'Mean: {0:.0f} $\pm$ {1:.1f}'.format(meantilt, stdetilt))
+        ax3.axhline(
+            mediantilt,
+            ls=':',
+            label='Median: {0:.0f}'.format(mediantilt))
+        for i, (co, p, d, t) in enumerate(zip(coh, ph, date, tilt)):
+            mcoh = np.max(co)
             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])
+            ax4.plot(d, mcoh, 'o', c=colors[i])
+        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')
-        plt.tight_layout()
+        ax1.set_xlabel('Angle clockwise from H1 (deg)')
+        ax2.set_ylabel('Phase Z-H1')
+        ax2.set_xlabel('Angle clockwise from H1 (deg)')
+        ax3.legend(loc='best', fontsize=8)
+        ax3.set_xticklabels([])
+        ax3.set_ylabel('Tilt dir. (deg)')
+        ax3.set_ylim(-10, 190)
+        ax4.set_ylabel('Coherence')
+        ax4.set_ylim(-0.1, 1.1)
 
     else:
         plt.figure()

obstools/atacr/utils.py

@@ -304,7 +304,7 @@ def get_event(eventpath, tstart, tend):
     # Find out how many events from Z.SAC files
     eventfiles = list(eventpath.glob('*Z.SAC'))
     if not eventfiles:
-        raise(Exception("No event found in folder "+str(eventpath)))
+        raise Exception("No event found in folder "+str(eventpath))
 
     # Extract events from time stamps
     prefix = [file.name.split('.') for file in eventfiles]
@@ -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):
 
@@ -452,6 +451,8 @@ def calculate_tilt(ft1, ft2, ftZ, ftP, f, goodwins, tiltfreqs):
         tilt += 180.
     if tilt > 360.:
         tilt -= 360.
+    if tilt > 180.:
+        tilt -= 180.
 
     # print('Tilt corrected = ', tilt)
 
@@ -486,6 +487,12 @@ 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.
+    if tilt > 180.:
+        tilt -= 180.
 
     # Now calculate spectra at tilt direction
     ftH = rotate_dir(ft1, ft2, tilt)
@@ -614,16 +621,14 @@ def phase(Gxy):
 
 def rotate_dir(x, y, theta):
 
-    d = -theta*np.pi/180.
-    rot_mat = [[np.cos(d), -np.sin(d)],
-               [np.sin(d), np.cos(d)]]
+    d = theta*np.pi/180.
+    rot_mat = [[np.cos(d), np.sin(d)],
+               [-np.sin(d), np.cos(d)]]
 
     vxy = [x, y]
     vxy_rotated = np.tensordot(rot_mat, vxy, axes=1)
-    tr_2 = vxy_rotated[0]
-    tr_1 = vxy_rotated[1]
 
-    return tr_1
+    return vxy_rotated[1]
 
 
 def ftest(res1, pars1, res2, pars2):

obstools/examples/figures/Figure_4.png

@@ -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.csv')
 
         if fileavst.exists():
             if not args.ovr:
@@ -383,14 +384,19 @@ def main(args=None):
         # Initialize StaNoise object
         stanoise = StaNoise()
 
+        # Date + tilt list
+        date_list = []
+        tilt_list = []
+        coh_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 +407,9 @@ 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)
+                coh_list.append(daynoise.rotation.coh_value)
                 stanoise += daynoise
             else:
                 t1 += 3600.*24.
@@ -526,21 +535,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 +569,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 +586,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 +602,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 +619,24 @@ 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 and coherence as function of time saved to: ")
+        print("*   "+str(filetilt))
+        print()
+        fid = open(filetilt, 'w')
+        fid.writelines("Date, Tilt dir. (deg from H1), Max coherehce\n")
+        for i in range(len(tilt_list)):
+            line1 = "{0},{1:.0f},{2:.2f}\n".format(
+                date_list[i], tilt_list[i], coh_list[i])
+            fid.writelines(line1.replace(" ", ""))
+        fid.close()
+
 
 if __name__ == "__main__":
 

obstools/scripts/atacr_clean_spectra.py

@@ -414,7 +414,7 @@ def main(args=None):
         while t2 <= tend:
 
             # Time stamp
-            tstamp = str(t1.year).zfill(4)+'.'+str(t1.julday).zfill(3)+'.'
+            tstamp = str(t1.year).zfill(4)+'.'+str(t1.julday).zfill(3)
 
             print("\n"+"*"*60)
             print("* Downloading day-long data for key {0} and day {1}.{2:03d}".format(
@@ -441,7 +441,7 @@ def main(args=None):
                     t2 += dt
                     continue
 
-            print("*   "+tstamp +"*.SAC")
+            print("*   " + tstamp + ".*.SAC")
             # Get waveforms from client, one channel at a time
             try:
                 cha = sta.channel.upper() + '1'