Merge pull request #48 from krisvanneste/variable_win_len

Variable window length

CHANGELOG.md

@@ -73,6 +73,8 @@ previous versions. You will need to upgrade your old database manually or using
 - For weights computed from spectral S/N ratio (noise weighting), set to zero
   all the weights below 20% of the maximum weight, so that these weakly
   constrained parts of the spectrum are ignored in the inversion
+- Possibility of using variable signal window lengths for each station
+  as a function of the travel time of the P or S wave (see [#48])
 
 ### Inversion
 
@@ -153,6 +155,8 @@ previous versions. You will need to upgrade your old database manually or using
 - New parameter `plot_map_api_key` to provide a Stadia Maps
   api key for Stamen Terrain basemap
 - New option for the parameter `plot_coastline_resolution`: `no_coastline`
+- New config parameter `variable_win_length_factor` to specify window
+  length as a fraction of the travel time of the P/S wave (see [#48])
 
 ### Bugfixes
 
@@ -706,4 +710,5 @@ Initial Python port.
 [#40]: https://github.com/SeismicSource/sourcespec/issues/40
 [#43]: https://github.com/SeismicSource/sourcespec/issues/43
 [#44]: https://github.com/SeismicSource/sourcespec/issues/44
+[#48]: https://github.com/SeismicSource/sourcespec/issues/48
 [#49]: https://github.com/SeismicSource/sourcespec/issues/49

sourcespec/config_files/configspec.conf

@@ -126,15 +126,21 @@ NLL_time_dir = string(default=None)
 p_arrival_tolerance = float(min=0, default=4.0)
 s_arrival_tolerance = float(min=0, default=4.0)
 
-# Start time (in seconds) of the noise window, respect to the P arrival time
-noise_pre_time = float(default=6.0)
+# Start time (in seconds) of the noise window, respect to the P window
+# If None, it will be set to the length of the signal (P or S) window plus
+# the value of "signal_pre_time" (see below)
+noise_pre_time = float(min=0.01, default=6.0)
 
 # Start time (in seconds) of the signal window, respect to the P or S arrival
 # times (see "wave_type" below)
 signal_pre_time = float(default=1.0)
 
 # Length (in seconds) for both noise and signal windows
 win_length = float(min=0, default=5.0)
+# Variable window length factor (fraction of travel time):
+#   win_length = max(win_length, variable_win_length_factor * travel_time)
+# Set to None to disable variable window length.
+variable_win_length_factor = float(min=0.01, default=None)
 # -------- TIME WINDOW PARAMETERS
 
 

sourcespec/ssp_plot_spectra.py

@@ -539,9 +539,10 @@ def _params_text(spec, ax, color, path_effects, stack_plots):
         sep = ' '
     params_text = (
         f'{Mo_text} {Mw_text}\n'
-        f'{fc_text}{sep}{t_star_text}\n'
-        f'{Er_text}'
+        f'{fc_text}{sep}{t_star_text}'
     )
+    if Er_text:
+        params_text += f'\n{Er_text}'
     ax.text(
         0.05, params_text_ypos, params_text,
         horizontalalignment='left',

sourcespec/ssp_plot_traces.py

@@ -333,8 +333,8 @@ def _set_ylim(axes):
 
 def _trim_traces(config, st):
     for trace in st:
-        t1 = trace.stats.arrivals['P'][1] - config.noise_pre_time
-        t2 = trace.stats.arrivals['S'][1] + 3 * config.win_length
+        t1 = trace.stats.arrivals['N1'][1]
+        t2 = trace.stats.arrivals['S2'][1] + 2 * config.win_length
         trace.trim(starttime=t1, endtime=t2)
     # compute time offset for correctly aligning traces when plotting
     min_starttime = min(tr.stats.starttime for tr in st)

sourcespec/ssp_process_traces.py

@@ -292,22 +292,44 @@ def _define_signal_and_noise_windows(config, trace):
             'Using (Ts-Tp)/2 instead')
     t1 = s_arrival_time - s_pre_time
     t1 = max(trace.stats.starttime, t1)
-    t2 = t1 + config.win_length
+    tt_s = trace.stats.travel_times['S']
+    # manage case where variable_win_length_factor is None:
+    variable_win_length_factor = config.variable_win_length_factor or 0
+    win_length_s = max(config.win_length, variable_win_length_factor * tt_s)
+    t2 = t1 + win_length_s
+    t2 = min(trace.stats.endtime, t2)
+    win_length_s = t2 - t1
     trace.stats.arrivals['S1'] = ('S1', t1)
     trace.stats.arrivals['S2'] = ('S2', t2)
     # Signal window for spectral analysis (P phase)
     t1 = p_arrival_time - config.signal_pre_time
     t1 = max(trace.stats.starttime, t1)
-    t2 = t1 + min(config.win_length, s_minus_p + s_pre_time)
+    tt_p = trace.stats.travel_times['P']
+    win_length_p = max(config.win_length, variable_win_length_factor * tt_p)
+    t2 = t1 + min(win_length_p, s_minus_p + s_pre_time)
+    t2 = min(trace.stats.endtime, t2)
+    win_length_p = t2 - t1
     trace.stats.arrivals['P1'] = ('P1', t1)
     trace.stats.arrivals['P2'] = ('P2', t2)
     # Noise window for spectral analysis
-    t1 = max(trace.stats.starttime, p_arrival_time - config.noise_pre_time)
-    t2 = t1 + config.win_length
-    if t2 >= p_arrival_time:
-        logger.warning(f'{trace.id}: noise window ends after P-wave arrival')
-        # Note: maybe we should also take into account signal_pre_time here
-        t2 = p_arrival_time
+    if config.wave_type[0] == 'P':
+        win_length = win_length_p
+    elif config.wave_type[0] == 'S':
+        win_length = win_length_s
+    if config.variable_win_length_factor:
+        logger.info(f'{trace.id}: window length {win_length:.3f} seconds')
+    if config.noise_pre_time is None:
+        noise_pre_time = win_length + config.signal_pre_time
+    else:
+        noise_pre_time = config.noise_pre_time
+    t1 = max(
+        trace.stats.starttime,
+        p_arrival_time - noise_pre_time
+    )
+    t2 = t1 + win_length
+    if t2 > (p_arrival_time - config.signal_pre_time):
+        logger.warning(f'{trace.id}: noise window ends after P-window start')
+        t2 = p_arrival_time - config.signal_pre_time
         t1 = min(t1, t2)
     trace.stats.arrivals['N1'] = ('N1', t1)
     trace.stats.arrivals['N2'] = ('N2', t2)