How to import my own SAC data and calculate SRF

  1. I need to understand how should I import my own data and make the figures
  2. I have a set of calculated and selected P-RF. Can I import the .eqr files and just plot the figures using RF software?
  1. Please have look here:
  1. Probably not, rf only supports waveforms which can be read by ObsPy.

Hey this is me again. I am trying to calculate P-RF of my selected dataset. I could run the python codes given in notebook 1 and got my results but when I am trying to run the iterative deconvolution , I am not getting the desired result.
I have used the following codes
Step 1.
import matplotlib.pyplot as plt
from rf import read_rf, rfstats
stream = read_rf(’/home/arpita/Documents/LADAKH/select/DRBK*’, debug_headers=True)
print(stream)
print(’\nStats:\n’, stream[0].stats)
stream[:3].plot()

step 2
rfstats(stream)
stream.filter(‘bandpass’, freqmin=0.4, freqmax=1)
stream.trim2(5, 95, ‘starttime’)
print(stream)
stream[:3].plot(type=‘relative’, reftime=stream[0].stats.onset)

step3
stream.rf()
stream.moveout()
stream.trim2(-5, 22, ‘onset’)
print(stream)
stream.select(component=‘L’).plot_rf()
stream.select(component=‘Q’).plot_rf()
plt.show()

step 4
import matplotlib.pyplot as plt
from rf import read_rf, rfstats
stream = read_rf(’/home/arpita/Documents/LADAKH/select/DRBK*’, debug_headers=True)
rfstats(stream)
stream.rf(deconvolve=‘iterative’, gauss=2.5, itmax=1000, normalize=0, response_components=None, winsrc=‘P’)
stream.plot_rf()
plt.show()

I am attaching the resultant rf that I am getting which is not matching with my desired one. P.S I have already calculated the P-RF of this set of SAC files using Liggorias iterative deconvolution and that is how I know how the final figure should look like.

Can you please assist me here to understand where I am doing the steps wrong to get my P-RF?

Hi - have you tried filtering your data after you re-import it in step 4? Unless the data you’re reading in have already been filtered elsewhere that might be a good place to start for figuring this out.

import matplotlib.pyplot as plt
from rf import read_rf, rfstats
stream = read_rf(’/home/arpita/Documents/LADAKH/select/DRBK*’, debug_headers=True)
rfstats(stream)
stream.filter(‘bandpass’, freqmin=0.4, freqmax=1)
stream.rf(deconvolve=‘iterative’, gauss=2.5, itmax=1000, normalize=0, response_components=None, winsrc=‘P’)
stream.plot_rf()
plt.show()

So I added the 5th line for filtering the data but still found no much change in the outcome.

I’m not familiar with Liggorria’s code so it’s hard to tell what might cause you to get different results, but you might check how the gaussian parameter is defined in that code versus this one. The value you input (2.5) might not mean exactly the same thing in the calculation here as it does in the other code. Also, depending on the length of the traces you’re using, it’s possible that trimming them before the deconvolution would help, though that effect is probably small since you’ve set the number of iterations fairly high at 1000.

Thanks for your response. I have tried to play with the number of iterations and gaussian value also I have trimmed the data but still not finding much change in the figure. All I am trying to get is the figures that are displayed in notebook 2 example but with my own sac data. I can obtain the same images given in notebook 1 using the commands. Can you please guide me to understand how should I change the commands given in notebook 2 to obtain results with my own sac files?

Indeed, those receiver functions look messy, with some peaks on the source component going up and some going down. In step 3 the rfs look as expected? Can you try to feed only 3 traces from one event to the rf method? By the way, you can leave away the kwargs normalize, response_components, winsrc which you only set to the default values.

hello trichter . Finally I can plot my own sac files in RF with a slight change in some codes .
I can plot the time domain waterlevel deconvolution using this . But as I am changing the ‘time’ parameter into ‘freq’ ‘iter’ or ‘multi’ some error messages are appearing. How should I make it flexible to calculate both iterative and waterlevel deconvolution?
if calculcate_RF==“PRF”:
RFstream.rf(deconvolve=‘time’,rotate=‘NE->RT’,solve_toeplitz=‘scipy’)

It is difficult to help you without more information. Can you provide a short self-contained example?

Yes I can send a sample data. Should I mail it to you?

Thanks for the waveforms, but I can only only iterate:

It is difficult to help you without more information.

For example, can you specify which error messages are appearing.