paudetseis / rfpy Goto Github PK

While using rfpy_ccp PBA.pkl --start=92.74,11.66 --end=92.75,11.65 --snr=10 --snrh=5 --no-outlier --load --prep --prestack --ccp --gccp --linear --figure , I am getting this error.

Traceback (most recent call last):
File "c:\users\uplc\anaconda3\envs\rfpy\lib\runpy.py", line 194, in _run_module_as_main
return _run_code(code, main_globals, None,
File "c:\users\uplc\anaconda3\envs\rfpy\lib\runpy.py", line 87, in run_code
exec(code, run_globals)
File "C:\Users\UPLC\anaconda3\envs\rfpy\Scripts\rfpy_ccp.exe_main.py", line 7, in
File "c:\users\uplc\anaconda3\envs\rfpy\lib\site-packages\rfpy\scripts\rfpy_ccp.py", line 623, in main
ccpimage.save(prep_file)
File "c:\users\uplc\anaconda3\envs\rfpy\lib\site-packages\rfpy\ccp.py", line 637, in save
if not ".pkl" in title:
TypeError: argument of type 'WindowsPath' is not iterable

However, the message also shows this before throwing error:

While running the command rfpy_recalc --align=ZRT --fmax=0.5 --gfilt=1. AC07.pkl for gaussian filter with rfpy_recalc script, it throws following error:

Traceback (most recent call last):
File "c:\users\uplc\anaconda3\envs\rfpy\lib\runpy.py", line 194, in _run_module_as_main
return _run_code(code, main_globals, None,
File "c:\users\uplc\anaconda3\envs\rfpy\lib\runpy.py", line 87, in run_code
exec(code, run_globals)
File "C:\Users\UPLC\anaconda3\envs\rfpy\Scripts\rfpy_recalc.exe_main.py", line 7, in
File "c:\users\uplc\anaconda3\envs\rfpy\lib\site-packages\rfpy\scripts\rfpy_recalc.py", line 377, in main
rfdata.deconvolve(
File "c:\users\uplc\anaconda3\envs\rfpy\lib\site-packages\rfpy\rfdata.py", line 851, in deconvolve
rfL, rfQ, rfT = _decon(trL, trQ, trT, trNl, nn, method)
File "c:\users\uplc\anaconda3\envs\rfpy\lib\site-packages\rfpy\rfdata.py", line 752, in _decon
gauss = _gauss_filt(dt, npad, gfilt)
File "c:\users\uplc\anaconda3\envs\rfpy\lib\site-packages\rfpy\rfdata.py", line 669, in _gauss_filt
gauss[nft21:] = np.flip(gauss[1:nft21-1])
ValueError: could not broadcast input array from shape (724,) into shape (725,)

I am getting following error. How to resolve?
Please help.
I tried changing directory path too but didn,t work.

FileNotFoundError Traceback (most recent call last)
Input In [2], in <cell line: 2>()
1 from rfpy import RFData
----> 2 rfdata = RFData('demo')

File ~\anaconda3\envs\rfpy\lib\site-packages\rfpy-0.1.1-py3.9.egg\rfpy\rfdata.py:182, in RFData.init(self, sta)
179 import os
180 import pickle
181 sta = pickle.load(
--> 182 open(os.path.join(
183 os.path.dirname(file),
184 "examples/data", "MMPY.pkl"), 'rb'))['NY.MMPY']
186 # Attributes from parameters
187 self.sta = sta

FileNotFoundError: [Errno 2] No such file or directory: 'C:\Users\UPLC\anaconda3\envs\rfpy\lib\site-packages\rfpy-0.1.1-py3.9.egg\rfpy\examples/data\MMPY.pkl'

Hi!
Thanks for sharing this repository.
It would be great if you can share example jupyter notebook on demo data covering all the analyses.

Thanks

Currently, RFData.deconvolve._decon uses the entire time window of the parent for deconvolution. It might be beneficial to select a shorter wavelet e.g. based on the signal envelope or a user-defined time window.

I suggest to include a new private function to RFData.deconvolve, e.g.

def _wavelet(method='complete', envelope_threshold=0.1, time=5):

    """
    Select wavelet from the parent function for deconvolution using method.
    method: (str)
        'complete' use complete, tapered parent signal (current implementation)
        'envelope' use only the part of the parent signal where envelope > envelope_threshold*max(envelope)
        'time' use only the first window seconds after signal`
    envelope_threshold: (float) fraction of the envelope that defines the wavelet (for method='envelope')
    time: (float) window (seconds) that defines the wavelet (for method='time')
    """

In rfpy.Meta.__init__(), snr, snrh, cc, ttime, slow and inc get initialized as None. However, when making quality checks such as:

if meta.snr > snr_threshold:
    continue

undetermined values raise a TypeError

It would be better to initialize them as np.nan. Then, the above tests would evaluate as False, as expected. The internal tests for availability of the values could be re-written as:

if self.meta.snr:  # False if None (for backward compatibility)
    if np.isfinite(self.meta.snr):  # False if nan (None throws error)
        print("SNR already calculated - continuing")
        return

I have already implemented this, but did not push it into the future branch, as it may break backward compatibility for scripts that rely on the implicit evaluation of if self.meta.snr as False.

Hi Pascal,
np.float was deprecated in NumPy ≥1.20 and this will result in the weight not being calculated correctly when stacking receiver functions into (baz or slow) bins. This should also be one of the reasons why the rfpy_plot doesn't work.

for tr in rf_tmp[0]:
IndexError: list index out of range

I'm just curious why the denominator in the Wiener deconvolution doesn't include preevent noise spectra from the S components as in Audet (2010).

Hii expert, I have calculated moho depth below each station.Now i want to do a image plot of the moho depth on the basemap,how can it be done can you please suggest.

Hi Pascal!

I was following the given tutorial to compute PRF, I was able to get results until Back-azimuth panel.
My command was the following:

rfpy_plot.py --no-outlier --bp=0.05,0.5 --nbaz=36 --normalize --trange=-2.,30. MMPY.pkl

The error is in the plot script.

Traceback (most recent call last):
File "/home/qian/anaconda3/envs/rfpy/bin/rfpy_plot.py", line 588, in
main()
File "/home/qianl/anaconda3/envs/rfpy/bin/rfpy_plot.py", line 523, in main
for tr in rf_tmp[0]:
IndexError: list index out of range

Do you know what should be the problem? Let me know if you require more information

I have other question.
I plan to apply PRF for two networks, one with short-period instruments (1 s corner frequency) and
other with broad -band sensors (30 s corner frequency).

Do you think that it is best to process them separately? For example band-passing the short period seismograms from 0.5 Hz to let's say 1.5 Hz and do the same for the broad-band seismograms ?
Or is it better to process them using the same band-pass filter for both networks? lets say a adjust to a more or less common frequency between 0.05 Hz to 1.5 Hz

Thanks in advance!

As a new user to receiver functions, can you please provide/suggest some code for the computation of the rf of a offline station from scratch. I have only raw data recorded in a station saved in a directory. How can i compute rf using the same data sets. Moreover I have stla,stlo,elevation information with me.

• Accept either StDb or obspy.inventory object (or seamlessly switch between them using StDb tools)
• Include QC thresholding using hierarchical clustering of correlation matrix
• Convert core functions (H-k, CCP) to FORTRAN to improve speed

when i try to plot after rfpy_calc, this error shown:
Traceback (most recent call last):
File "C:\Users\andre\anaconda3\envs\rfpy\lib\runpy.py", line 197, in _run_module_as_main
return _run_code(code, main_globals, None,
File "C:\Users\andre\anaconda3\envs\rfpy\lib\runpy.py", line 87, in run_code
exec(code, run_globals)
File "C:\Users\andre\anaconda3\envs\rfpy\Scripts\rfpy_plot.exe_main.py", line 7, in
File "C:\Users\andre\anaconda3\envs\rfpy\lib\site-packages\rfpy\scripts\rfpy_plot.py", line 523, in main
for tr in rf_tmp[0]:
IndexError: list index out of range

Hi!
I didnt find the command in 'rfpy_hk' script to convert the harmonic decomposition from time to depth domain.
Can u tell how it is done in Audet, P. (2015), Layered crustal anisotropy around the San Andreas Fault near Parkfield, California, J. Geophys. Res. Solid Earth, 120, 3527–3543, doi:10.1002/2014JB011821 ?

utils.download_data() useses obspy.client.get_waveforms(..., attach_response=False). With this option, traces can't be restituted to true ground displacement. What is the rationale behind this choice? I think RfPy should make an attempt to retrieve the instrument response and remove it, before returning the traces.

The final receiver functions are fftshifted just before they're returned, but I can only reproduce the receiver function calculation demo if I remove the fftshift.

I am not sure this is a true issue, but I am just looking for clarification about:

1. How should I interpret the time axis after the fftshift?
2. Why should we apply the fftshift to the time-domain receiver function? Looking at the vertical component receiver function (rfp), we should expect a delta function (to 0th order). I think the delta function should be at elapsed time t=0, but it plots at the midpoint of the time axis due to the fftshift.

Thank you, and please let me know if I can clarify my questions!

I have this error while calculating RF in terminal.
Please help.

Thanks.

When the HkStack and Harmonics classes are initialized with 'demo', they load receiver functions from demo_streams.pkl. The receiver functions in demo_streams.pkl were created (dc0eb66) before the _add_rfstats nested function of rfdata's to_stream method included taxis (and other stats like phase, cc, etc.) (0fe006e). The following lines throw an AttributeError because taxis does not exist.

I think demo_streams.pkl just needs to be recreated to include these stats.