Hi there,

Working through the Chile tutorial for a lab I'm putting together and hit a snag here:

pstream = read_rf(profilefile)
pstream.trim2(-5, 20, 'onset')
pstream.select(channel='??Q').normalize().plot_profile(scale=1.5, top='hist')
plt.gcf().set_size_inches(15, 10)
plt.show()

File ~/.local/lib/python3.11/site-packages/rf/imaging.py:318, in plot_profile(profile, fname, figsize, dpi, scale, fillcolors, trim, top, moveout_model)
    316 pd = model.calculate_delay_times(phase=phase, slowness=slowness)
    317 ax2 = ax.twinx()
--> 318 ax.get_shared_y_axes().join(ax, ax2)
    319 dkm = 50
    320 if profile[0].stats.endtime - profile[0].stats.onset > 50:

AttributeError: 'GrouperView' object has no attribute 'join'

matplotlib 3.8.2 so I suppose this is probably the issue.

Replacing line 318 with ax.sharey(ax2) solves it for me.

Sorry Sir, may I know the references you have used from the time deconvolution method? Thank you.

Applying a RF calculation using the ZRT coordinate system gives wrong result.

For example, the following code computes the RF for a ZRT rotation on the first set of the example stream loaded with rf_read():

from rf import read_rf, rfstats

stream = read_rf()
stream = stream[:3]
rfstats(stream)

stream.detrend()

stream.filter("bandpass", freqmin=0.01, freqmax=1)

stream.rf(rotate="NE->RT", method='P')
stream.select(component="R").plot(type='relative',
                                  reftime=stream[0].stats.onset)

We expect something like what Ammon (1991) shows in this Figure:

But instead, we get the following RF in the radial direction:

......../home/sph/seismo/rf/rf/simple_model.py:83: RuntimeWarning: invalid value encountered in sqrt qp = np.sqrt(self.vp ** (-2) - slowness ** 2)
It happened when I test moveout in test_simple_modle.py. How can I fix this warning? What is the limit of slowness?
Appreciate your job very much!

The iterator iter_event_data does not gracefully handle when the Origin of an Event is missing an expected field. Namely, if any of {latitude, longitude, depth, magnitude, or time} are missing, then the call to stats.update(obj2stats(event=event, station=station)) in function rf_stats can raise an unhandled exception. This is difficult for the user to handle in a recoverable way, since the iterator will be left before the user code has a chance to handle it. It should probably be handled (with a warning message) in rf_stats, and then return an empty object. This would allow the iterator iter_event_data to continue yielding events.

To reproduce the issue, use a valid Event with valid Origin, and del the depth field of the Origin. Then pass such data to iter_event_data to iterate over.

A real example of a typical, valid Origin:

Origin
resource_id: ResourceIdentifier(id="smi:ISC/origid=609833474")
time: UTCDateTime(2017, 9, 11, 17, 35, 9, 620000) [uncertainty=2.24]
longitude: 142.3373
latitude: 23.9529
depth: 35000.0
quality: OriginQuality(used_phase_count=166, standard_error=1.1, azimuthal_gap=19.0, minimum_distance=3.132)
origin_uncertainty: OriginUncertainty(min_horizontal_uncertainty=9.9, max_horizontal_uncertainty=12.2, azimuth_max_horizontal_uncertainty=258.0, preferred_description='uncertainty ellipse')
creation_info: CreationInfo(author='NEIC')

A real example of a bad Origin, missing the depth field, taken from Event id "smi:ISC/evid=611767244":

Origin
resource_id: ResourceIdentifier(id="smi:ISC/origid=609672372")
time: UTCDateTime(2018, 3, 7, 14, 46, 12, 200000)
longitude: 57.65
latitude: 28.03
quality: OriginQuality(azimuthal_gap=27.0)
creation_info: CreationInfo(author='CSEM')

As a workaround, users could remove bad events from a Catalog before passing it to iter_event_data.

Sorry Sir, may I know the references you have used from the time deconvolution method? Thank you.

Hi all,
at the moment I'm trying to calculate iterative RFs by using the PREM velocity model.
For some distances I get 'empty' receiver functions, it is happening when doing: rfstats(stream, phase=phase, dist_range=dist_range, tt_model=tt_model). My assumption is that there is a problem with the dist_range and the phase. I set P phase with a distance range of 30-95°, Pdiff for 95-105°, PP for 105-140° and PKP for 140-165°, which worked perfectly fine for the RF calculation using the iasp91 model.
Now using the prem model, I get problems with for example distances of 134° and 139°, 134° needs to be associated with a PKP, but it is not working to associate 139° with the PKP phase. I checked the arrivals using the taupy model, where I get PKP arrivals and I also checked the ray paths using cake. I also tried the phases P, PP, PKP and Pdiff (just to try, when it is working), but the rfstream is always getting empty after using rfstats.
So, do you have any idea how to solve this problem? Did I miss anything? How can I see, which distance range rf wants me to use for which phase, if it is not conform with the taupy model?

Any help is very appreciated! Many thanks! :)

When a stream is passed to rf.imaging.plot_rf(), the code assumes that all traces in the stream have the same time points, and the times array is generated from the first trace and assumed to be valid for all traces.

However it is possible occasionally for a trace to have different times if it has truncated data that makes the data array smaller than the time window of the other streams. If you try to plot such data, which should in principle be plottable, the function throws an error plotting the truncated trace, because the times points array is longer than the data array for that trace.

It would be better if the times array were computed per trace when plotting individual traces.

I'm not sure how the computing and plotting of the stacked trace should be handled in this situation. Possibly exclude the truncated trace from the stack, or else skip the stack plotting altogether.

There are links of two example notebooks in the README.md, but I can't find them in the documentation.

It would be better to list them in the documentation page, since not all users will read README in this repository.

openjournals/joss-reviews#1808

I was wondering if rf is offering some tools to generate synthetic Receiver Functions from a 1D velocity model. If not, is there any tool for doing it? Are we planning to include this kind of functionality in rf?

https://nbviewer.jupyter.org/github/trichter/notebooks/blob/master/receiver_function_profile_chile.ipynb

Issues:

1. Need to create a "data" directory before running the script
2. I got following errors, possibly an issue of matplotlib or obspy?

Traceback (most recent call last):
  File "/Users/seisman/.anaconda/envs/rfenv/lib/python3.7/site-packages/matplotlib/colors.py", line 109, in _sanitize_extrema
    ret = ex.item()
AttributeError: 'str' object has no attribute 'item'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "test.py", line 26, in <module>
    inventory.plot(label=False)
  File "/Users/seisman/.anaconda/envs/rfenv/lib/python3.7/site-packages/obspy/core/inventory/inventory.py", line 749, in plot
    title=None, show=False, fig=fig, **kwargs)
  File "/Users/seisman/.anaconda/envs/rfenv/lib/python3.7/site-packages/obspy/imaging/maps.py", line 768, in plot_map
    return plot_cartopy(*args, **kwargs)
  File "/Users/seisman/.anaconda/envs/rfenv/lib/python3.7/site-packages/obspy/imaging/maps.py", line 553, in plot_cartopy
    scal_map = ScalarMappable(norm=Normalize(min_color, max_color),
  File "/Users/seisman/.anaconda/envs/rfenv/lib/python3.7/site-packages/matplotlib/colors.py", line 909, in __init__
    self.vmin = _sanitize_extrema(vmin)
  File "/Users/seisman/.anaconda/envs/rfenv/lib/python3.7/site-packages/matplotlib/colors.py", line 111, in _sanitize_extrema
    ret = float(ex)
ValueError: could not convert string to float: '#b15928'

openjournals/joss-reviews#1808

On platforms where the from toeplitz import sto_sl import fails due to Toeplitz library not being available, consider reverting to using scipy since rf is already dependent on scipy. See
https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve_toeplitz.html

Both RTZ and LQT coordinate systems are commonly used in receiver function studies. When I read this section, I thought this package only supports LQT coordinate system, then I realized that's not true after reading the API of the rf function.

It would be better to clarify that this package supports both.

openjournals/joss-reviews#1808

Sorry I am new in python and never deal with Station XML format
If I do not have any Station XML format
Only have sacpz file like this

ZEROS 6
-90.0000  0.0000
-160.7000  0.0000
-3108.0000  0.0000
POLES 7
-0.0385  0.0366
-0.0385  -0.0366
-178.0000  0.0000
-135.0000  160.0000
-135.0000  -160.0000
-671.0000  1154.0000
-671.0000  -1154.0000
CONSTANT 1.480864e+14

Is it possible to use this rf repository ？Or how to rewrite conf.json ?
Thank you a lot to develop this useful rf project !

I have a proposed change that I have found very helpful for analyzing the position of arrival multiples in RF plots when there are many plots on one axes. I have added pale dotted vertical lines by turning on grid for the xaxis in RF plots (optional, off by default).

I shall generate a pull request off a branch for the change to be reviewed.

I'm using macOS. I followed the install instructions (create an environment, and replace gfortran_linux-64 with gfortran_osx-64), and had trouble installing the toeplitz package.

Please see the attached logfile errors.log.

openjournals/joss-reviews#1808

Dear all,
I got this error in one specific case: UnboundLocalError: local variable 'sumsq' referenced before assignment
on row 541 of deconvole.py : lsfit[c] = 100 - (sumsq * 100)

it occurs when the error is small (in my case) and doe not go inside the loop on row 513 (deconvolce.py):

shouldn't be lsfit[c] = 100 - (sumsq_i * 100) ?? in this case, at least, the UnboundLocalError goes away and there is a Runtime warning (there is a divide by zero).

I am sorry to get a bit unclear, but I am a newby on python. I hope this helps

Function signature of function rf.imaging.plot_rf() will be updated to pass kwargs on to fig.savefig(). This will allow more control over the image file parameters such as dpi resolution.

I am running this code on python notebook 3.10.
I am getting following two errors,

AttributeError: can't set attribute 'collections'.

ImportError: cannot import name 'Iterable' from 'collections' (/home/seismology11/anaconda3/envs/myenv/lib/python3.10/collections/init.py)

Please tell me, how can i resolve this issue.

In function rf.deconvolve.deconvolve(), when the response data is gathered for passing on to the underlying deconvolution function (e.g. deconvt or deconvf), it is put into a list of numpy.array, as at line 113:
rsp_data = [tr.data for tr in rsp]

In my use case of a custom deconvolution function this works fine, however there are output metadata of the deconvolution process that ideally I would like to store with the traces. If the deconvolution function were passed RFTrace objects rather than numpy.arrays then I could add the metadata to the RFTrace.stats container.

Hi! I quite new to python and I'm having a hard time understanding how to use ppoint to find piercing points for SAC file with stream because I get an error "key error slowness"...

Can you help me out ? thanks

Hi Tom!

I had some trouble to convert an obspy stream to rfstream. The "Q" file written from RFStream cannot be read by seismic handler, when the COMMENT is filled with a json dictionary.

Thus, as workaround I converted all Q-Files written from RFStream to be readable with seismic handler by replacing sh.COMMENT with "" (empty string) and writing station_lat, station_lon into DCVREG,DCVINCI header variables:

for tr in st:
    tr.stats.sh.DCVREG   = tr.stats.station_latitude
    tr.stats.sh.DCVINCI  = tr.stats.station_longitude
    tr.stats.sh.COMMENT  = ""
st.write(SH_RFFILE, format="Q")

Those Q-FIles can be read by seismic handler.
However reading those Q-Files again with obspy and trying to convert them back to RFStream-objects fails, since the COMMENT is empty.

Error message:

json.decoder.JSONDecodeError: Expecting value: line 1 column 1 (char 0)

I could fix this by replacing

            elif format == 'sh' and head_format == 'COMMENT':
                st.update(json.loads(value))
                continue

with

            elif format == 'sh' and head_format == 'COMMENT':
                if len(value)>0:
                    st.update(json.loads(value))
                continue

in _read_format_specific_header in rfstream.py

Great job, but moveout for Sp doesn't work properly:

In [177]: new_t = np.interp(t, t0, t1, left=0, right=0)
In [178]: new_t
Out[178]: array([ 0., 0., 0., ..., 0., 0., 0.])

For Sp you have to change right to None:
new_t = np.interp(t, t0, t1, left=0, right=None)

If the array is empty, the function should simply return 0. Instead, after np.searchsorted returns 0, if tries to index into the empty array.

I was wondering if this is appropriate for S RFs since the time shift is muting the signal of interest which arrives before onset. I could get meaningful S RFs by setting the time shift to 0.

Hello,
I installed rf as it is given in the documentation. Any help will be of great help.

Hi!
I am getting this error while attempting to run the codes.
Please help

How to solve this onset problem?

I noticed the commented code here that preserves the event ID when iterating an event catalog. It would be helpful if this event ID was passed through from the original catalog data, if available.

Hi there,

In the equation for the Gaussian low pass filter
gauss = np.exp(np.maximum(-(0.5 * 2 * pi * freq / gauss) ** 2, -700.) -
1j * tshift * 2 * pi * freq)

I believe the 0.5 should be outside of the brackets: -0.5*( 2 * pi * freq / gauss) ** 2
per https://en.wikipedia.org/wiki/Gaussian_filter#Definition equation 4.

Also a related question, what's behind the choice of the Gauss parameter being 2pi times the standard deviation of the resulting Gaussian?

Cheers

For example if n is int and n == 3 and len(s) = 2, then the logic will continue through to len(s) in n, but n being int is not iterable, so this raises TypeError: TypeError: argument of type 'int' is not iterable

The logic on this line needs to check for iterability before applying in to n.

During iter_event_data, RF catches all exception types, meaning you can't even stop execution using keyboard interrupt (CTRL-C) since this is of exception type KeyboardInterrupt.

Consider narrowing the except catching statements to except Exception:, as this is a base class for all the exceptions I think would be intended to be caught here. See Python exception documentation.

This would then allow CTRL-C to interrupt the execution when the user wants it to.

Hello
I want to calculate receiver function by using rf function, but have problem.

When I execute code like your tutorial example ('rf_iter = stream.copy().rf(deconvolve='iterative', rotate='NE->RT', gauss=1.0') with my broadband and geophone data, I receive error message 'IndexError: index 2001 is out of bounds for axis 0 with size 2001'. The number of '2001' may be npts.

But, When I practice that function with your tutorial example, there are no error.
I already checked 'stats' by comparing the headers between tutorial example data and my data.

I have searched why I receive that error message using google.
I think the 'Numpy' has problem...

I'm python beginner, so don't know how can i treat this problem...

please help me...=(
I want to calculate receiver function....

Thank you!!

don't forget to add the recipe you found in trichter/toeplitz#4
to the travisCI tests here ;-)