Merge pull request #36 from nfsi-canada/scripts

Script verbosity, unit tests, windowed fft

Author: paudetseis

Diff for: .github/workflows/deploy.yml

@@ -17,18 +17,21 @@ on:
 jobs:
   # This workflow contains a single job called "build"
   build:
+    name: Python ${{ matrix.python-version }}
     # The type of runner that the job will run on
     runs-on: ${{ matrix.os }}
     strategy:
       matrix:
         os: [ubuntu-latest]
-        python-version: [3.7, 3.8, 3.9]
+        python-version: ['3.7', '3.8', '3.9']
 
     steps:
       # Checks-out your repository under $GITHUB_WORKSPACE, so your job can access it
       - uses: actions/checkout@v2
       - uses: conda-incubator/setup-miniconda@v2
         with:
+          auto-update-conda: true
+          python-version: ${{ matrix.python-version }}
           activate-environment: test
           auto-activate-base: false
 
@@ -44,6 +47,7 @@ jobs:
       - name: Tests
         shell: bash -l {0}
         run: |
+          export MPLBACKEND=Agg
           mkdir empty
           cd empty
           conda install pytest-cov

Diff for: docs/atacr.rst

@@ -731,7 +731,7 @@ M08A and send the prompt to a logfile
 
 .. code-block::
 
-   $ query_fdsn_stdb.py -N 7D -C ?H? -S M08A M08A > logfile
+   $ query_fdsn_stdb -N 7D -C ?H? -S M08A M08A > logfile
 
 .. note::
 
@@ -743,11 +743,11 @@ M08A and send the prompt to a logfile
     enclose the `?H?` in single or double quotes (e.g., 
     `query_fdsn_stdb.py -N 7D -C "?H?" -S M08A M08A > logfile`) 
 
-To check the station info for M08A, use the program ``ls_stdb.py``:
+To check the station info for M08A, use the program ``ls_stdb``:
 
 .. code-block::
 
-   $ ls_stdb.py M08A.pkl
+   $ ls_stdb M08A.pkl
     Listing Station Pickle: M08A.pkl
     7D.M08A
     --------------------------------------------------------------------------
@@ -1144,7 +1144,7 @@ Vanuatu earthquake (be conservative with the options), type in a terminal:
     |  Found     1 possible events                  |
      
     ****************************************************
-    * #1 (2/1):  20120309_070953
+    * (1/1):  20120309_070953
     *   Origin Time: 2012-03-09 07:09:53
     *   Lat: -19.22; Lon:  169.75
     *   Dep:  33.70; Mag: 6.6

Diff for: obstools/__init__.py

@@ -109,7 +109,7 @@
 
 """
 
-__version__ = '0.1.2'
+__version__ = '0.1.3'
 
 __author__ = 'Pascal Audet & Helen Janiszewski'
 

Diff for: obstools/atacr/classes.py

@@ -22,7 +22,7 @@
 
 
 import sys
-from scipy.signal import spectrogram, detrend
+from scipy.signal import spectrogram, stft, detrend
 from scipy.linalg import norm
 import matplotlib.pyplot as plt
 import numpy as np
@@ -349,6 +349,8 @@ def QC_daily_spectra(self, pd=[0.004, 0.2], tol=1.5, alpha=0.05,
         psdP = None
         f, t, psdZ = spectrogram(
             self.trZ.data, self.fs, window=wind, nperseg=ws, noverlap=ss)
+        self.f = f
+        print(f)
         if self.ncomp == 2 or self.ncomp == 4:
             f, t, psdP = spectrogram(
                 self.trP.data, self.fs, window=wind, nperseg=ws, noverlap=ss)
@@ -640,20 +642,37 @@ def average_daily_spectra(self, calc_rotation=True, fig_average=False,
         ws = int(self.window/self.dt)
 
         # Number of points in step
-        ss = int(self.window*(1.-self.overlap)/self.dt)
+        ss = int(self.window*self.overlap/self.dt)
+
+        # hanning window
+        hanning = np.hanning(2*ss)
+        wind = np.ones(ws)
+        wind[0:ss] = hanning[0:ss]
+        wind[-ss:ws] = hanning[ss:ws]
 
         ft1 = None
         ft2 = None
         ftZ = None
         ftP = None
-        ftZ, f = utils.calculate_windowed_fft(self.trZ, ws, ss)
+
+        _f, _t, ftZ = stft(
+            self.trZ.data, self.fs, return_onesided=False, boundary=None, padded=False, 
+            window=wind, nperseg=ws, noverlap=ss)
+        ftZ = ftZ.T
         if self.ncomp == 2 or self.ncomp == 4:
-            ftP, f = utils.calculate_windowed_fft(self.trP, ws, ss)
+            _f, _t, ftP = stft(
+                self.trP.data, self.fs, return_onesided=False, boundary=None, padded=False, 
+                window=wind, nperseg=ws, noverlap=ss)
+            ftP = ftP.T
         if self.ncomp == 3 or self.ncomp == 4:
-            ft1, f = utils.calculate_windowed_fft(self.tr1, ws, ss)
-            ft2, f = utils.calculate_windowed_fft(self.tr2, ws, ss)
-
-        self.f = f
+            _f, _t, ft1 = stft(
+                self.tr1.data, self.fs, return_onesided=False, boundary=None, padded=False, 
+                window=wind, nperseg=ws, noverlap=ss)
+            _f, _t, ft2 = stft(
+                self.tr2.data, self.fs, return_onesided=False, boundary=None, padded=False, 
+                window=wind, nperseg=ws, noverlap=ss)
+            ft1 = ft1.T
+            ft2 = ft2.T
 
         # Extract good windows
         c11 = None
@@ -662,18 +681,18 @@ def average_daily_spectra(self, calc_rotation=True, fig_average=False,
         cPP = None
         cZZ = np.abs(
             np.mean(ftZ[self.goodwins, :]*np.conj(ftZ[self.goodwins, :]),
-                    axis=0))[0:len(f)]
+                    axis=0))[0:len(self.f)]
         if self.ncomp == 2 or self.ncomp == 4:
             cPP = np.abs(
                 np.mean(ftP[self.goodwins, :]*np.conj(ftP[self.goodwins, :]),
-                        axis=0))[0:len(f)]
+                        axis=0))[0:len(self.f)]
         if self.ncomp == 3 or self.ncomp == 4:
             c11 = np.abs(
                 np.mean(ft1[self.goodwins, :]*np.conj(ft1[self.goodwins, :]),
-                        axis=0))[0:len(f)]
+                        axis=0))[0:len(self.f)]
             c22 = np.abs(
                 np.mean(ft2[self.goodwins, :]*np.conj(ft2[self.goodwins, :]),
-                        axis=0))[0:len(f)]
+                        axis=0))[0:len(self.f)]
 
         # Extract bad windows
         bc11 = None
@@ -683,18 +702,18 @@ def average_daily_spectra(self, calc_rotation=True, fig_average=False,
         if np.sum(~self.goodwins) > 0:
             bcZZ = np.abs(np.mean(
                 ftZ[~self.goodwins, :]*np.conj(ftZ[~self.goodwins, :]),
-                axis=0))[0:len(f)]
+                axis=0))[0:len(self.f)]
             if self.ncomp == 2 or self.ncomp == 4:
                 bcPP = np.abs(np.mean(
                     ftP[~self.goodwins, :]*np.conj(ftP[~self.goodwins, :]),
-                    axis=0))[0:len(f)]
+                    axis=0))[0:len(self.f)]
             if self.ncomp == 3 or self.ncomp == 4:
                 bc11 = np.abs(np.mean(
                     ft1[~self.goodwins, :]*np.conj(ft1[~self.goodwins, :]),
-                    axis=0))[0:len(f)]
+                    axis=0))[0:len(self.f)]
                 bc22 = np.abs(np.mean(
                     ft2[~self.goodwins, :]*np.conj(ft2[~self.goodwins, :]),
-                    axis=0))[0:len(f)]
+                    axis=0))[0:len(self.f)]
 
         # Calculate mean of all good windows if component combinations exist
         c12 = None
@@ -705,27 +724,27 @@ def average_daily_spectra(self, calc_rotation=True, fig_average=False,
         cZP = None
         if self.ncomp == 3 or self.ncomp == 4:
             c12 = np.mean(ft1[self.goodwins, :] *
-                          np.conj(ft2[self.goodwins, :]), axis=0)[0:len(f)]
+                          np.conj(ft2[self.goodwins, :]), axis=0)[0:len(self.f)]
             c1Z = np.mean(ft1[self.goodwins, :] *
-                          np.conj(ftZ[self.goodwins, :]), axis=0)[0:len(f)]
+                          np.conj(ftZ[self.goodwins, :]), axis=0)[0:len(self.f)]
             c2Z = np.mean(ft2[self.goodwins, :] *
-                          np.conj(ftZ[self.goodwins, :]), axis=0)[0:len(f)]
+                          np.conj(ftZ[self.goodwins, :]), axis=0)[0:len(self.f)]
         if self.ncomp == 4:
             c1P = np.mean(ft1[self.goodwins, :] *
-                          np.conj(ftP[self.goodwins, :]), axis=0)[0:len(f)]
+                          np.conj(ftP[self.goodwins, :]), axis=0)[0:len(self.f)]
             c2P = np.mean(ft2[self.goodwins, :] *
-                          np.conj(ftP[self.goodwins, :]), axis=0)[0:len(f)]
+                          np.conj(ftP[self.goodwins, :]), axis=0)[0:len(self.f)]
         if self.ncomp == 2 or self.ncomp == 4:
             cZP = np.mean(ftZ[self.goodwins, :] *
-                          np.conj(ftP[self.goodwins, :]), axis=0)[0:len(f)]
+                          np.conj(ftP[self.goodwins, :]), axis=0)[0:len(self.f)]
 
         # Store as attributes
         self.power = Power(c11, c22, cZZ, cPP)
         self.cross = Cross(c12, c1Z, c1P, c2Z, c2P, cZP)
         bad = Power(bc11, bc22, bcZZ, bcPP)
 
         if fig_average:
-            plot = plotting.fig_average(f, self.power, bad, self.goodwins,
+            plot = plotting.fig_average(self.f, self.power, bad, self.goodwins,
                                         self.ncomp, key=self.key)
             if save:
                 fname = self.key + '.' + self.tkey + '.' + 'average.' + form
@@ -739,7 +758,7 @@ def average_daily_spectra(self, calc_rotation=True, fig_average=False,
         if calc_rotation and self.ncomp >= 3:
             cHH, cHZ, cHP, coh, ph, direc, tilt, coh_value, phase_value = \
                 utils.calculate_tilt(
-                    ft1, ft2, ftZ, ftP, f, self.goodwins)
+                    ft1, ft2, ftZ, ftP, self.f, self.goodwins)
             self.rotation = Rotation(
                 cHH, cHZ, cHP, coh, ph, tilt, coh_value, phase_value, direc)
 
@@ -1949,12 +1968,15 @@ def correct_data(self, tfnoise):
         ft2 = None
         ftZ = None
         ftP = None
-        ftZ, f = utils.calculate_windowed_fft(trZ, ws, hann=False)
+        ftZ = np.fft.fft(trZ, n=ws)
         if self.ncomp == 2 or self.ncomp == 4:
-            ftP, f = utils.calculate_windowed_fft(trP, ws, hann=False)
+            ftP = np.fft.fft(trP, n=ws)
         if self.ncomp == 3 or self.ncomp == 4:
-            ft1, f = utils.calculate_windowed_fft(tr1, ws, hann=False)
-            ft2, f = utils.calculate_windowed_fft(tr2, ws, hann=False)
+            ft1 = np.fft.fft(tr1, n=ws)
+            ft2 = np.fft.fft(tr2, n=ws)
+
+        # Use one-sided frequency axis to match spectrogram
+        f = np.fft.rfftfreq(ws, d=self.dt)
 
         if not np.allclose(f, tfnoise.f):
             raise(Exception('Frequency axes are different: ', f, tfnoise.f))

Diff for: obstools/scripts/atacr_clean_spectra.py

@@ -317,8 +317,7 @@ def main(args=None):
         sta.location = tlocs
 
         # Update Display
-        print()
-        print("|===============================================|")
+        print("\n|===============================================|")
         print("|===============================================|")
         print("|                   {0:>8s}                    |".format(
             sta.station))
@@ -385,8 +384,7 @@ def main(args=None):
 
             # Load file if it exists
             if filespec.exists():
-                print()
-                print("*"*60)
+                print("\n"+"*"*60)
                 print('* Calculating noise spectra for key ' +
                       stkey+' and day '+year+'.'+jday)
                 print("*   -> file "+str(filespec)+" found - loading")

Diff for: obstools/scripts/atacr_daily_spectra.py

@@ -321,15 +321,13 @@ def main(args=None):
         # Path where data are located
         datapath = Path('DATA') / stkey
         if not datapath.is_dir():
-            print()
-            print("Path to "+str(datapath)+" doesn`t exist - continuing")
+            print("\nPath to "+str(datapath)+" doesn`t exist - continuing")
             continue
 
         # Path where spectra will be saved
         specpath = Path('SPECTRA') / stkey
         if not specpath.is_dir():
-            print()
-            print("Path to "+str(specpath)+" doesn`t exist - creating it")
+            print("\nPath to "+str(specpath)+" doesn`t exist - creating it")
             specpath.mkdir(parents=True)
 
         # Path where plots will be saved
@@ -365,8 +363,7 @@ def main(args=None):
         sta.location = tlocs
 
         # Update Display
-        print()
-        print("|===============================================|")
+        print("\n|===============================================|")
         print("|===============================================|")
         print("|                   {0:>8s}                    |".format(
             sta.station))
@@ -403,8 +400,7 @@ def main(args=None):
             year = str(trZ.stats.starttime.year).zfill(4)
             jday = str(trZ.stats.starttime.julday).zfill(3)
 
-            print()
-            print("*"*60)
+            print("\n"+"*"*60)
             print("* Calculating noise spectra for key " +
                   stkey+" and day "+year+"."+jday)
             tstamp = year+'.'+jday+'.'

Diff for: obstools/scripts/atacr_download_data.py

@@ -318,8 +318,7 @@ def main(args=None):
         # Define path to see if it exists
         datapath = Path('DATA') / Path(stkey)
         if not datapath.is_dir():
-            print()
-            print('Path to '+str(datapath)+' doesn`t exist - creating it')
+            print('\nPath to '+str(datapath)+' doesn`t exist - creating it')
             datapath.mkdir(parents=True)
 
         # Establish client
@@ -354,8 +353,7 @@ def main(args=None):
         sta.location = tlocs
 
         # Update Display
-        print()
-        print("|===============================================|")
+        print("\n|===============================================|")
         print("|===============================================|")
         print("|                   {0:>8s}                    |".format(
             sta.station))
@@ -389,9 +387,7 @@ def main(args=None):
             # Time stamp
             tstamp = str(t1.year).zfill(4)+'.'+str(t1.julday).zfill(3)+'.'
 
-            print()
-            print(
-                "***********************************************************")
+            print("\n"+"*"*60)
             print("* Downloading day-long data for key "+stkey +
                   " and day "+str(t1.year)+"."+str(t1.julday))
             print("*")