added coherence subplot in tilt figure, cleaned up code base, and updated documentation

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

@@ -267,31 +267,37 @@ def fig_coh_ph(coh, ph, direc, tilt, date):
         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])
+        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(
-            np.mean(tilt),
+            meantilt,
             ls='--',
-            label='Mean: {0:.0f} $\pm$ {1:.1f}'.format(meantilt, stdetilt))
+            label=r'Mean: {0:.0f} $\pm$ {1:.1f}'.format(meantilt, stdetilt))
         ax3.axhline(
-            np.median(tilt),
+            mediantilt,
             ls=':',
-            label='Median: {0:.0f}'.format(np.median(tilt)))
+            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])
             ax3.plot(d, t, 'o', c=colors[i])
-        ax3.legend(loc='best', fontsize=8)
+            ax4.plot(d, mcoh, 'o', c=colors[i])
         ax1.set_ylabel('Coherence Z-H1')
         ax1.set_ylim((0, 1.))
-        ax2.set_ylabel('Phase Z-H1')
-        ax3.set_ylabel('Tilt direction (deg)')
         ax1.set_xlabel('Angle clockwise from H1 (deg)')
+        ax2.set_ylabel('Phase Z-H1')
         ax2.set_xlabel('Angle clockwise from H1 (deg)')
-        plt.tight_layout()
+        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]
@@ -451,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)
 
@@ -489,6 +491,8 @@ def calculate_tilt(ft1, ft2, ftZ, ftP, f, goodwins, tiltfreqs):
         tilt -= 360.
     elif tilt < 0:
         tilt += 360.
+    if tilt > 180.:
+        tilt -= 180.
 
     # Now calculate spectra at tilt direction
     ftH = rotate_dir(ft1, ft2, tilt)
@@ -617,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

@@ -166,7 +166,7 @@ def get_cleanspec_arguments(argv=None):
         action="store_true",
         dest="fig_tilt",
         default=False,
-        help="Plot Coherence, phase and tilt direction figure. " +
+        help="Plot coherence, phase and tilt direction figure. " +
         "[Default does not plot figure]")
     FigureGroup.add_argument(
         "--figCross",
@@ -349,7 +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')
+        filetilt = avstpath / (dstart+dend+'tilt.csv')
 
         if fileavst.exists():
             if not args.ovr:
@@ -387,6 +387,7 @@ def main(args=None):
         # Date + tilt list
         date_list = []
         tilt_list = []
+        coh_list = []
 
         # Loop through each day withing time range
         while t1 < tend:
@@ -408,6 +409,7 @@ def main(args=None):
                 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.
@@ -624,13 +626,14 @@ def main(args=None):
 
         # Write out events
         print()
-        print("* Tilt direction as function of time saved to: ")
+        print("* Tilt direction and coherence as function of time saved to: ")
         print("*   "+str(filetilt))
         print()
         fid = open(filetilt, 'w')
-        fid.writelines("Date, Tilt azimuth (deg from H1)\n")
+        fid.writelines("Date, Tilt dir. (deg from H1), Max coherehce\n")
         for i in range(len(tilt_list)):
-            line1 = "{0},{1:.0f}\n".format(date_list[i], tilt_list[i])
+            line1 = "{0},{1:.0f},{2:.2f}\n".format(
+                date_list[i], tilt_list[i], coh_list[i])
             fid.writelines(line1.replace(" ", ""))
         fid.close()
 

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'