How to get ML from waveform without complete station/inventory data? noob

I hope this is the right forum to ask for help about instrument response and using waveforms for magnitude estimation.
The problem description might be long-winded.

Background: one of the local meteorological agencies in charge of earthquake data used several different local magnitude, ML, equations that were not calibrated locally. A researcher derived a local calibration for ML, but the agency didn’t adjust previous records. There are multiple communication issues across the board, so I am on my own (with some data files) at this point. I started using ObsPy just because the tutorials are helpful, convenient, and I don’t want to spend so much time on formats.

My situation: Would like to estimate MLs (and if possible Mw) from several separate acceleration and velocity mSeed waveforms. For ML, I inferred, and confirmed with a friend in another country in a similar department, that I would need the equivalent maximum displacement amplitude from a simulated Wood-Anderson seismomograph. My friend said I would need to remove the instrument response from my waveforms, then add instrument response of a simulated Wood-Anderson type instrument, then pick the maximum amplitude from the waveform… in nm. Use that amplitude in the locally calibrated ML equation (along with the agency’s listed hypocenter information).

I am running into perhaps several problems. One is the station data and it’s relationship to instrument response. I was provided with “RESP” station files for acceleration and velocity instruments. Example station data is in a text file named “RESP.KS.USN2…HGE”:

q330hrs(6084) + ES-T-A(7400)

B050F03 Station: USN2
B050F16 Network: KS
B052F03 Location: ??
B052F04 Channel: HGE
B052F22 Start date: 2016,151,00:00:00
B052F23 End date: 3000,001,23:59:59

B053F03 Transfer function type: A
B053F04 Stage sequence number: 1
B053F05 Response in units lookup: M/S**2 - Velocity in Meters Per Second
B053F06 Response out units lookup: V - Volts
B053F07 A0 normalization factor: +2.459570e+13
B053F08 Normalization frequency: +1.00000e+00
B053F09 Number of zeroes: 0
B053F14 Number of poles: 4

Complex zeroes:

i real imag real_error imag_error

Complex poles:

i real imag real_error imag_error

B053F15-18 0 -9.810000e+02 +1.009000e+03 +0.00000e+00 +0.00000e+00
B053F15-18 1 -9.810000e+02 -1.009000e+03 +0.00000e+00 +0.00000e+00
B053F15-18 2 -3.290000e+03 +1.263000e+03 +0.00000e+00 +0.00000e+00
B053F15-18 3 -3.290000e+03 -1.263000e+03 +0.00000e+00 +0.00000e+00

SENSOR

B058F03 Stage sequence number: 1
B058F04 Sensitivity: +4.070438e+00
B058F05 Frequency of sensitivity: +1.00000e+00
B058F06 Number of calibrations: 0

RECORDER

B058F03 Stage sequence number: 2
B058F04 Sensitivity: +4.194300e+05
B058F05 Frequency of sensitivity: +1.00000e+00
B058F06 Number of calibrations: 0

SENSOR*RECORDER

B058F03 Stage sequence number: 0
B058F04 Sensitivity: +1.707264e+06
B058F05 Frequency of sensitivity: +1.00000e+00
B058F06 Number of calibrations: 0

=========================================================================================

I noticed there might be some missing or mislabled information, but the worker at the agency identified and assured us these files are related to the acceleration aspect of the station recordings. The worker informed us that he has no idea what instrument types are installed (make/model). Another worker told us to multiply the count in the mSeed file by the inverse of the product of the sensitivities for the SENSOR and RECORDER (e.g. 1/+1.707264e+06) to get acceleration in m/s**2. I do not see “GAIN” anywhere in these RESP/text files.
I also tried to remove instrument response following the ObsPy tutorial at Seismometer Correction/Simulation — ObsPy 1.4.0 documentation, I ended up getting “Exception: Failed to autodetect sampling rate of channel from response stages. Please manually specify parameter sampling_rate” when I tried to inv.plot_response(1).

After this, I am lost. Which function(s) do I use to transform the acceleration waveform to a Wood-Anderson equivalent? Do I have to work in velocity? I noticed some specific functions, i.e. “estimate_wood_anderson_amplitude()”, “estimate_wood_anderson_amplitude_using_response()”, “estimate_magnitude()”, but am not sure how I should do the inputs give the RESP file shown above. Additionally, there was some online code:

# Note Wood anderson sensitivity is 2080 as per Uhrhammer & Collins 1990
PAZ_WA={'poles': [-6.283 + 4.7124j, -6.283 - 4.7124j],
            'zeros': [0 + 0j], 'gain': 1.0, 'sensitivity': 2080}
from obspy.signal import seisSim
# De-trend data
trace.detrend('simple')
# Simulate Wood Anderson
if PAZ:
    trace.data=seisSim(trace.data, trace.stats.sampling_rate, paz_remove=PAZ,\
                       paz_simulate=PAZ_WA, water_level=water_level,\
                       remove_sensitivity=True)
elif seedresp:
    trace.data=seisSim(trace.data, trace.stats.sampling_rate, paz_remove=None,\
                       paz_simulate=PAZ_WA, water_level=water_level,\
                       seedresp=seedresp)
else:
    UserWarning('No response given to remove, will just simulate WA')
    trace.data=seisSim(trace.data, trace.stats.samplng_rate, paz_remove=None,\
                        water_level=water_level)
return trace

=========================================================================================
but is it applicable to my situation? the gains, sensitivity, poles, and zeroes? the seisSim() function?
Any suggestions on the Mw front?

Any help is appreciated!
TYIA

Bump. Need 20 characters.

You might find the code here helpful for your conversion and picking. You might need to adapt to your use case. The docs for that function are here.

Converting from resp to stationxml: I don’t have a good way of doing this. If anyone else knows a good way of doing this sing out! Otherwise pyrocko has a conversion tool, but I haven’t tried it.

We can read RESP :smile:

read_inventory("/path/to/RESP.AZ.DHL..BS1").write("...", format="STATIONXML")

Doh! I don’t know why I didn’t know/forgot that! Too easy. No need to convert either then. You should just be able to read as @megies said and then use that.

Sorry for the (added) confusion and thank you @megies for correcting me.

1 Like

Also, rather than EQcorrscan’s code that I linked earlier, you might find it helpful to go through the obspy advanced exercise that runs you through the steps and provides solutions as well.

I don’t have any advice for Mw computation.