Rfstats function returns empty string

Ali_Aslanli · October 17, 2024, 1:38pm

Hey fellas,

Fairly new to the seismology domain, trying to calculate receiver functions using RF package (making some experiments with local high frequency events), however, I get an empty array after executing rfstats() function. Tried to analyse the source code, but could not understand what I am missing. Maybe I am lacking some header information? Below are the stats for one trace:

       network: HP
      station: SERG
     location: 
      channel: HHE
    starttime: 2011-09-27T06:15:52.002500Z
      endtime: 2011-09-27T06:16:52.002500Z
sampling_rate: 100.0
        delta: 0.01
         npts: 6001
        calib: 1.0
      _format: SAC
  event_depth: 10.52

event_latitude: 38.3777
event_longitude: 22.0568
event_magnitude: 3.43
event_time: 2011-09-27T06:16:08.999499Z
sac: AttribDict({‘delta’: 0.01, ‘depmin’: -5680819.0, ‘depmax’: 5710203.0, ‘scale’: 1.0, ‘b’: 0.0005, ‘e’: 60.0005, ‘o’: 16.9975, ‘stla’: 38.4133, ‘stlo’: 22.0566, ‘evla’: 38.3777, ‘evlo’: 22.0568, ‘evdp’: 10.52, ‘mag’: 3.43, ‘depmen’: -6413.086, ‘nzyear’: 2011, ‘nzjday’: 270, ‘nzhour’: 6, ‘nzmin’: 15, ‘nzsec’: 52, ‘nzmsec’: 2, ‘nvhdr’: 6, ‘npts’: 6001, ‘iftype’: 1, ‘iztype’: 9, ‘leven’: 1, ‘lpspol’: 1, ‘lovrok’: 1, ‘lcalda’: 0, ‘kstnm’: ‘SERG’, ‘kcmpnm’: ‘HHE’, ‘knetwk’: ‘HP’})
station_latitude: 38.4133
station_longitude: 22.0566

trichter · October 19, 2024, 10:28pm

Hi Ali!

Please check out the documentation of the rfstats function here:

Returns: Stats object with event and station attributes, distance, back_azimuth, inclination, onset and slowness or None if epicentral distance is not in the given interval. Stream instance if stream was specified instead of stats.

In your case, station and event are too close to each other, therefore rfstats returns None.

Ali_Aslanli · October 28, 2024, 11:07am

Thanks Tom, it is more clear now.

My objective is to make use of local microseismicity, meaning for sure the events will be significantly closer to the station compared to teleseismic events.

Tried testing it while specifying the range I need but now it seems that Taupy can’t calculate the P arrival phase. Could it be that iasp91 model is not suitable for the task?

Can you suggest any workaround?

rfstats(stream, phase = ‘P’, dist_range=(0.01,80))

Exception Traceback (most recent call last)
Cell In[84], line 5

```
 1 #print(stream)*
```

 2 #print(stream[0].stats.station_latitude)*

 3 #rfstats(stream, phase='P', dist_range=(10,80))*

```
 4 #print(rfstats)*
```

----> 5 rfstats(stream, phase = ‘P’, dist_range=(0.01,80))

File ~\anaconda3\Lib\site-packages\rf\rfstream.py:683, in rfstats(obj, event, station, phase, dist_range, tt_model, pp_depth, pp_phase, model)

681 traces = *
682 for tr in stream:*
*–> 683 if rfstats(tr.stats, *kwargs) is not None:
684 traces.append(tr)*
685 stream.traces = traces*

File ~\anaconda3\Lib\site-packages\rf\rfstream.py:704, in rfstats(obj, event, station, phase, dist_range, tt_model, pp_depth, pp_phase, model)

702 arrivals = tt_model.get_travel_times(stats.event_depth, dist, (phase,))*
703 if len(arrivals) == 0:*
→ 704 raise Exception(‘TauPy does not return phase %s at distance %s’ %
705 (phase, dist))*
706 if len(arrivals) > 1:*
707 msg = ('TauPy returns more than one arrival for phase %s at '*
708 ‘distance %s → take first arrival’)*

Exception: TauPy does not return phase P at distance 0.03553634879672242

trichter · October 29, 2024, 9:30am

Not sure if rf will be working for such a local setting. With such small distances the P wave will not be going downwards from the source, but upwards. After nomenclature, the phase is p, not P. Therefore, you can try out setting phase='p' instead of phase='P' in the rftsats call.

Please check out how to correctly format code example and tracebacks. This will enhance readability. Thanks.

Ali_Aslanli · October 30, 2024, 11:12am

Your suggestion worked! Setting phase argument to ‘p’ solved the arrival problem, although by simple visual checking it is obvious that the picking is not accurate (i guess it is related to iasp91 model). But where did you check it? Is this nomenclature related to Taupy library?

Another thing if I may, while calculating RFs I have ValueError: Invalid number of source components. 0 not equal to one.

How to implement the parameters for source_components in rf function? Tried checking the source but couldn’t find info on that.

trichter · November 4, 2024, 1:07pm

No, this is standard seismological convention. I double-checked with ObsPy’s TauPy module, which is used by rf:

from obspy.taup import TauPyModel
model = TauPyModel()
print(model.get_ray_paths(10, 0.1))

17 arrivals
        p phase arrival at 2.577 seconds
        sP phase arrival at 4.342 seconds
        s phase arrival at 4.449 seconds
        PcP phase arrival at 509.544 seconds
        ScP phase arrival at 720.440 seconds
        ...

Yes, the relative time error using the iasp91 model will be greater for waves traveling a shorter distance. You can also manually adjust the metadata using manual picks, but that should usually not be necessary.

To help you with the ValuError, I would need to see some code and maybe some output from print(stream).

I have never really used the rf method on local data. I do not know if this makes sense.

Ali_Aslanli · November 5, 2024, 8:43am

Thanks Tom for all your time.

The code is pretty simple from RF’s minimal tutorial. With teleseismic data I can make it work without issues. Having a difficulty to adapt it to my local data, trying to experiment with its high freq. content, bear with me

stream = read_rf(path)
print(stream.str(extended=True))

21 Trace(s) in Stream:
HP.SERG..HHE | 2011-09-27T06:15:52.002500Z - 2011-09-27T06:16:52.002500Z | 100.0 Hz, 6001 samples
HP.SERG..HHN | 2011-09-27T06:15:52.000400Z - 2011-09-27T06:16:52.000400Z | 100.0 Hz, 6001 samples
HP.SERG..HHZ | 2011-09-27T06:15:52.004300Z - 2011-09-27T06:16:52.004300Z | 100.0 Hz, 6001 samples
HP.SERG..HHE | 2011-09-08T12:22:47.002500Z - 2011-09-08T12:23:47.002500Z | 100.0 Hz, 6001 samples
HP.SERG..HHN | 2011-09-08T12:22:47.003800Z - 2011-09-08T12:23:47.003800Z | 100.0 Hz, 6001 samples
HP.SERG..HHZ | 2011-09-08T12:22:47.000500Z - 2011-09-08T12:23:47.000500Z | 100.0 Hz, 6001 samples
HP.SERG..HHE | 2014-06-10T03:58:11.002700Z - 2014-06-10T03:59:11.002700Z | 100.0 Hz, 6001 samples
HP.SERG..HHN | 2014-06-10T03:58:11.002700Z - 2014-06-10T03:59:11.002700Z | 100.0 Hz, 6001 samples
HP.SERG..HHZ | 2014-06-10T03:58:11.004700Z - 2014-06-10T03:59:11.004700Z | 100.0 Hz, 6001 samples
HP.SERG..HHE | 2014-06-21T14:42:37.004700Z - 2014-06-21T14:43:37.004700Z | 100.0 Hz, 6001 samples
HP.SERG..HHN | 2014-06-21T14:42:37.003000Z - 2014-06-21T14:43:37.003000Z | 100.0 Hz, 6001 samples
HP.SERG..HHZ | 2014-06-21T14:42:37.002900Z - 2014-06-21T14:43:37.002900Z | 100.0 Hz, 6001 samples
HP.SERG..HHE | 2014-06-27T00:47:05.000000Z - 2014-06-27T00:48:05.000000Z | 100.0 Hz, 6001 samples
HP.SERG..HHN | 2014-06-27T00:47:05.000600Z - 2014-06-27T00:48:05.000600Z | 100.0 Hz, 6001 samples
HP.SERG..HHZ | 2014-06-27T00:47:05.000000Z - 2014-06-27T00:48:05.000000Z | 100.0 Hz, 6001 samples
HP.SERG..HHE | 2017-01-07T15:00:51.000000Z - 2017-01-07T15:01:51.000000Z | 100.0 Hz, 6001 samples
HP.SERG..HHN | 2017-01-07T15:00:51.000000Z - 2017-01-07T15:01:51.000000Z | 100.0 Hz, 6001 samples
HP.SERG..HHZ | 2017-01-07T15:00:51.000000Z - 2017-01-07T15:01:51.000000Z | 100.0 Hz, 6001 samples
HP.SERG..HNE | 2022-06-05T19:23:24.000000Z - 2022-06-05T19:24:24.000000Z | 100.0 Hz, 6001 samples
HP.SERG..HNN | 2022-06-05T19:23:24.000000Z - 2022-06-05T19:24:24.000000Z | 100.0 Hz, 6001 samples
HP.SERG..HNZ | 2022-06-05T19:23:24.000000Z - 2022-06-05T19:24:24.000000Z | 100.0 Hz, 6001 samples

rfstats(stream, phase = ‘p’, dist_range=(0.01,80))
print(stream.str(extended=True))

rf_iter = stream.copy().rf(deconvolve=‘iterative’,rotate=‘NE->RT’).moveout()

ValueError Traceback (most recent call last)
Cell In[12], line 1
----> 1 rf_iter = stream.copy().rf(deconvolve=‘iterative’,rotate=‘NE->RT’).moveout()

*File ~\anaconda3\Lib\site-packages\decorator.py:232, in decorate..fun(*args, *kw)

230 if not kwsyntax:*
231 args, kw = fix(args, kw, sig)*
*–> 232 return caller(func, *(extras + args), *kw)

*File ~\anaconda3\Lib\site-packages\rf\util.py:243, in _add_processing_info(func, *args, *kwargs)

240 arguments = [‘%s=%s’ % (k, repr(v)) if not isinstance(v, str) else*
241 “%s=‘%s’” % (k, v) for k, v in kw.items()]*
242 info = info % ‘::’.join(sorted(arguments))*
*–> 243 stream = func(*args, *kwargs)
244 try:*
245 for tr in stream:*

*File ~\anaconda3\Lib\site-packages\rf\rfstream.py:340, in RFStream.rf(self, method, filter, trim, downsample, rotate, deconvolve, source_components, *kwargs)

338 for stream3c in iter3c(self):*
339 kwargs.setdefault(‘winsrc’, method)*
→ 340 stream3c.deconvolve(method=deconvolve,
341 source_components=source_components,*
342 *kwargs)
343 # Mirrow Q/R and T component at 0s for S-receiver method for a better*
344 # comparison with P-receiver method (converted Sp wave arrives before*
345 # S wave, but converted Ps wave arrives after P wave)*
346 if method == ‘S’:*

*File ~\anaconda3\Lib\site-packages\rf\rfstream.py:258, in RFStream.deconvolve(self, *args, *kwargs)

251 def deconvolve(self, *args, *kwargs):
252 “”"*
253 Deconvolve source component of stream.*
254 *
255 All args and kwargs are passed to the function*
256 ~rf.deconvolve.deconvolve().*
257 “”"*
*–> 258 rsp = deconvolve(self, *args, *kwargs)
259 self.traces = rsp.traces*
260 return self*

*File ~\anaconda3\Lib\site-packages\decorator.py:232, in decorate..fun(*args, *kw)

230 if not kwsyntax:*
231 args, kw = fix(args, kw, sig)*
*–> 232 return caller(func, *(extras + args), *kw)

*File ~\anaconda3\Lib\site-packages\rf\util.py:243, in _add_processing_info(func, *args, *kwargs)

240 arguments = [‘%s=%s’ % (k, repr(v)) if not isinstance(v, str) else*
241 “%s=‘%s’” % (k, v) for k, v in kw.items()]*
242 info = info % ‘::’.join(sorted(arguments))*
*–> 243 stream = func(*args, *kwargs)
244 try:*
245 for tr in stream:*

*File ~\anaconda3\Lib\site-packages\rf\deconvolve.py:83, in deconvolve(stream, method, func, source_components, response_components, winsrc, *kwargs)

```
81 if len(src) != 1:*
```

82     msg = 'Invalid number of source components. %d not equal to one.'*

—> 83 raise ValueError(msg % len(src))

```
84 src = src[0]*
```

85 rsp = [tr for tr in stream if response_components is None or*

86        tr.stats.channel[-1] in response_components]*

ValueError: Invalid number of source components. 0 not equal to one.

trichter · November 5, 2024, 9:41am

Hi, this looks all good.

Not sure, what is the problem. I can have a look at the data if you like. It would be enough if you sent me data from a single event by pm.

By the way, please, use code markers ``` to post code (readability).

trichter · November 5, 2024, 12:29pm

I checked your data.

The rf() method iterates through the stream. It uses the stats.onset attribute to find out which traces belong to the same event.

The problem with your data is, that the event origin time is not the same for each three traces that belong to one event. The onsets then also differ slightly.

One straight-forward way to fix this in your script is:

for i in range(len(stream) // 3):
    for tr in stream[3*i:3*i+3]:
        tr.stats.event_time = stream[3*i].stats.event_time

Ali_Aslanli · November 5, 2024, 1:30pm

Amazing, thank you!

I’ll keep you posted once I have my results for the project. Have a nice day!

trichter · November 5, 2024, 1:49pm

You are welcome
Have a nice day, too!