How to use trim and multitaper deconvolution for SRF calculation

ObsPy Related Projects Related: rf
1

Hello,
I tried to use your rf python code to compute the S receiver functions to decipher the Moho and LAB beneath my study region. I downloaded your code and installed all the necessary packages and it is working well. Moreover, I have one doubt regarding the onset time because every time it is taking P as the onset time while using the stream.trim2 after storing the S travel time in “a” variable of SAC header. I am unable to find out where I have to change to consider S time as the onset time while using the trim2 command, therefore, I request you to help me to solve this issue. And while using the multitaper deconvolution, following error is showing.

I am running the following code:

from matplotlib import pyplot as plt
from rf import read_rf 
from rf import rfstats
stream = read_rf('DATA/PCH/S/CUT/*.SAC')
rfstats(stream)
stream.filter('bandpass', freqmin=0.03, freqmax=0.5)
stream.trim2(-25, 75, 'onset')
rf_mult = stream.copy().rf(method='S', deconvolve='multitaper', rotate='ZNE->LQT').moveout(phase='Sp')

Error:

>>> stream.copy().rf(method='S', deconvolve='multitaper', rotate='ZNE->LQT', winsrc='Sp').moveout(phase='Sp')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<decorator-gen-54>", line 2, in rf
  File "/home/omsai/.local/lib/python3.8/site-packages/rf/util.py", line 243, in _add_processing_info
    stream = _func_(*args, **kwargs)
  File "/home/omsai/.local/lib/python3.8/site-packages/rf/rfstream.py", line 335, in rf
    stream3c.deconvolve(method=deconvolve,
  File "/home/omsai/.local/lib/python3.8/site-packages/rf/rfstream.py", line 253, in deconvolve
    rsp = deconvolve(self, *args, **kwargs)
  File "<decorator-gen-53>", line 2, in deconvolve
  File "/home/omsai/.local/lib/python3.8/site-packages/rf/util.py", line 243, in _add_processing_info
    stream = _func_(*args, **kwargs)
  File "/home/omsai/.local/lib/python3.8/site-packages/rf/deconvolve.py", line 123, in deconvolve
    src.trim(onset + winsrc[0], onset + winsrc[1], pad=True, fill_value=0.)
  File "/usr/local/lib/python3.8/dist-packages/obspy/core/utcdatetime.py", line 1011, in __add__
    return UTCDateTime(ns=self._ns + int(round(value * 1e9)))
TypeError: can't multiply sequence by non-int of type 'float'
2

Hi, welcome to the forum.

I will have a look.
Can you please additionally show the contents of stream[0].stats?

3

Thank you.
The contents of stream[0].stats as follows

Stats({'sampling_rate': 100.0, 'delta': 0.01, 'starttime': UTCDateTime(2000, 11, 18, 5, 41, 20, 707000), 'endtime': UTCDateTime(2000, 11, 18, 5, 47, 20, 687000), 'npts': 35999, 'calib': 1e+09, 'network': '', 'station': 'PCH', 'location': '', 'channel': 'BHE', 'sac': AttribDict({'delta': 0.0099999998, 'depmin': -342.80176, 'depmax': 540.90839, 'scale': 1e+09, 'b': 0.0, 'e': 359.97998, 'o': -1125.608, 'a': 179.983, 'stla': 10.53, 'stlo': 76.339996, 'stel': nan, 'evla': -5.5900002, 'evlo': 152.88, 'evdp': 33.0, 'dist': 8667.3984, 'az': 281.97604, 'baz': 98.068489, 'gcarc': 77.862747, 'depmen': 0.0027554086, 'nzyear': 2000, 'nzjday': 323, 'nzhour': 5, 'nzmin': 41, 'nzsec': 20, 'nzmsec': 707, 'nvhdr': 6, 'norid': 0, 'nevid': 0, 'npts': 35999, 'iftype': 1, 'iztype': 9, 'leven': 1, 'lpspol': 0, 'lcalda': 1, 'unused23': 0, 'kstnm': 'PCH', 'kuser1': 'S', 'kcmpnm': 'BHE', 'kinst': '', 'kevnm': ''}), '_format': 'SAC', 'station_latitude': 10.53, 'station_longitude': 76.339996, 'station_elevation': nan, 'event_latitude': -5.5900002, 'event_longitude': 152.88, 'event_depth': 33.0, 'event_time': UTCDateTime(2000, 11, 18, 5, 22, 35, 98968), 'onset': UTCDateTime(2000, 11, 18, 5, 34, 29, 844333), 'phase': 'P', 'distance': 77.94393806498742, 'back_azimuth': 98.04045484118696, 'inclination': 16.838115313547863, 'slowness': 5.553392462022277, 'processing': ["ObsPy 1.2.2: filter(options={'freqmin': 0.03, 'freqmax': 0.5}::type='bandpass')"]})
4

When calling rfstats you have to specify which phase to use. If you are interested in Srf you should call rfstats with the keyword phase='S'. rfstats writes the following entries: 'distance', 'back_azimuth', 'inclination', 'onset', 'slowness', 'phase' and others. In your case rfstats overwrites the onset from the SAC file with the onset calculated with taupy (P onsets because you did not specify the phase). If you want to use your own onsets, you can save these before calling rfstats and write them to the stats dictionary afterwards.

Note that Sp converted waves arrive before the S wave. Therefore, I think it is preferable to not trim too much data before the S onset.

The error occurs, because onset < starttime.

5

Thank you very much, now the error has been rectified.

6

Hello,

The program is running without any errors, however, for some waveforms, program is giving following error, can you help me to solve this issue.

rf_time = stream.copy().rf(method='S', deconvolve='time', rotate='ZNE->LQT', solve_toeplitz='scipy').moveout(phase='Sp')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<decorator-gen-54>", line 2, in rf
  File "/home/omsai/.local/lib/python3.8/site-packages/rf/util.py", line 243, in _add_processing_info
    stream = _func_(*args, **kwargs)
  File "/home/omsai/.local/lib/python3.8/site-packages/rf/rfstream.py", line 335, in rf
    stream3c.deconvolve(method=deconvolve,
  File "/home/omsai/.local/lib/python3.8/site-packages/rf/rfstream.py", line 253, in deconvolve
    rsp = deconvolve(self, *args, **kwargs)
  File "<decorator-gen-53>", line 2, in deconvolve
  File "/home/omsai/.local/lib/python3.8/site-packages/rf/util.py", line 243, in _add_processing_info
    stream = _func_(*args, **kwargs)
  File "/home/omsai/.local/lib/python3.8/site-packages/rf/deconvolve.py", line 83, in deconvolve
    raise ValueError(msg % len(src))
ValueError: Invalid number of source components. 0 not equal to one.
7

Can you check the stream object, e.g. with print(stream)?