Use trace.remove_response() for instrument correction

[claudiodsf]

This PR modernizes the instrument correction routine remove_instr_response() which goes from using trace.simulate() and poles and zeros to trace.remove_response() and full response, as defined in the Inventory() object.

In case users provide a PAZ file, this is converted to an Inventory() object using a custom routine.

Tested and validated on all my examples.

@krisvanneste, I don't know if you use instrument correction from SourceSpec. If that's the case, I'll wait on your feedback before merging. Thanks!

[krisvanneste]

Claudio,

I will try it.
Kris

[krisvanneste]

Unfortunately, I have a problem with my obspy version. I'm still on 1.2.2, but now 1.3.0 is required.
Do you use some specific functionality in 1.3.0 that was not there in 1.2.2?
It would of course be good if I update obspy, but I tried it before and the conda updater never found a solution to sort out the dependencies. I will try to clone my environment and try forcing the update...

[claudiodsf]

The reason why I moved to ObsPy 1.3.0 is the option 'DEF' in trace.remove_response() (this line), so that I convert the trace to its nominal units (before doing the integration in ssp_build_spectra()).

A workaround could be to retrieve nominal units from inv.get_response(trace.id, trace.stats.starttime).instrument_sensitivity.input_units and set output to 'VEL', 'DISP' or 'ACC' accordingly.

I can do that! ;-)

[claudiodsf]

Ok, done!

[krisvanneste]

Thanks for the quick fix!
It runs now.
Now I get the following error message for all stations:
BE.BEBN..HHZ broadb: Unable to remove instrument response: skipping trace
However, this may be due to the way I run sourcespec. I will investigate.

[krisvanneste]

I finally fixed it. Instead of _add_paz_and_coords, I had to call _add_inventory, _add_coords and _check_instrtype. I also had to turn off clipping detection if I pre-restitute my traces before running sourcespec.
For the example I used, I obtain moment magnitudes that differ less than 0.1 (with 1 exception). The differences can best be evaluated based on the spectra.
Spectra obtained after restitution inside sourcespec:

Spectra obtained after pre-restitution to displacement:

The largest difference can be observed for station BE.UCCS. I think this is due to more narrow pre-filtering on my side.

[krisvanneste]

The example also reminds me that I should ask you if it is possible to constrain Qo to positive values...

[claudiodsf]

The example also reminds me that I should ask you if it is possible to constrain Qo to positive values...

Yes, you can either provide limits in t_star or in Qo.

From the Configuration File:

# Initial value for t_star (seconds)
t_star_0 = 0.045
# Try to invert for t_star_0.
# If the inverted t_star_0 is non-positive, then fixed t_star_0 will be used
invert_t_star_0 = False
# Allowed variability around inverted t_star_0 in the main inversion
# (expressed as a fraction of t_star_0, between 0 and 1).
# If the inverted t_star_0 is non-positive, then t_star_min_max is used
# (see below).
t_star_0_variability = 0.1

and

# t_star_min_max does not superseed t_star_0_variability
t_star_min_max = None
# optional : Qo bounds (converted into t_star bounds in the code).
# (comment out or use None to indicate no bound)
# Note: if you want to explore negative t_star values, you have to specify
# -Qo_min, Qo_min. This beacause t_star is proportional to 1/Qo.
# Example, for searching only positive t_star values:
#   Qo_min_max = 10, 1000
# If you want to search also negative t_star values:
#   Qo_min_max = -10, 10
Qo_min_max = None

I generally set t_star_min_max

[claudiodsf]

For the example I used, I obtain moment magnitudes that differ less than 0.1 (with 1 exception). The differences can best be evaluated based on the spectra.

What's your metadata format? StationXML, datalessSEED, PAZ, other?

[krisvanneste]

Thanks for the pointer!

[krisvanneste]

For the example I used, I obtain moment magnitudes that differ less than 0.1 (with 1 exception). The differences can best be evaluated based on the spectra.

What's your metadata format? StationXML, datalessSEED, PAZ, other?

It's StationXML in this case; in other cases it may also be dataless SEED.

[claudiodsf]

It's StationXML in this case; in other cases it may also be dataless SEED.

Ok, thanks.

Do you think we can validate this, or do you prefer to do more testing?

[krisvanneste]

I don't think I can test more from my side. I always work with full obspy station inventories, not with PAZ.

[claudiodsf]

Ok, thanks. So, I'm going to merge this and wait for bug reports ::😜

[krisvanneste]

OK!