Added a parameter to better control time-domain normalization

Author: geojunky

seismic/xcorqc/correlator.py

@@ -43,8 +43,8 @@ def process(data_source1, data_source2, output_path,
             fmin=None, fmax=None, netsta_list1='*', netsta_list2='*', pairs_to_compute=None,
             start_time='1970-01-01T00:00:00', end_time='2100-01-01T00:00:00',
             instrument_response_inventory=None, instrument_response_output='vel', water_level=50,
-            clip_to_2std=False, whitening=False, whitening_window_frequency=0,
-            one_bit_normalize=False, location_preferences=None, ds1_zchan=None, ds1_nchan=None,
+            time_domain_norm='zero_mean_one_std', whitening=False, whitening_window_frequency=0,
+            location_preferences=None, ds1_zchan=None, ds1_nchan=None,
             ds1_echan=None, ds2_zchan=None, ds2_nchan=None, ds2_echan=None, corr_chan=None,
             envelope_normalize=False, ensemble_stack=False, subset_stacker=None, apply_simple_stacking=True,
             restart=False, dry_run=False, no_tracking_tag=False, scratch_folder=None):
@@ -113,8 +113,7 @@ def outputConfigParameters():
             f.write('%35s\t\t\t: %s\n' % ('--instrument-response-output', instrument_response_output))
             f.write('%35s\t\t\t: %s\n' % ('--corr-chan', corr_chan))
             f.write('%35s\t\t\t: %s\n' % ('--water-level', water_level))
-            f.write('%35s\t\t\t: %s\n' % ('--clip-to-2std', clip_to_2std))
-            f.write('%35s\t\t\t: %s\n' % ('--one-bit-normalize', one_bit_normalize))
+            f.write('%35s\t\t\t: %s\n' % ('--time-domain-norm', time_domain_norm))
             f.write('%35s\t\t\t: %s\n' % ('--envelope-normalize', envelope_normalize))
             f.write('%35s\t\t\t: %s\n' % ('--whitening', whitening))
             if(whitening):
@@ -309,8 +308,8 @@ def get_loccha(cha1, cha2):
                                baz_netsta1, baz_netsta2,
                                resample_rate, taper_length, read_ahead_window_seconds,
                                interval_seconds, window_seconds, window_overlap,
-                               window_buffer_length, fmin, fmax, clip_to_2std, whitening,
-                               whitening_window_frequency, one_bit_normalize, envelope_normalize,
+                               window_buffer_length, fmin, fmax, time_domain_norm, whitening,
+                               whitening_window_frequency, envelope_normalize,
                                ensemble_stack, subset_stacker, apply_simple_stacking, output_path, 2,
                                time_tag, scratch_folder, git_hash)
         # end for
@@ -389,19 +388,21 @@ def get_loccha(cha1, cha2):
               help="Water-level in dB to limit amplification during instrument response correction "
                    "to a certain cut-off value. Note, this parameter has no effect if instrument "
                    "response correction is not performed.")
-@click.option('--clip-to-2std', is_flag=True,
-              help="Clip data in each window to +/- 2 standard deviations. Note that the default time-domain normalization "
-                   "is N(0,1), i.e. 0-mean and unit variance")
+@click.option('--time-domain-norm', type=click.Choice(['zero_mean_one_std', 'clip_to_two_std',
+                                                       'one_bit_norm', 'none']),
+              default='zero_mean_one_std',
+              help="Time domain normalization: "
+                   "zero_mean_one_std: the default time-domain normalization, N(0,1), i.e. 0-mean and unit variance"
+                   "clip_to_two_std: clips data in each window to +/- 2 standard deviations. "
+                   "one_bit_norm: apply one-bit normalization to data in each window"
+                   "none: no time domain normalization")
 @click.option('--whitening', is_flag=True,
               help="Apply spectral whitening")
 @click.option('--whitening-window-frequency', type=float, default=0,
               help="Window frequency (Hz) determines the half-window length (of averaging window) used for smoothing weights "
                    "that scale the spectral amplitudes of the waveform being spectrally whitened. The default value of 0 "
                    "implies no smoothing of weights. Note that this parameter has no effect unless whitening is activated with "
                    "'--whitening'")
-@click.option('--one-bit-normalize', is_flag=True,
-              help="Apply one-bit normalization to data in each window.  Note that the default time-domain normalization "
-                   "is N(0,1), i.e. 0-mean and unit variance")
 @click.option('--location-preferences', default=None,
               type=click.Path(exists=True),
               help="A comma-separated two-columned text file containing location code preferences for "
@@ -465,7 +466,7 @@ def main(data_source1, data_source2, output_path, window_seconds, window_overlap
          stacking_interval_seconds, window_buffer_length, resample_rate, taper_length, nearest_neighbours,
          pair_min_dist, pair_max_dist, fmin, fmax, station_names1, station_names2, pairs_to_compute,
          start_time, end_time, instrument_response_inventory, instrument_response_output, water_level,
-         clip_to_2std, whitening, whitening_window_frequency, one_bit_normalize, location_preferences,
+         time_domain_norm, whitening, whitening_window_frequency, location_preferences,
          ds1_zchan, ds1_nchan, ds1_echan, ds2_zchan, ds2_nchan, ds2_echan, corr_chan, envelope_normalize,
          ensemble_stack, subset_stack, restart, dry_run, no_tracking_tag, scratch_folder):
     """
@@ -538,7 +539,7 @@ def main(data_source1, data_source2, output_path, window_seconds, window_overlap
             window_buffer_length, read_ahead_window_seconds, resample_rate, taper_length, nearest_neighbours,
             pair_min_dist, pair_max_dist, fmin, fmax, station_names1, station_names2, pairs_to_compute,
             start_time, end_time, instrument_response_inventory, instrument_response_output, water_level,
-            clip_to_2std, whitening, whitening_window_frequency, one_bit_normalize, location_preferences,
+            time_domain_norm, whitening, whitening_window_frequency, location_preferences,
             ds1_zchan, ds1_nchan, ds1_echan, ds2_zchan, ds2_nchan, ds2_echan, corr_chan, envelope_normalize,
             ensemble_stack, subset_stacker, apply_simple_stacking, restart, dry_run, no_tracking_tag,
             scratch_folder)

seismic/xcorqc/xcorqc.py

@@ -115,8 +115,8 @@ def xcorr2(tr1, tr2, sta1_inv=None, sta2_inv=None,
            instrument_response_output='vel', water_level=50.,
            window_seconds=3600, window_overlap=0.1, window_buffer_length=0,
            interval_seconds=86400, taper_length=0.05, resample_rate=None,
-           flo=None, fhi=None, clip_to_2std=False, whitening=False,
-           whitening_window_frequency=0, one_bit_normalize=False, envelope_normalize=False,
+           flo=None, fhi=None, time_domain_norm='zero_mean_one_std', whitening=False,
+           whitening_window_frequency=0, envelope_normalize=False,
            apply_simple_stacking=True, verbose=1, logger=None):
 
     # Length of window_buffer in seconds
@@ -346,34 +346,30 @@ def xcorr2(tr1, tr2, sta1_inv=None, sta2_inv=None,
                 # end if
 
                 # STEP 7: time-domain normalization
+                # Apply Rhys Hawkins-style default time domain normalization
+                if (time_domain_norm == 'zero_mean_one_std'):
+                    # 0-mean
+                    tr1_d -= np.mean(tr1_d)
+                    tr2_d -= np.mean(tr2_d)
+
+                    # unit-std
+                    tr1_d /= np.std(tr1_d)
+                    tr2_d /= np.std(tr2_d)
                 # clip to +/- 2*std
-                if clip_to_2std:
+                elif time_domain_norm == 'clip_to_two_std':
                     std_tr1 = np.std(tr1_d)
                     std_tr2 = np.std(tr2_d)
                     clip_indices_tr1 = np.fabs(tr1_d) > 2 * std_tr1
                     clip_indices_tr2 = np.fabs(tr2_d) > 2 * std_tr2
 
                     tr1_d[clip_indices_tr1] = 2 * std_tr1 * np.sign(tr1_d[clip_indices_tr1])
                     tr2_d[clip_indices_tr2] = 2 * std_tr2 * np.sign(tr2_d[clip_indices_tr2])
-                # end if
-
                 # 1-bit normalization
-                if one_bit_normalize:
+                elif time_domain_norm == 'one_bit_norm':
                     tr1_d = np.sign(tr1_d)
                     tr2_d = np.sign(tr2_d)
                 # end if
 
-                # Apply Rhys Hawkins-style default time domain normalization
-                if (clip_to_2std == 0) and (one_bit_normalize == 0):
-                    # 0-mean
-                    tr1_d -= np.mean(tr1_d)
-                    tr2_d -= np.mean(tr2_d)
-
-                    # unit-std
-                    tr1_d /= np.std(tr1_d)
-                    tr2_d /= np.std(tr2_d)
-                # end if
-
                 # STEP 8: taper
                 if adjusted_taper_length > 0:
                     tr1_d = taper(tr1_d, int(np.round(adjusted_taper_length*tr1_d.shape[0])))
@@ -520,8 +516,9 @@ def IntervalStackXCorr(refds, tempds,
                        read_ahead_window_seconds=864000, interval_seconds=86400,
                        window_seconds=3600, window_overlap=0.1, window_buffer_length=0,
                        flo=None, fhi=None,
-                       clip_to_2std=False, whitening=False, whitening_window_frequency=0,
-                       one_bit_normalize=False, envelope_normalize=False,
+                       time_domain_norm='zero_mean_one_std',
+                       whitening=False, whitening_window_frequency=0,
+                       envelope_normalize=False,
                        ensemble_stack=False,
                        subset_stacker: SubsetStacker=None,
                        apply_simple_stacking=True,
@@ -584,15 +581,14 @@ def IntervalStackXCorr(refds, tempds,
     :param flo: Lower frequency for Butterworth bandpass filter
     :type fhi: float
     :param fhi: Upper frequency for Butterworth bandpass filter
-    :type clip_to_2std: bool
-    :param clip_to_2std: Clip data in each window to +/- 2 standard deviations
+    :type time_domain_norm: str
+    :param time_domain_norm: type of time domain normalization: 'zero_mean_one_std', 'clip_to_two_std',
+                             'one_bit_norm', 'none'
     :type whitening: bool
     :param whitening: Apply spectral whitening
     :type whitening_window_frequency: float
     :param whitening_window_frequency: Window frequency (Hz) used to determine length of averaging window \
                                        for smoothing spectral amplitude
-    :type one_bit_normalize: bool
-    :param one_bit_normalize: Apply one-bit normalization to data in each window
     :type envelope_normalize: bool
     :param envelope_normalize: Envelope via Hilbert transforms and normalize
     :type ensemble_stack: bool
@@ -628,11 +624,6 @@ def IntervalStackXCorr(refds, tempds,
         if resample_rate < 2*fhi:
             raise ValueError('Resample-rate should be >= 2*fmax')
 
-    if clip_to_2std and one_bit_normalize:
-        raise ValueError('Mutually exclusive parameterization: clip_to_2std and one-bit-normalizations'
-                         'together is redundant')
-    # end if
-
     # get preferred location codes, or use the first available
     ref_loc, ref_cha = ref_loc_cha.split('.')
     temp_loc, temp_cha = temp_loc_cha.split('.')
@@ -778,10 +769,9 @@ def IntervalStackXCorr(refds, tempds,
                    resample_rate=resample_rate,
                    taper_length=taper_length,
                    flo=flo, fhi=fhi,
-                   clip_to_2std=clip_to_2std,
+                   time_domain_norm=time_domain_norm,
                    whitening=whitening,
                    whitening_window_frequency=whitening_window_frequency,
-                   one_bit_normalize=one_bit_normalize,
                    envelope_normalize=envelope_normalize,
                    apply_simple_stacking=apply_simple_stacking,
                    verbose=verbose, logger=logger)
@@ -975,9 +965,7 @@ def IntervalStackXCorr(refds, tempds,
                   'window_buffer_length': window_buffer_length,
                   'bandpass_fmin': flo if flo else -999,
                   'bandpass_fmax': fhi if fhi else -999,
-                  'clip_to_2std': int(clip_to_2std),
-                  'one_bit_normalize': int(one_bit_normalize),
-                  'zero_mean_1std_normalize': int(clip_to_2std is False and one_bit_normalize is False),
+                  'time_domain_norm': time_domain_norm,
                   'spectral_whitening': int(whitening),
                   'envelope_normalize': int(envelope_normalize),
                   'ensemble_stack': int(ensemble_stack),

tests/test_seismic/xcorqc/test_correlator.py

@@ -74,9 +74,9 @@ def test_correlator():
                 86400, 3600, 0.1,
                 0.05, 86400, 4, 0.05, -1, None, None, 0.002, 2, netsta1,
                 netsta2, None, start_time, end_time, None, 'vel',
-                50, False, True, 0.02, True, loc_pref,
+                50, 'one_bit_norm', True, 0.02, loc_pref,
                 '*Z', '*N', '*E', '*Z', '*N', '*E', 'z', False, False,
-                None, True, False, False, True, None)
+                None, True, False, False, True)
 
 
         # Read result
@@ -107,7 +107,7 @@ def test_correlator():
                 86400, 3600, 0.1,
                 0.05, 86400, 4, 0.05, -1, None, None, 0.002, 2, netsta1,
                 netsta2, None, start_time, end_time, None, 'vel',
-                50, False, True, 0.02, True, loc_pref,
+                50, 'one_bit_norm', True, 0.02, loc_pref,
                 '*Z', '*N', '*E', '*Z', '*N', '*E', 'z', False, False,
                 None, False, False, False, True, None)