Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Response removal prevented unneccessarily for response list type response in some cases #2988

Merged
merged 6 commits into from Mar 4, 2022

Conversation

megies
Copy link
Member

@megies megies commented Mar 4, 2022

In the case of Responses containing (or well, sometimes solely built of) a Response list type stage (Frequency-Amplitude-Phase tuples), we are forcefully raising an exception if an extrapolation would have to occur, i.e. the given response list stage does not span the full range of frequencies needed for response removal.

if min_f < min_f_avail or max_f > max_f_avail:
msg = (
"Cannot calculate the response as it contains a "
"response list stage with frequencies only from "
"%.4f - %.4f Hz. You are requesting a response from "
"%.4f - %.4f Hz.")
raise ValueError(msg % (min_f_avail, max_f_avail, min_f,
max_f))

File .../lib/python3.9/site-packages/obspy/core/inventory/response.py:1385, in Response._call_eval_resp_for_frequencies(self, frequencies, output, start_stage, end_stage, hide_sensitivity_mismatch_warning)
   1379 if min_f < min_f_avail or max_f > max_f_avail:
   1380     msg = (
   1381         "Cannot calculate the response as it contains a "
   1382         "response list stage with frequencies only from "
   1383         "%.4f - %.4f Hz. You are requesting a response from "
   1384         "%.4f - %.4f Hz.")
-> 1385     raise ValueError(msg % (min_f_avail, max_f_avail, min_f,
   1386                             max_f))
   1388 amp = scipy.interpolate.InterpolatedUnivariateSpline(
   1389     f, amp, k=3)(frequencies)
   1390 phase = scipy.interpolate.InterpolatedUnivariateSpline(
   1391     f, phase, k=3)(frequencies)

ValueError: Cannot calculate the response as it contains a response list stage with frequencies only from 0.1000 - 80.4000 Hz. You are requesting a response from 0.0038 - 125.0000 Hz.

When this was implemented this was probably chosen just as the easiest and safest route possible, however it prevents response removal in some very valid scenarios. Even if the response list stage does not cover all frequencies technically needed, it may well cover the frequency range of interest (a lot of times a bandpass is applied later on anyway) and in any case it is very easy to ensure there are no problems by using pre_filt (or water_level) to suppress out-of-range frequencies.

I think we should change this exception into a warning. If we want, we could even make the warning conditional and check in the high level routines (probably just need it in Trace.remove_response) if water_level/pre_filt seem sufficient to prevent problems.

This is a (simplified) example response to illustrate the issue, green box with red borders is the frequency range defined in the response list stage:

resp_list

When interested in e.g. 1-20 Hz, there really isn't a problem why the response shouldn't be usable, a simple pre_filt=(0.1, 0.3, 50, 80) can make sure there are no problems in the deconvolution.

PR Checklist

  • Correct base branch selected? master for new features, maintenance_... for bug fixes
  • This PR is not directly related to an existing issue (which has no PR yet).
  • If the PR is making changes to documentation, docs pages can be built automatically.
    Just add the "build_docs" tag to this PR.
  • If all tests including network modules (e.g. clients.fdsn) should be tested for the PR,
    just add the "test_network" tag to this PR.
  • All tests still pass.
  • Any new features or fixed regressions are covered via new tests.
  • Any new or changed features are fully documented.
  • Significant changes have been added to CHANGELOG.txt .
  • First time contributors have added your name to CONTRIBUTORS.txt .
  • If the changes affect any plotting functions you have checked that the plots
    from all the CI builds look correct. Add the "upload_plots" tag so that plotting
    outputs are attached as artifacts.
  • Add the "ready for review" tag when you are ready for the PR to be reviewed.

@megies megies added .core issues affecting our functionality at the very core .signal issues related to our signal processing functionalities .core.inventory issues related to our Inventory functionality labels Feb 24, 2022
@trichter
Copy link
Member

I agree. I'd suggest to just extrapolate the amplitude and phase values constantly beyond the smallest and largest frequencies in the reponse list and raise a warning as you suggest. I guess, the real amplitude of the response at undefined frequencies is lower than the amplitudes given at the lowest/highest frequency in most cases. Therefore, after response removal, amplitudes will be lower than real amplitudes for the frequencies outside of the defined range and no distortions are expected. As you said with using pre-filt you are even more on the safe side.

@megies
Copy link
Member Author

megies commented Feb 25, 2022

I'd suggest to just extrapolate the amplitude and phase values constantly beyond the smallest and largest frequencies in the reponse list

Yes, that's pretty much what is done currently if you just comment out that block that raises the exception, the spline interpolation just extrapolates outside of the given range, like shown in the response plot above (that is with the "raise" statement commented out)

response list stages can be confined to quite narrow frequency bands,
e.g. when coming from industry standard calibration protocols. in real
world cases response calculation for deconvolution on time series will
for technical reasons likely have to resort on extrapolation outside of
the frequency range defined by the response list stage. while this is
problematic it is easy to counter by appropriate `pre_filt` in response
removal, which essentially zeros out unwanted parts of the signal
spectrum anyway.

currently it is not possible to do this as we forcefully raise an
exception. this commit changes that and converts it into a warning
message that gives back control to the user.
for the warning message

just warn about the actual first/last frequency defined and don't tamper
with these boundaries for the warning
not raising an exception on extrapolation anymore but showing a warning.
only saw there already was a test case after implementing the new one,
still deleting the old test as the new one is slightly faster, not going
through read_inventory() but setting up a response list stage from
scratch
that warning probably wasnt shown before because of that obscure
calculation for how much extrapolation was deemed acceptable without
raising an exception
@megies
Copy link
Member Author

megies commented Mar 4, 2022

Changed raising an exception into showing a warning.
I'd like to have this in 1.3.0, for absolutely selfish reasons ;-)

@megies megies added this to the 1.3.0 milestone Mar 4, 2022
@megies megies added the ready for review PRs that are ready to be reviewed to get marked ready to merge label Mar 4, 2022
@megies megies merged commit ce932f6 into master Mar 4, 2022
@megies megies deleted the response_list_stage branch March 4, 2022 11:33
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
.core.inventory issues related to our Inventory functionality .core issues affecting our functionality at the very core ready for review PRs that are ready to be reviewed to get marked ready to merge .signal issues related to our signal processing functionalities
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants