Variable window length

[krisvanneste]

Claudio,

Here is a first implementation of variable window length as a function of travel time (or equivalently distance).

It contains the following modifications:

• Constrain end time of P and S windows to end time of waveform data: this may not be strictly necessary as exceeding the end time of the record will result in zero padding, but it allows limiting the noise window to the same length as the signal (P or S) window
• Take into account signal_pre_time for definition of noise window (there was a note that this should probably be done)
• Do not warn if noise window ends exactly at start of P window and updated warning message accordingly
• Set noise_pre_time to noise window length if it is not specified or zero (so it can be different for each record if variable window lengths are used)
• Added variable_win_length_factor configuration parameter and use it to calculate window length from travel time of considered phase, while keeping win_length as the minimum value
• Trim traces in trace plot to real start of noise window and real end of S window (plus some margin)
• Write variable window length to log file

There are also 2 commits that are not directly related to variable window length, but they address issues I ran into while testing:

• If the radiated energy can't be computed for a record, then plotting of the spectra results in an obscure matplotlib error (maybe because my version is too old). I made a modification to prevent a trailing newline in the label text when the Er label is empty, which avoids this error.
• When examining the reason why no Er could be computed, I found out that computation of the spectral integral is very sensitive to the fmax parameter (default None). In cases where the noise spectrum lies above the signal spectrum at the highest frequencies, noise_integral becomes quickly higher than signal_integral. I fixed this by setting fmax from spectral_snratio_fmax (available in the header of the signal spectrum) when fmax is not specified or zero. I can document this problem with a plot and manual calculations of noise_integral for different fmax values. I can include this in a next post or open a separate issue if you like.

Finally, I also considered your suggestion to pad the windows to the same length for all records. However, that is not so easy to do as it requires knowing the available lengths beforehand. I'm not sure this is really necessary, because there is already the spectral_win_length configuration parameter, which can be set to ensure padding/trimming of all windows to the same length.

I tested this variable window length for a case where I only have long-distance records (> 400 km). Using different fixed window lengths, the magnitude starts stabilizing for window lengths between 40 and 50 s, which more or less corresponds to using a variable_win_length_factor of 0.4 (for the S wave). Obviously, this would be much too long if I would include short-distance records and use a fixed window length.

Can you have a look at these modifications and let me know what you think?

Kris

[claudiodsf]

Thanks Kris!

I'll look at them next week 😉

[krisvanneste]

OK, no hurry.

[krisvanneste]

Documentation of issue with radiated energy and fmax:

This plot shows that radiated energy could not be computed for stations GR.GRA1 and GR.GRC1. The log shows that this is because noise energy is larger than signal energy: skipping spectrum. However, the noise spectrum is everywhere below the signal spectrum, except at the highest frequencies (> 9 Hz).

I did some manual computations of the spectral integral of these spectra:

• signal spectrum (fmax = None): spectral_integral = 5725.75
• signal spectrum (fmax = 9): spectral_integral = 5567.65
• noise spectrum (fmax = None): spectral_integral = 16886.83
• noise spectrum (fmax = 9): spectral_integral = 0.75

This problem can be avoided by setting fmax to spectral_snratio_fmax from the header of the signal spectrum, which is 8.77 Hz in this case.

[claudiodsf]

Kris, can you move the radiated energy problem to a separate issue?

[krisvanneste]

OK!

[claudiodsf]

I will come back to this PR, once we fix #49 😉

[krisvanneste]

Claudio,
I added a commit in which the explanation for the noise_pre_time option in configspec.conf is updated, but I can't push it to the variable_win_length branch anymore. Maybe it is easiest if I provide a patch:
variable_win_len_patch1.diff

[claudiodsf]

Kris, you're probably having problems pushing to the branch because I rebased and force-pushed it. I also slightly modified some commits.

You can try "force pull" or backup your local branch and recreate it from the remote.

Otherwise, I'll take care of adding this commit 😉

[krisvanneste]

Claudio,

Thanks, it worked with force pull.

[claudiodsf]

Thanks!

Force-pushed again after some slight modification: I prefer to use None values instead of 0 when parameters are not used (here: noise_pre_time and variable_win_length_factor)

[claudiodsf]

Going to test (and hopefully merge 😉) this afternoon!

[claudiodsf]

Hi Kris, coming back to the analysis of this PR (sorry for the delay).

I don't understand this extra - config.signal_pre_time. In fact, because of it, I get different results from the previous version, when using a fixed window length.

If I remove this extra term, then the results are identical to the previous version, when using a fixed window length.

[krisvanneste]

Claudio,
This is to ensure that the noise window ends at the start of the P window (which is equal to p_arrival_time - config.signal_pre_time) in the case that config.signal_pre_time is set to zero or None and hence defaults to the length of the signal window (but note that this is in the subsequent commit).
It may indeed give different results to before (I hope only slightly different), but in the previous version there was this note:

        # Note: maybe we should also take into account signal_pre_time here

I think the noise and signal windows are consistent this way.

Alternatively, we could leave it as before, and modify the next commit as follows:

if config.noise_pre_time == 0:
    noise_pre_time = win_length - signal_pre_time

This will give the same result when window length is variable, and should not change anything when it is fixed.

[claudiodsf]

Hi Kris, sorry again for late reaction.

Does the last commit I (force) pushed make sense to you? I tried to implement your idea (except that there should be a + and not a - in the formula 😉)

[krisvanneste]

Claudio, I'm in Germany for a meeting at the moment, so it's not so easy to look at your commit. But you are right, it should be +.

[claudiodsf]

No worries, everyone is busy! We'll find the time to merge this PR next week 😉

[krisvanneste]

Claudio, I had a look at your commit and I also tested it. Everything looks fine, so I hope it can be merged now!

[claudiodsf]

This is an old line, but it's maybe worth discussing here: it implies that t1 can be equal to t2 and that the noise window length is zero. Is this acceptable?

[krisvanneste]

True, but otherwise t2 might be < t1.
There is already a warning that the noise window ends after the start of the P window. Maybe we can add another test for the length of the noise window and warn if it is zero (or shorter than some threshold value)?

[claudiodsf]

Claudio, I had a look at your commit and I also tested it. Everything looks fine, so I hope it can be merged now!

Great! Let me rebase on current main and force-push it.

I will then ask you one last check, since I added two fixes (normally unrelated to this PR) in the current main.

[claudiodsf]

Ok, done. One last test please 🙏

[krisvanneste]

I'm not sure if I will manage. I run in 2 problems now:

• a merge conflict in ssp_inversion.py (due to SpectrumClass), which I temporarily resolved manually;
• an error which is related to the earlier fix for radiated energy calculation:

  File "C:\Users\kris\Documents\Python\cloned_repos\sourcespec\sourcespec\ssp_build_spectra.py", line 443, in _build_weight_from_inv_frequency
    snr_fmin = getattr(spec.stats, 'spectral_snratio_fmin', None)

  File "C:\Users\kris\Documents\Python\cloned_repos\sourcespec\sourcespec\spectrum.py", line 45, in __getattr__
    return self[name]

KeyError: 'spectral_snratio_fmin'

I need to check how to resolve this. Unless you know?

[claudiodsf]

The first one seems odd, since there is no SpectrumClass is not in main...

The second error seems related to the first one.

Maybe try a git pull --force? Or restart from a fresh clone?

[krisvanneste]

It's probably because I just switched to the new_spectral_class branch before merging. I will delete my local branch and try again.

[krisvanneste]

OK, I managed to get back to the variable_win_len branch and force-merge it.
The results are exactly the same and the signal windows on the trace plots look fine.

[claudiodsf]

Great! One last request: a CHANGELOG entry referring to this PR 🙏

[krisvanneste]

Only in the Config file section or in other sections as well?

[claudiodsf]

Great! Added a missing link at the end of the CHANGELOG. All is green for merging!