'slowness' error while using stream.ppoints()

Hello! I’m trying to determine peircing points using rf. I’m new to python and I’m not quite sure I understood the doc correctly. Can anyone help me out?

Here’s my code:

from rf import read_rf
stream = read_rf('ev1.SAC')
stream.ppoints(50,'P')

Here’s the error:

C:\Users\Maureen\anaconda3\lib\site-packages\obspy\core\trace.py:220: UserWarning: Calibration factor set to 0.0!
  warnings.warn(msg, UserWarning)
Traceback (most recent call last):

  File "C:\Users\Maureen\Documents\M1 Stage\Python\PyTauP_test.py", line 11, in <module>
    stream.ppoints(50,'P')

  File "C:\Users\Maureen\anaconda3\lib\site-packages\rf\rfstream.py", line 391, in ppoints
    model.ppoint(tr.stats, pp_depth, phase=pp_phase)

  File "C:\Users\Maureen\anaconda3\lib\site-packages\rf\simple_model.py", line 229, in ppoint
    dr = self.ppoint_distance(depth, stats['slowness'], phase=phase)

  File "C:\Users\Maureen\anaconda3\lib\site-packages\obspy\core\trace.py", line 235, in __getitem__
    return super(Stats, self).__getitem__(key, default)

  File "C:\Users\Maureen\anaconda3\lib\site-packages\obspy\core\util\attribdict.py", line 74, in __getitem__
    return self.__dict__[name]

KeyError: 'slowness'

Hello!

To calculate the piercing point the slowness of the ray has to be known.
Can you upload your file ev1.SAC or post the output of print(stream[0].stats)?

Thanks for your answer. Ok, I thought that the module was using as input iasp97 (or something like that ) as a velocity model :thinking: I’ve uploaded the file. But it seems that I have a problem because now I get an error “no module named rf” when I run the previous code… I’m the right environment but something’s wrong and I can’t understand why…When I run conda activate rfenv, sets path for LIB and CPATH or does this mean that I missed a step in the installation process? I’m using spyder as python interpreter.

(Attachment ev1.SAC is missing)

I couldn’t upload the SAC file “unauthorized”…

Yes the iasp model is used by default, but to calculate piercing points you need to know something about the incoming ray.

I think when you want to use rf in spyder, you need to install rf into the same environment spyder is installed in. You can check the installed packages with conda list.

I couldn’t upload the SAC file “unauthorized”…

What a pitty, maybe we can improve on that, see Upload attachments to post.

Edit: Uploading waveform files should work now.

Also, please have a look at the rfstats function. It can calculate the slowness from station and event coordinates.

It’s weird because it was working fine with spyder until now… Anyway is the slowness vp/vs? I used TauP once and I don’t remember entering any slowness parameter…

This is the list from conda list - backcall 0.2.0 pyh9f0ad1d_0 conda-forge - backports 1.0 py_2 conda-forge - backports.functools_lru_cache 1.6.1 py_0 conda-forge - blas 1.0 mkl - boost-cpp 1.72.0 h89d28cc_2 conda-forge - brotlipy 0.7.0 py38h1e8a9f7_1000 conda-forge - bzip2 1.0.8 hfa6e2cd_2 conda-forge - ca-certificates 2020.6.20 hecda079_0 conda-forge - cairo 1.16.0 hc79e3cf_1005 conda-forge - cartopy 0.18.0 py38hde6099d_1 conda-forge - certifi 2020.6.20 py38h32f6830_0 conda-forge - cffi 1.14.1 py38hba49e27_0 conda-forge - cfitsio 3.470 hbbe6aef_6 conda-forge - chardet 3.0.4 py38h32f6830_1006 conda-forge - clangdev 5.0.0 flang_3 [flang] conda-forge - cloudpickle 1.5.0 py_0 conda-forge - colorama 0.4.3 py_0 conda-forge - cryptography 3.0 py38hba49e27_0 conda-forge - curl 7.71.1 h4b64cdc_3 conda-forge - cycler 0.10.0 py_2 conda-forge - decorator 4.4.2 py_0 conda-forge - expat 2.2.9 he025d50_2 conda-forge - flang 5.0.0 he025d50_20180525 conda-forge - flang_win-64 5.0.0 20180525 conda-forge - fortran-compiler 1.1.1 he67e8c7_0 conda-forge - freetype 2.10.2 hd328e21_0 conda-forge - freexl 1.0.5 hd288d7e_1002 conda-forge - future 0.18.2 py38h32f6830_1 conda-forge - gdal 3.1.2 py38h5f218aa_1 conda-forge - geographiclib 1.50 py_0 conda-forge - geos 3.8.1 he025d50_0 conda-forge - geotiff 1.6.0 h09e6dc1_1 conda-forge - gettext 0.19.8.1 hb01d8f6_1002 conda-forge - glib 2.65.0 he4de6d7_0 conda-forge - h5py 2.10.0 nompi_py38ha551744_104 conda-forge - hdf4 4.2.13 hf8e6fe8_1003 conda-forge - hdf5 1.10.6 nompi_ha405e13_100 conda-forge - icc_rt 2019.0.0 h0cc432a_1 - icu 67.1 h33f27b4_0 conda-forge - idna 2.10 pyh9f0ad1d_0 conda-forge - intel-openmp 2019.4 245 - ipykernel 5.3.4 py38h5ca1d4c_0 conda-forge - ipython 7.16.1 py38h23f93f0_0 conda-forge - ipython_genutils 0.2.0 py_1 conda-forge - jedi 0.17.2 py38h32f6830_0 conda-forge - jpeg 9d he774522_0 conda-forge - jupyter_client 6.1.6 py_0 conda-forge - jupyter_core 4.6.3 py38h32f6830_1 conda-forge - kealib 1.4.13 h3b59ab9_1 conda-forge - kiwisolver 1.2.0 py38heaebd3c_0 conda-forge - krb5 1.17.1 hc04afaa_1 conda-forge - libblas 3.8.0 14_mkl conda-forge - libcblas 3.8.0 14_mkl conda-forge - libcurl 7.71.1 h4b64cdc_3 conda-forge - libffi 3.2.1 h6538335_1007 conda-forge - libflang 5.0.0 h6538335_20180525 conda-forge - libgdal 3.1.2 h586bd2e_1 conda-forge - libiconv 1.15 hfa6e2cd_1006 conda-forge - libkml 1.3.0 h7e985d0_1011 conda-forge - liblapack 3.8.0 14_mkl conda-forge - libnetcdf 4.7.4 nompi_h256d12c_105 conda-forge - libpng 1.6.37 hfe6a214_1 conda-forge - libpq 12.3 hd9aa61d_0 conda-forge - libsodium 1.0.17 h2fa13f4_0 conda-forge - libspatialite 4.3.0a h9ed2f40_1039 conda-forge - libssh2 1.9.0 hb06d900_4 conda-forge - libtiff 4.1.0 h885aae3_6 conda-forge - libwebp-base 1.1.0 hfa6e2cd_3 conda-forge - libxml2 2.9.10 h5d81f1c_2 conda-forge - libxslt 1.1.33 h579f668_1 conda-forge - lxml 4.5.2 py38he3d0fc9_0 conda-forge - lz4-c 1.9.2 h62dcd97_1 conda-forge - m2w64-expat 2.1.1 2 - m2w64-gcc-libgfortran 5.3.0 6 - m2w64-gcc-libs 5.3.0 7 - m2w64-gcc-libs-core 5.3.0 7 - m2w64-gettext 0.19.7 2 - m2w64-gmp 6.1.0 2 - m2w64-libiconv 1.14 6 - m2w64-libwinpthread-git 5.0.0.4634.697f757 2 - m2w64-xz 5.2.2 2 - matplotlib-base 3.3.0 py38hfb9ee82_1 conda-forge - mkl 2019.4 245 - mkl-service 2.3.0 py38hfa6e2cd_0 conda-forge - msys2-conda-epoch 20160418 1 - numpy 1.19.1 py38h72c728b_0 conda-forge - obspy 1.2.2 py38h4b8e87e_0 conda-forge - olefile 0.46 py_0 conda-forge - openjpeg 2.3.1 h57dd2e7_3 conda-forge - openmp 5.0.0 vc14_0 conda-forge - openssl 1.1.1g he774522_0 conda-forge - owslib 0.20.0 py_0 conda-forge - parso 0.7.1 pyh9f0ad1d_0 conda-forge - pcre 8.44 h6538335_0 conda-forge - pickleshare 0.7.5 py38h32f6830_1001 conda-forge - pillow 7.2.0 py38h7011068_1 conda-forge - pip 20.1.1 py_1 conda-forge - pixman 0.38.0 hfa6e2cd_1003 conda-forge - poppler 0.89.0 h5d62644_1 conda-forge - poppler-data 0.4.9 1 conda-forge - postgresql 12.3 he14cc48_0 conda-forge - proj 7.1.0 h7d85306_1 conda-forge - prompt-toolkit 3.0.5 py_1 conda-forge - pycparser 2.20 pyh9f0ad1d_2 conda-forge - pyepsg 0.4.0 py_0 conda-forge - pygments 2.6.1 py_0 conda-forge - pyopenssl 19.1.0 py_1 conda-forge - pyparsing 2.4.7 pyh9f0ad1d_0 conda-forge - pyproj 2.6.1.post1 py38h5a885fd_1 conda-forge - pyreadline 2.1 py38_1001 conda-forge - pyshp 2.1.0 py_0 conda-forge - pysocks 1.7.1 py38h32f6830_1 conda-forge - python 3.8.5 h5fd99cc_1_cpython conda-forge - python-dateutil 2.8.1 py_0 conda-forge - python_abi 3.8 1_cp38 conda-forge - pytz 2020.1 pyh9f0ad1d_0 conda-forge - pywin32 227 py38hfa6e2cd_0 conda-forge - pyyaml 5.3.1 py38h9de7a3e_0 conda-forge - pyzmq 19.0.1 py38h77b9d75_0 conda-forge - requests 2.24.0 pyh9f0ad1d_0 conda-forge - scipy 1.5.0 py38h9439919_0 - setuptools 49.2.0 py38h32f6830_0 conda-forge - shapely 1.7.0 py38hbf43935_3 conda-forge - six 1.15.0 pyh9f0ad1d_0 conda-forge - spyder-kernels 1.9.3 py38h32f6830_0 conda-forge - sqlalchemy 1.3.18 py38h1e8a9f7_0 conda-forge - sqlite 3.32.3 he774522_1 conda-forge - tbb 2020.1 he980bc4_0 conda-forge - tiledb 2.0.6 h0b90766_0 conda-forge - tk 8.6.10 hfa6e2cd_0 conda-forge - tornado 6.0.4 py38hfa6e2cd_0 conda-forge - tqdm 4.48.0 pyh9f0ad1d_0 conda-forge - traitlets 4.3.3 py38h32f6830_1 conda-forge - urllib3 1.25.10 py_0 conda-forge - vc 14.1 h869be7e_1 conda-forge - vs2015_runtime 14.16.27012 h30e32a0_2 conda-forge - wcwidth 0.2.5 pyh9f0ad1d_1 conda-forge - wheel 0.34.2 py_1 conda-forge - win_inet_pton 1.1.0 py38_0 conda-forge - wincertstore 0.2 py38_1003 conda-forge - xerces-c 3.2.3 ha925a31_1 conda-forge - xz 5.2.5 h62dcd97_1 conda-forge - yaml 0.2.5 he774522_0 conda-forge - zeromq 4.3.2 ha925a31_3 conda-forge - zlib 1.2.11 h2fa13f4_1006 conda-forge - zstd 1.4.5 h1f3a1b7_2 conda-forge

The spyder kernels are installed… I’ve sent you the sac file via e-mail.

Hi,

I do not see rf in the list of installed packages.

Everything does nicely work with your provided file. You just have to calculate the slowness with the rfstats function:

from rf import read_rf, rfstats
s = read_rf('ev1.SAC')
rfstats(s)
s.ppoints(50)

Please take some time and find out what this slowness is about. It’s not the vp/vs ratio, but the ray parameter. If you have no resource at hand, in the rf doc you will find a link to a dissertation with some information about receiver functions.

OK, you were right, I reinstalled and now it works! Thank you so much for your help😀! Just a quick question, is it possible to make loop in order to have all the events done at the same time?

Sure, that is possible. You can also load everything into one stream. The code block above works also for a stream of multiple traces.

Sorry to bother you again but I tried to put a pathname to read all the sac files but I get an error “permission denied”. I tried to find a way to join my sac files into one through ObsPy but I can’t find it either… I’m a little lost…

Something like read_rf('ev1.SAC') + read_rf('ev2.SAC') or using a glob expression with read_rf('*.SAC') should work.

OK thank you very much for your patience and assistance :blush:

Hi! it’s me again: I tried to import all my files with glob using two methods but none of them seems to be working. One returns only 9 piercing points out of
77 and with the wrong coordinates and the other returns an empty array… Here’s the codes for both methods:
method A:

e=[]
from rf import read_rf, rfstats
import os, glob
folder_path='Users\Maureen\Documents\M1 Stage\Python\LN14'
for filename in glob.glob(os.path.join(folder_path,'*.SAC')):
    with open(filename, 'r') as f:
        s =read_rf()
        rfstats(s)
        e.append(s.ppoints(50,'P'))
print(e)

return : []

method B:

from rf import read_rf, rfstats
import os, glob
for file in glob.glob('*.SAC'):
    s= read_rf()
    rfstats(s)
    s.ppoints(50,'P')
    print(s.ppoints(50,'P'))

return: [[-20.82703869 -69.64815651]
[-20.82703869 -69.64815651]
[-20.82703869 -69.64815651]
[-21.26771754 -69.34484486]
[-21.26771754 -69.34484486]
[-21.26771754 -69.34484486]
[-20.77253233 -69.6304343 ]
[-20.77253233 -69.6304343 ]
[-20.77253233 -69.6304343 ]]

I know those coordinates are wrong because the station is located in Tanzania and those points are not… I’m leaning towards a problem using glob: the program works for one event… Do you know if I can use an ascii file or a .txt file or use pandas? I built a txt file with pandas that contains backazimut and the distance between the event and the station, and I seem to remember that TauP was capable of finding the piercing points with those parameters…

@Sarenrenegade you can see here for how to properly format code blocks (and other things). :slight_smile:
Edited your above post. :wink:

You are not using the variable f in method A and file in method B.

Even when I use file and filename it doesn’t work… It’s problem with the pathname… The quickest solution was to work in the directory😊
Thanks for your help :grinning: