Adding tests for waveform_analytics.py

seismic/inventory/response.py

@@ -113,7 +113,7 @@ def __init__(self, pzTransferFunctionType='LAPLACE (RADIANS/SECOND)',
                      stageGainFreq=1e-2,
                      poles=[0 + 0j],
                      zeros=[0 + 0j]):
-            self.m_base = u'''<?xml version="1.0" encoding="UTF-8" standalone="yes"?>
+            self.m_base = '''<?xml version="1.0" standalone="yes"?>
             <FDSNStationXML
                 xmlns="http://www.fdsn.org/xml/station/1" schemaVersion="1">
                 <Source>-i</Source>
@@ -205,7 +205,6 @@ def __init__(self, pzTransferFunctionType='LAPLACE (RADIANS/SECOND)',
                 </Network>
             </FDSNStationXML>
             '''
-
             self.m_response = None
             self.m_pzTransferFunctionType = pzTransferFunctionType
             self.m_normFactor = normFactor
@@ -223,7 +222,6 @@ def __init__(self, pzTransferFunctionType='LAPLACE (RADIANS/SECOND)',
             # Fetch the response object and adapt its parameters based on user-input.
             self.m_response = inv.get_response('IU.ANMO.00.BHZ', datetime)
 
-            self.m_response.instrument_sensitivity = None  # Get rid of redundant information
             self.m_response.response_stages[0].pz_transfer_function_type = self.m_pzTransferFunctionType
             self.m_response.response_stages[0].normalization_factor = self.m_normFactor
             self.m_response.response_stages[0].normalization_frequency = self.m_normFreq

tests/test_seismic/ASDFdatabase/test_waveform_analytics.py

@@ -0,0 +1,127 @@
+#!/bin/env python
+"""
+Description:
+    Tests various aspects of the functionality in waveform_analytics.py
+
+References:
+
+CreationDate:   08/04/24
+Developer:      [email protected]
+
+Revision History:
+    LastUpdate:     08/04/24   RH
+    LastUpdate:     dd/mm/yyyy  Who     Optional description
+"""
+
+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.inventory.response import ResponseFactory
+from shutil import rmtree
+from scipy.interpolate import interp1d
+import sys
+
+path = os.path.dirname(os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__)))))
+path = os.path.join(path, 'seismic/inventory/')
+resp_file = os.path.join(path, 'RESP.COMPACT120S.MINIMUS.txt')
+tempdir = tempfile.mkdtemp()
+
+def test_fast_psd():
+    network = 'AA'
+    station = 'BB'
+    location = ''
+    channel = 'BHZ'
+    sampling_rate = 10
+    scaling_factor = 1
+
+    def get_time_range_func(net, sta, loc, cha):
+        return UTCDateTime(0), UTCDateTime(0)
+    # end func
+
+    def get_waveforms_func(net, sta, loc, cha, st, et):
+        np.random.seed(10)
+        t = Trace(data=np.random.random(86400 * sampling_rate) * scaling_factor,
+              header={'network':net, 'station':sta, 'location':loc,
+                      'channel':cha, 'sampling_rate':sampling_rate,
+                      'starttime':UTCDateTime('2007-01-01')})
+        return Stream([t])
+    # end func
+
+    # create a flat response
+    rf = ResponseFactory()
+    rf.CreateFromPAZ('flat_response', 'LAPLACE (RADIANS/SECOND)', normFactor=1,
+                      normFreq=1,
+                      stageGain=1,
+                      stageGainFreq=1,
+                      poles=[0 + 1j],
+                      zeros=[0 + 1j])
+    resp = rf.getResponse('flat_response')
+
+    # generate spectrum for fast_ppsd
+    manager = Manager()
+    prog_tracker = ProgressTracker(manager)
+    sa = StationAnalytics(get_time_range_func,
+                          get_waveforms_func,
+                          prog_tracker,
+                          network,
+                          station,
+                          location,
+                          channel,
+                          sampling_rate,
+                          resp,
+                          output_folder = tempdir,
+                          start_time = None,
+                          end_time = None,
+                          nproc=1)
+
+    is_windows = sys.platform.startswith('win')
+
+    # create an interpolation object for the spectrum from fast_ppsd
+    fast_ppsd_io = None
+    for start_time, end_time in zip(sa.st_list, sa.et_list):
+        output_fn_stem = '{}.{}.{}.{}.{}.{}'.format(sa.network,
+                                                    sa.station,
+                                                    sa.location,
+                                                    sa.channel,
+                                                    start_time,
+                                                    end_time)
+
+        if(is_windows): output_fn_stem = output_fn_stem.replace(':', '__')
+        output_fn_npz = os.path.join(sa.output_folder, output_fn_stem + '.npz')
+
+        results = np.load(output_fn_npz)
+        spec = results['sparse_spec']
+
+        fast_ppsd_io = interp1d(sa.sparse_periods, spec)
+    # end for
+
+    # generate spectrum using obspy ppsdf, with a flat response
+    stream = get_waveforms_func(network, station, location, channel, None, None)
+    ppsd = PPSD(stream[0].stats, {'sensitivity': 1.0,
+                                  'gain': 1.0,
+                                  'poles': [0 + 1j],
+                                  'zeros': [0 + 1j]},
+                db_bins=(-150, 50, 1.))
+    ppsd.add(stream)
+
+    periods, spec = ppsd.get_mean()
+    ppsd_io = interp1d(periods, spec)
+
+    # check conformity of spectra
+    common_periods = np.linspace(np.max([sa.sparse_periods[0], periods[0]]),
+                             np.min([sa.sparse_periods[-1], periods[-1]]), 100)
+
+    corr = np.corrcoef(fast_ppsd_io(common_periods), ppsd_io(common_periods))[0,1]
+    assert corr > 0.9
+
+    rmtree(tempdir)
+# end func