diff --git a/CHANGELOG.md b/CHANGELOG.md index 28c8b751..de9b3fdd 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/sourcespec/config_files/configspec.conf b/sourcespec/config_files/configspec.conf index 3e7ad102..d4c73d86 100644 --- a/sourcespec/config_files/configspec.conf +++ b/sourcespec/config_files/configspec.conf @@ -126,8 +126,10 @@ 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) @@ -135,6 +137,10 @@ 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 diff --git a/sourcespec/ssp_plot_spectra.py b/sourcespec/ssp_plot_spectra.py index d1e6ca8a..ffc50ac5 100644 --- a/sourcespec/ssp_plot_spectra.py +++ b/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', diff --git a/sourcespec/ssp_plot_traces.py b/sourcespec/ssp_plot_traces.py index 965f9bb6..58c896b5 100644 --- a/sourcespec/ssp_plot_traces.py +++ b/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) diff --git a/sourcespec/ssp_process_traces.py b/sourcespec/ssp_process_traces.py index 76cebfdd..1c64ca33 100644 --- a/sourcespec/ssp_process_traces.py +++ b/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)