Fixed computation of spectral deviation

seismic/ASDFdatabase/analytics/sun.py

@@ -178,6 +178,6 @@ def get_overlap(a, b):
             # end if
     # end for
 
-    return day_coverage_seconds / day_length_seconds, \
-           night_coverage_seconds / night_length_seconds
+    return day_coverage_seconds / day_length_seconds if day_coverage_seconds > 0 else 0, \
+           night_coverage_seconds / night_length_seconds if night_length_seconds > 0 else 0
 # end func

seismic/ASDFdatabase/analytics/waveform_analytics.py

@@ -134,7 +134,6 @@ def __init__(self,
         self.freqs[0] *= -1
         self.periods = 1. / self.freqs
         self.sparse_periods = None
-        self.sparse_period_indices = None
         self.nm_periods, self.lnm = get_nlnm()
         _, self.hnm = get_nhnm()
         self.AMP_DB_MIN = -200
@@ -147,7 +146,6 @@ def __init__(self,
                                                        (1, 1, 0),
                                                        (0, 1, 0)], N=5)
         self.cmap_norm = matplotlib.colors.Normalize(vmin=0, vmax=100)
-        self.overlap_denominator = None
 
         # flip noise-model values and restrict them to data range
         self.nm_periods = self.nm_periods[::-1]
@@ -192,7 +190,7 @@ def __init__(self,
         self.progress_tracker.initialize(len(self.st_list))
 
         # launch parallel computations
-        if (0):
+        if (1):
             p = Pool(ncpus=self.nproc)
             p.map(self._generate_psds, proc_st_list, proc_et_list)
         else:
@@ -233,14 +231,14 @@ def _setup_period_bins(self):
 
         # initialize nearest-neighbour indices for mapping:
         # i) raw PSD periods to sparse-periods
-        # ii) sparse-periods to noise-model periods
+        # ii) noise-model periods to sparse-periods
         self.dense_to_sparse_indices = np.fabs(self.sparse_periods[:, None] -
                                                self.periods[None, :]).argmin(axis=-1)
-        self.sparse_to_nm_indices = np.fabs(self.nm_periods[None, :] -
-                                            self.sparse_periods[:, None]).argmin(axis=-1)
-        self.overlap_denominator = self.sparse_periods[
-            np.where((self.sparse_periods >= self.nm_periods[0]) & \
-                     (self.sparse_periods <= self.nm_periods[-1]))].shape[0]
+        self.nm_to_sparse_indices = np.fabs(self.sparse_periods[:, None] -
+                                            self.nm_periods[None, :]).argmin(axis=-1)
+
+        self.sparse_nm_overlap_indices = np.where((self.sparse_periods >= self.nm_periods[0]) &
+                                                  (self.sparse_periods <= self.nm_periods[-1]))[0]
     # end func
 
     def _generate_psds(self, start_time_list, end_time_list):
@@ -327,14 +325,20 @@ def fft_taper(data):
             # compare spec to low- and high-noise-model
             sparse_spec = spec[self.dense_to_sparse_indices]
 
-            sparse_spec_deviation = np.where((self.lnm[self.sparse_to_nm_indices] > sparse_spec) | \
-                                             (self.hnm[self.sparse_to_nm_indices] < sparse_spec))[0]
+            sparse_spec_deviation_indices = np.where(
+                (self.lnm[self.nm_to_sparse_indices] > sparse_spec) | \
+                (self.hnm[self.nm_to_sparse_indices] < sparse_spec))[0]
+
+            # Only include overlapping periods between nm and sparse_spec while calculating
+            # spectral deviation
+            overlap_deviation_indices = np.intersect1d(sparse_spec_deviation_indices,
+                                                       self.sparse_nm_overlap_indices)
+            deviation_fraction = overlap_deviation_indices.shape[0] /\
+                                 self.sparse_nm_overlap_indices.shape[0]
 
             ############################################
             # Plot results and write output npz
             ############################################
-            deviation_fraction = len(sparse_spec_deviation) / self.overlap_denominator
-
             output_fn_stem = '{}.{}.{}.{}.{}.{}'.format(self.network,
                                                         self.station,
                                                         self.location,
@@ -367,18 +371,10 @@ def fft_taper(data):
                     ax.semilogx(self.nm_periods, self.lnm, 'k', lw=1)
                     ax.semilogx(self.nm_periods, self.hnm, 'k', lw=1)
 
-                    # ax.semilogx(self.sparse_periods[sparse_spec_deviation],
-                    #            sparse_spec[sparse_spec_deviation], 'b',
-                    #            linestyle=None, marker='+')
-
                     ax.set_ylim(self.AMP_DB_MIN, self.AMP_DB_MAX)
                     ax.set_xlim(self.PER_MIN, np.max(self.periods))
                     ax.grid(True, which="both", ls="-", lw=0.2)
                     ax.tick_params(labelsize=6)
-                    # ax.set_xlabel("Period [s]", fontsize=5)
-                    # ax.set_ylabel("Amp. [m2/s4][dB]", fontsize=5)
-                    # ax.xaxis.set_label_coords(0.5, -0.04)
-                    # ax.yaxis.set_label_coords(1.05, 0.5)
 
                     ax.text(0.015, -180,
                             'Coverage: {:.2f} %'.format(

seismic/ASDFdatabase/utils.py

@@ -186,6 +186,8 @@ def __init__(self, mseed_folder, pattern):
         self.tree = defaultdict(lambda: defaultdict(lambda: defaultdict(lambda: defaultdict(list))))
         self.stream_cache = MseedIndex.StreamCache()
         self.mseed_files = np.array(sorted(glob(os.path.join(self.mseed_folder, pattern))))
+        self.sr_coverage_dict = defaultdict(float) # dict keyed by unique sampling rates found, with coverage as values
+        self.ch_coverage_dict = defaultdict(float) # dict keyed by unique channels found, with coverage as values
 
         fc = len(self.mseed_files)
         if(fc > 0):
@@ -209,14 +211,19 @@ def __init__(self, mseed_folder, pattern):
             # end try
 
             for tr in st:
-                nc, sc, lc, cc, st, et = \
+                nc, sc, lc, cc, st, et, sr = \
                     tr.stats.network, tr.stats.station, tr.stats.location, \
                         tr.stats.channel, tr.stats.starttime.timestamp, \
-                        tr.stats.endtime.timestamp
+                        tr.stats.endtime.timestamp, tr.stats.sampling_rate
 
                 # skip bogus traces
                 if(nc == sc == lc == cc == ''): continue
                 self.meta_list.append([i, nc, sc, lc, cc, st, et])
+
+                # store coverage for unique channels and sampling rates
+                cov = float(et - st) # coverage in seconds
+                self.sr_coverage_dict[sr] += cov
+                self.ch_coverage_dict[cc] += cov
             # end for
             # if (i > 0): break
         # end for
@@ -233,6 +240,14 @@ def __init__(self, mseed_folder, pattern):
             # end for
     # end func
 
+    def get_sampling_rate_coverage(self) -> defaultdict:
+        return self.sr_coverage_dict
+    # end func
+
+    def get_channel_coverage(self) -> defaultdict:
+        return self.ch_coverage_dict
+    # end func
+
     def __getstate__(self):
         #print('pickling..')
         return self.__dict__

tests/test_seismic/ASDFdatabase/test_waveform_analytics.py

@@ -14,17 +14,13 @@
 """
 
 import os
-import pytest
-from ordered_set import OrderedSet as set
-import numpy as np
 import tempfile
-import sqlite3
 import numpy as np
 from obspy.core import Trace, Stream
 from obspy.signal.spectral_estimation import PPSD
 from obspy import UTCDateTime
 from multiprocessing import Manager
-from seismic.ASDFdatabase.waveform_analytics import ProgressTracker, StationAnalytics
+from seismic.ASDFdatabase.analytics.waveform_analytics import ProgressTracker, StationAnalytics
 from seismic.inventory.response import ResponseFactory
 from shutil import rmtree
 from scipy.interpolate import interp1d
@@ -44,7 +40,7 @@ def test_fast_psd():
     scaling_factor = 1
 
     def get_time_range_func(net, sta, loc, cha):
-        return UTCDateTime(0), UTCDateTime(0)
+        return UTCDateTime('2007-01-01'), UTCDateTime('2007-01-02')
     # end func
 
     def get_waveforms_func(net, sta, loc, cha, st, et):