pyglimer.rf#

pyglimer.rf.create#

A toolset to create RFs and RF classes

Database management and overview for the PyGLImER database.

copyright:

The PyGLImER development team (makus@gfz-potsdam.de).

license:

GNU Lesser General Public License, Version 3 (https://www.gnu.org/copyleft/lesser.html)

author:

Peter Makus (makus@gfz-potsdam.de)

Created: Friday, 12th February 2020 03:24:30 pm Last Modified: Wednesday, 19th October 2022 11:28:13 am

The file is split and has a second copyright disclaimer Some parts of this code are modified versions of the rf module by Tom Eulenfeld.

class pyglimer.rf.create.RFStream(traces=None)[source]#

Bases: Stream

Class providing the necessary functions for receiver function calculation. :param traces: list of traces, single trace or stream object To initialize a RFStream from a Stream object use >>> rfstream = RFStream(stream) To initialize a RFStream from a file use >>> rfstream = read_rf(‘test.SAC’) Format specific headers are loaded into the stats object of all traces.

bootstrap(b: int = 1000, vmodel_file: str = 'iasp91.dat') → List[ndarray][source]#

Monte Carlo Case resampling algorithm that creates b numbers of replacement resampled stacks of the moveout corrected receiver function.

The returned list can be used to compute statistical measures like standard deviation, standard error, median, etc.

Parameters:

  • b (int, optional) – Number of iterations that the algorithm is supposed to perform, defaults to 1000.

  • vmodel_file (str) – The velocity model to use for the depth migration. Defaults to ‘iasp91.dat’

Returns:

A list of the bootstraped stacks.

Return type:

List[np.ndarray]

dirty_ccp_stack(dlon: float = 1.0, z_res: float = 1.0, extent=[-180.0, 180.0, -90.0, 90.0], maxz=750, vmodel_file='iasp91.dat')[source]#

This is the simplest way of creating quick not really accurate CCP stacks.

Warning

There are some things that the user should be aware of. The bins right at the poles may or may not have very different areas compared to the other bins, which is a side effect of the iterative computation of bin area correction. I have an idea on how to fix it, but it is not really necessary right now. If you are working with stations at the pole, and you want to CCP stack, I’d strongly advise against using this function anyways. The bins are just too narrow. This really is just a ‘dirty’ CCP stack and should only be used as a first order check if thing are ok.

Parameters:

  • dlon (float, optional) – stetp size in longitude, by default 1.0

  • z_res (float, optional) – depth resolution, by default 1.0

  • extent (list, optional) – list of bounds [minlon, maxlon, minlat, maxlat], by default [-180.0, 180.0, -90.0, 90.0]

  • maxz (int, optional) – maxz, by default 750

  • vmodel_file (str, optional) – velocity model file. IASP91 1-D model used as standard since the assumption of rectangular bins is already a bit rough, by default ‘iasp91.dat’

Returns:

containing, vectors outlining the mesh illumination etc.

Return type:

tuple

property method#

Property for used rf method, ‘P’ or ‘S’

moveout(vmodel: str = 'iasp91.dat', multiple: bool = False, latb: Optional[tuple] = None, lonb: Optional[tuple] = None, taper: bool = True)[source]#

Depth migration of all receiver functions given Stream. Also calculates piercing points and adds them to RFTrace.stats.

Parameters:

  • vmodel (str) – Velocity model located in /data/vmodel. Standard options are iasp91.dat and 3D (GYPSuM).

  • multiple (bool, optional) – Appends two multiple mode RFs to the returned RFStream. False by default.

  • latb (Tuple, optional) – Tuple in Form (minlat, maxlat). To save RAM on 3D raytraycing. Will remain unused for 1D RT, defaults to None

  • lonb (Tuple, optional) – Tuple in Form (minlon, maxlon), defaults to None

  • taper (bool, optional) – If True, the last 10km of the RF will be tapered, which avoids jumps in station stacks. Should be False for CCP stacks, defaults to True.

Returns:

1D np.ndarray containing depths and an object of type depth.

Return type:

1D np.ndarray, RFStream

plot(channel: str = 'PRF', lim: Optional[list] = None, epilimits: Optional[list] = None, scalingfactor: float = 2.0, ax: Optional[Axes] = None, line: bool = True, linewidth: float = 0.25, outputdir: Optional[str] = None, title: Optional[str] = None, show: bool = True, format: str = 'pdf', color: str = 'seismic', bold: bool = False)[source]#

Creates plot of a receiver function section as a function of epicentral distance or single plot if len(RFStream)==1.

Parameters:

  • lim (list or tuple or None) – y axis time limits in seconds (if self.stats.type==time) or depth in km (if self.stats==depth) (len(list)==2). If None full traces is plotted. Default None.

  • epilimits (list or tuple or None = None,) – y axis time limits in seconds (len(list)==2). If None from 30 to 90 degrees plotted. Default None.

  • scalingfactor (float) – sets the scale for the traces. Could be automated in future functions(Something like mean distance between traces) Defaults to 2.0

  • line (bool) – plots black line of the actual RF Defaults to True

  • linewidth (float) – sets linewidth of individual traces

  • ax (matplotlib.pyplot.Axes, optional) – Can define an axes to plot the RF into. Defaults to None. If None, new figure is created.

  • outputdir (str, optional) – If set, saves a pdf of the plot to the directory. If None, plot will be shown instantly. Defaults to None.

  • clean (bool) – If True, clears out all axes and plots RF only. Defaults to False.

  • color (str, optional) – Color-scale to use. Options are ‘seismic’, ‘pyglimer’, or ‘mono’. Defaults to ‘seismic’.

  • bold (bool, optional) – Print titles and labels larger

Returns:

ax

Return type:

matplotlib.pyplot.Axes

plot_distribution(nbins=50, phase='P', outputfile=None, format='pdf', dpi=300, title: Optional[str] = None, bold: bool = False)[source]#

Plot back azimuth and rayparameter distributions.

Parameters:

  • nbins (int) – Number of bins. Default

  • v (float) – assummed surface velocity for the computation of the incidence angle. Default 5.8 km/s.

  • outputfile (str, optional) – Path to savefile. If None plot is not saved just shown. Defaults to None.

  • format (str, optional) – outputfile format

  • dpi (int, optional) – only used if file format is none vector.

  • bold (bool, optional) – Print titles and labels larger

Return type:

None

ppoint(vmodel_file='iasp91.dat', latb=None, lonb=None)[source]#

Calculates piercing points for all receiver functions in given RFStream and adds them to self.stats in form of lat, lon and depth.

Parameters:

  • vmodel_file (str, optional) – velocity model located in /data/vmodel. The default is ‘iasp91.dat’.

  • latb (tuple, optional) – Tuple in Form (minlat, maxlat). To save RAM on 3D raytraycing. Will remain unused for 1D RT, defaults to None.

  • lonb (tuple, optional.) – Tuple in Form (minlon, maxlon), defaults to None

slice2(starttime=None, endtime=None, reftime=None, keep_empty_traces=False, **kwargs)[source]#

Alternative slice method accepting relative times. See slice() and trim2().

station_stack(vmodel_file='iasp91.dat', multiple=False) → tuple[source]#

Performs a moveout correction and stacks all receiver functions in Stream. Make sure that Stream only contains RF from one station!

Parameters:

  • vmodel_file (str, optional) – 1D velocity model in data/vmodels. The default is ‘iasp91.dat’.

  • multiple (bool or str, optional) – Should RFs for multiples be calculated? If yes, the parameter has to be a string, providing the stacking mode; i.e. ‘linear’, ‘zk’ for a Zhu & Kanamori approach, or ‘pws’ for a phase-weighted stacking. False by default.

Returns:

z : Depth Vector. stack : Object containing station stack. RF_mo : Object containing depth migrated traces.

Return type:

z : 1D np.ndarray stack : RFTrace RF_mo : RFStream

trim2(starttime=None, endtime=None, reftime=None, **kwargs)[source]#

Alternative trim method accepting relative times. See trim().

Parameters:

  • starttime (UTCDateTime or float, optional) – starttime as UTC or seconds relative to reftime

  • endtime (UTCDateTime or float, optional) – endtime as UTC or seconds relative to reftime

  • reftime (UTCDateTime, optional) – reference time, can be an UTCDateTime object or a string. The string will be looked up in the stats dictionary (e.g. ‘starttime’, ‘endtime’, ‘onset’).

property type#

Property for the type of stream, ‘rf’, ‘profile’ or None

write(filename, format, **kwargs)[source]#

Save stream to file including format specific headers. See Stream.write() <obspy.core.stream.Stream.write> in ObsPy.

class pyglimer.rf.create.RFTrace(data=None, header=None, trace=None)[source]#

Bases: Trace

Class providing the Trace object for receiver function calculation.

This object is a modified version of Tom Eulenfeld’s rf project’s RFStream class. License below.

The MIT License (MIT)

Copyright (c) 2013-2019 Tom Eulenfeld

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the “Software”), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

moveout(vmodel: str, multiple: bool = False, latb: Optional[tuple] = None, lonb: Optional[tuple] = None, taper: bool = True)[source]#

Depth migration of the receiver function. Also calculates piercing points and adds them to RFTrace.stats.

Parameters:

  • vmodel (str) – Velocity model located in /data/vmodel. Standard options are iasp91.dat and 3D (GYPSuM).

  • multiple (bool , optional) – Should RFs from multiples be computed and returned? False by default.

  • latb (Tuple, optional) – Tuple in Form (minlat, maxlat). To save RAM on 3D raytraycing. Will remain unused for 1D RT, defaults to None

  • lonb (Tuple, optional) – Tuple in Form (minlon, maxlon), defaults to None

  • taper (bool, optional) – If True, the last 10km of the RF will be tapered, which avoids jumps in station stacks. Should be False for CCP stacks, defaults to True.

Returns:

1D np.ndarray containing depths and an RFTrace object of type depth.

Return type:

1D np.ndarray, RFTrace

plot(lim: Optional[list] = None, depth: Optional[ndarray] = None, ax: Optional[Axes] = None, outputdir: Optional[str] = None, format: str = 'pdf', clean: bool = False, std: Optional[ndarray] = None, flipxy: bool = False, color: str = 'seismic', bold: bool = False)[source]#

Creates plot of a single receiver function

Parameters:

  • lim (list or tuple or None) – x axis time limits in seconds or km (depth) (len(list)==2). If None full trace is plotted. Default None.

  • depth (numpy.ndarray) – 1D array of depths

  • ax (matplotlib.pyplot.Axes, optional) – Can define an axes to plot the RF into. Defaults to None. If None, new figure is created.

  • outputdir (str, optional) – If set, saves a pdf of the plot to the directory. If None, plot will be shown instantly. Defaults to None.

  • clean (bool) – If True, clears out all axes and plots RF only. Defaults to False.

  • std (np.ndarray, optional) – Only if self.type == stastack. Plots the upper and lower limit of the standard deviation in the plot. Provide the std as a numpy array (can be easily computed from the output of bootstrap())

  • flipxy (bool, optional) – Plot Depth/Time on the Y-Axis and amplitude on the x-axis. Defaults to False.

  • color (str, optional) – Color-scale to use. Options are ‘seismic’, ‘pyglimer’, or ‘mono’. Defaults to ‘seismic’.

  • bold (bool, optional) – Print titles and labels larger

Returns:

ax

Return type:

matplotlib.pyplot.Axes

ppoint(vmodel, latb=None, lonb=None)[source]#

Calculates piercing points for receiver function in and adds them to self.stats in form of lat, lon and depth.

Parameters:

  • vmodel_file (str, optional) – velocity model located in /data/vmodel. The default is ‘iasp91.dat’.

  • latb (tuple, optional) – Tuple in Form (minlat, maxlat). To save RAM on 3D raytraycing. Will remain unused for 1D RT, defaults to None.

  • lonb (tuple, optional.) – Tuple in Form (minlon, maxlon), defaults to None

write(filename, format, **kwargs)[source]#

Save current trace into a file including format specific headers. See Trace.write() <obspy.core.trace.Trace.write> in ObsPy.

pyglimer.rf.create.createRF(st_in: Stream, phase: str, pol: str = 'v', onset: Optional[UTCDateTime] = None, method: str = 'it', trim: Optional[Tuple[float, float]] = None, event=None, station=None, info: Optional[dict] = None)[source]#

Creates a receiver function with the defined method from an obspy stream.

Parameters:

  • st_in (Stream) – Stream of which components are to be deconvolved.

  • phase (str) – Either “P” or “S”.

  • pol (str, optional) – Polarisation to be deconvolved from. Only for phase = P, defaults to ‘v’

  • onset (UTCDateTime, optional) – Is used to shift function rather than shift if provided. If info is provided, it will be extracted from info file, defaults to None

  • method (str, optional) – Deconvolution method, “waterlevel”, ‘dampedf’ for constant damping level, ‘it’ for iterative time domain deconvoltuion, ‘multit’ for multitaper or ‘fqd’ for frequency dependently damped spectral division. The default is iterative time domain (‘it’).

  • trim (tuple, optional) – taper/truncate. Given as tuple (a, b) in s - left,right. The default is None.

  • event (event, optional) – Event File used to extract additional information about arrival. Provide either event and station or info dict. The default is None.

  • station ('~obspy.core.station', optional) – Dictionary containing station information, defaults to None

  • info (dict, optional) – Dictionary containing information about the waveform, used to extract additional information for header. The default is None.

Raises:

  • Exception – If deconvolution method is unknown.

  • Exception – For wrong trim input.

Returns:

RFTrace object containing receiver function.

Return type:

RFTrace

pyglimer.rf.create.obj2stats(event=None, station=None)[source]#

Map event and station object to stats with attributes.

Parameters:

  • event (Event, optional) – Event File. The default is None.

  • station (Station, optional) – Station file with attributes lat, lon, and elevation. The default is None.

Returns:

Stats object with station and event attributes.

Return type:

AttribDict

pyglimer.rf.create.read_by_station(network: str, station: str, phase: str, rfdir: str)[source]#

Convenience method to read all available receiver functions of one station and a particular phase into one single stream. Subsequently, one could for instance do a station stack by using station_stack().

Parameters:

  • network (str) – Network shortcut, two digits, e.g. “IU”.

  • station (str) – Station shortcut, three digits, e.g. “HRV”

  • phase (str) – Primary phase (“S” or “P”)

  • rfdir (str) – Parental directory, in that the RF database is located

Returns:

RFStream object containing all receiver function of station x.

Return type:

RFStream

pyglimer.rf.create.read_rf(pathname_or_url: Optional[str] = None, format: Optional[str] = None, **kwargs)[source]#

Read waveform files into RFStream object. See read() in ObsPy.

Parameters:

  • pathname_or_url (str) – Path to file.

  • format (e.g. "MSEED" or "SAC". Will attempt to determine from ending if = None. The default is None.) – str, defaults to None

Returns:

RFStream object from file.

Return type:

RFStream

pyglimer.rf.create.rfstats(phase, info=None, starttime=None, event=None, station=None, tt_model='IASP91')[source]#

Creates a stats object for a RFTrace object. Provide an info dic and starttime or event and station object. Latter will take longer since computations are redone.

Parameters:

  • phase (str) – Either “P” or “S” for main phase.

  • info (dict) – Dictionary containing information about waveform. Can be None if event and station are provided. Has to be given in combination with starttime, defaults to None.

  • starttime (UTCDateTime, optional) – Starttime of first trace in stream, used as identifier in info dict. The default is None.

  • event (Event, optional) – Event file, provide together with station file if info=None. The default is None.

  • station (Station, optional) – Station inventory with attributes lat, lon, and elevation, defaults to None

  • tt_model (str, optional) – TauPy model to calculate arrival, only used if no info file is given. The default is “IASP91”.

Raises:

Exception – For unavailable phases

Returns:

Stats file for a RFTrace object with event and station attributes, distance, back_azimuth, onset and ray parameter.

Return type:

AttribDict

pyglimer.rf.deconvolve#

Various Deconvolution approaches used for the RF technique.

copyright:

The PyGLImER development team (makus@gfz-potsdam.de).

license:

GNU Lesser General Public License, Version 3 (https://www.gnu.org/copyleft/lesser.html)

author:

Peter Makus (makus@gfz-potsdam.de)

Created: Wednesday, 16th October 2019 02:24:30 pm Last Modified: Wednesday, 20th July 2022 02:55:20 pm

pyglimer.rf.deconvolve.it(P, H, dt, shift=0, width=2.5, omega_min=0.5, it_max=200)[source]#

Iterative deconvolution after Ligorria & Ammon (1999) Essentially, P will be deconvolved from H.

Parameters:

  • P (1D np.ndarray) – The source wavelet estimation / denominator.

  • H (1D np.ndarray) – The enumerator (impulse response * source wavelet estimation).

  • dt (float) – Sampling interval [s].

  • shift (float, optional) – Time shift of main arrival [s]. The default is 0.

  • width (float, optional) – Gaussian width parameter. The default is 2.5.

  • omega_min (float, optional) – Convergence control parameter (percentage improvement per iteration). The default is 0.5.

  • it_max (int, optional) – Maximal number of iterations before the algorithm interrupts, defaults to 400

Returns:

rf : The receiver function. it : Number of iterations until algorithm converged. IR : Estimation of the medium’s impulse response.

Return type:

rf : np.ndarray it : int IR : np.ndarray

pyglimer.rf.deconvolve.multitaper(P: ndarray, D: ndarray, dt: float, tshift: int, regul: str)[source]#

Output has to be filtered (Noise appears in high-frequency)! multitaper: takes Ved-Kathrins code and changes inputs to make it run like Helffrich algorithm, minus normalization and with questionable calculation of pre-event noise (have to double check this and compare it to the regular FFT based estimate without any tapers).

Findings: when P arrival is not in the center of the window, the amplitudes are not unity at the beginning and decreasing from there on. Instead they peak at the time shift which corresponds to the middle index in the P time window.

TB = time bandwidth product (usually between 2 and 4) NT = number of tapers to use, has to be <= 2*TB-1 TB = 4; NT = 7; #choice of TB = 4, NT = 3 is supposed to be optimal t0 = -5; t1 = max(time); # This defines the beginning and end of the lag times for the receiver function

changed by KS June 2016, added coherence and variance estimate; output RF in time domain (rf), RF in freq. domain (tmpp), and variance of RF (var_rf); can be used as input for multitaper_weighting.m the input “regul” defines which type of regularization is used, regul=’fqd’ defines frequency dependent regularisation from pre-event noise, regul=’con’ defines adding a constant value (here maximum of pre-event noise) as regularisation

Parameters:

  • P (np.ndarray) – Source time function estimation.

  • D (np.ndarray) – Component containing the converted wave.

  • dt (float) – Sampling interval [s].

  • tshift (float) – Time shift of primary arrival [s].

  • regul (str) – Regularization, either ‘fqd’ for frequency dependentor ‘con’ for constant.

Raises:

Exception – For unknown input.

Returns:

  • rf (np.ndarray) – The receiver function.

  • lrf (np.ndarray) – The source-time wavelet convolved by itself.

  • var_rf (np.ndarray) – The receiver function’s variance.

  • tmpp (np.ndarray) – Receiver function without variance weighting.

pyglimer.rf.deconvolve.spectraldivision(v, u, ndt, tshift, regul, phase, test=False)[source]#

Function spectraldivision(v,u,ndt,tshift,regul) is a standard spectral division frequency domain deconvolution.

Parameters:

  • v (1D np.ndarray) – Source wavelet estimation (denominator).

  • u (1D np.ndarray) – Impulse response convolved by v (enumerator).

  • ndt (float) – Sampling interval [s].

  • tshift (float) – Time shift of primary arrival [s].

  • regul (str) – Regularization, can be chosen by the variable “regul”, this can be ‘con’, ‘wat’, or ‘fqd’ for constant damping factor, waterlevel, or frequency-dependent damping, respectively

  • phase (str) – Phase either “P” for Ps or “S” for Sp.

Raises:

Exception – for uknown regulisations

Returns:

qrf : Receiver function. lrf : Output of deoncolution of source wavelet estimation from longitudinal component.

Return type:

1D np.ndarray

pyglimer.rf.moveout#

Contains functions for moveout correction and station stacking.

copyright:

The PyGLImER development team (makus@gfz-potsdam.de).

license:

GNU Lesser General Public License, Version 3 (https://www.gnu.org/copyleft/lesser.html)

author:

Peter Makus (makus@gfz-potsdam.de)

!The file is split and has a second copyright disclaimer!

Created: Tuesday, 7th April 2020 10:08:30 am Last Modified: Monday, 30th May 2022 02:56:38 pm

class pyglimer.rf.moveout.SimpleModel(z, vp, vs, n=None)[source]#

Bases: object

Simple 1D velocity model for move out and piercing point calculation. Calculated with Earth flattening algorithm as described by P. M. Shearer.

Parameters:

  • z – depths in km

  • vp – P wave velocities at provided depths in km/s

  • vs – S wave velocities at provided depths in km/s

  • n – number of support points between provided depths

All arguments can be of type numpy.ndarray or list. taken from the RF module based on obspy by Tom Eulenfeld

pyglimer.rf.moveout.dt_table(rayp: float, fname: str, phase: str, el: float, multiple: bool, debug: bool = False)[source]#

Creates a phase delay table and calculates piercing points for a specific ray parameter, using the equation as in Rondenay 2009. For SRF: The length of the output vectors is determined by the depth, at which Sp conversion is supercritical (vertical slowness = complex).

Parameters:

  • rayp (float) – ray-parameter given in s/deg.

  • fname (string) – 1D velocity model, located in data/vmodels.

  • phase (string) – Either “S” for Sp or “P” for Ps.

  • el (float) – station elevation in m.

  • multiple (bool) – Compute conversion time for PPS and PSS.

  • debug (bool, optional) – If True, A Cartesioan rather than a spherical Earth is assumed. By default False.

Returns:

  • htab (np.array) – Vector containing conversion depths.

  • dt (np.array) – Vector containing delay times between primary arrival and converted wave.

  • delta (np.array) – Vector containing angular distance of piercing point from the station at depth z.

pyglimer.rf.moveout.dt_table_3D(rayp: float, phase: str, lat, lon, baz, el, latb, lonb, multiple: bool, test=False)[source]#

Creates a phase delay table and calculates piercing points for a specific ray parameter, using the equation as in Rondenay 2009. For SRF: The length of the output vectors is determined by the depth, at which Sp conversion is supercritical (vertical slowness = complex).

Parameters:

  • rayp (float) – ray-parameter given in s/deg.

  • fname (string) – 1D velocity model, located in data/vmodels.

  • phase (string) – Either “S” for Sp or “P” for Ps.

  • el (float) – station elevation in m.

  • latb (Tuple) – Tuple in Form (minlat, maxlat). To save RAM on 3D raytraycing. Will remain unused for 1D RT.

  • lonb (Tuple) – Tuple in Form (minlon, maxlon)

  • test (Bool) – If True, the raysum 3D model is loaded.

Returns:

  • htab (np.array) – Vector containing conversion depths.

  • dt (np.array) – Vector containing delay times between primary arrival and converted wave.

  • delta (np.array) – Vector containing Euclidian distance of piercing point from the station at depth z.

pyglimer.rf.moveout.earth_flattening(maxz: int, z: ndarray, dz: ndarray, vp: ndarray, vs: ndarray)[source]#

Creates a flat-earth approximated velocity model down to maxz as in Peter M. Shearer

Parameters:

  • maxz (int) – Maximal interpolation depth [km]. Given at beginning of file.

  • z (np.array) – Depth vector [km].

  • dz (np.array) – Layer thicknesses in km.

  • vp (np.array) – P wave velocities [km/s].

  • vs (np.array) – S wave velocities [km/s].

Returns:

  • z (np.array) – Same as input, but shortened to maxz.

  • dz (np.array) – Same as input, but shortened to maxz.

  • zf (np.array) – Depths in a Cartesian coordinate system.

  • dzf (np.array) – Layer thicknesse in a Cartesian coordinate system.

  • vpf (np.array) – P wave velocities in a Cartesian coordinate system.

  • vsf (np.array) – S wave velocities in a Cartesian coordinate system.

pyglimer.rf.moveout.load_model(fname='iasp91.dat')[source]#

Load 1D velocity model from file. The model file should have 4 columns with depth, vp, vs, n. The model file for iasp91 starts like this:

#IASP91 velocity model
#depth  vp    vs   n
  0.00  5.800 3.360 0
  0.00  5.800 3.360 0
10.00  5.800 3.360 4
Parameters:

fname (string, optional) – filename of model in data or ‘iasp91’. The default is ‘iasp91.dat’.

Returns:

Returns SimpleModel instance.

Return type:

SimpleModel

pyglimer.rf.moveout.moveout(data: ndarray, st: Stats, fname: str, latb: tuple, lonb: tuple, taper: bool, multiple: bool = False)[source]#

Depth migration for RF. Corrects the for the moveout of a converted phase. Flips time axis and polarity of SRF.

Parameters:

  • data (np.ndarray) – Receiver Function.

  • st (obspy.core.trace.Stats) – Stream stats

  • fname (str) – 1D velocity model for moveout correction. Use ‘3D’ for a 3D raytracing.

  • latb (tuple) – Tuple in Form (minlat, maxlat). To save RAM on 3D raytraycing. Will remain unused for 1D RT.

  • lonb (tuple) – Tuple in Form (minlon, maxlon)

  • taper (bool) – If True, the last 10km of the RF will be tapered, which avoids jumps in station stacks. If False, the upper 5 km will be tapered. Should be False for CCP stacks.

  • multiple (bool, optional) – Either False (don’t include multiples), ‘linear’ for linear stack, ‘PWS’ (phase weighted stack or “zk” for a zhu & kanamori approach). False by default., by default False

Returns:

  • z2 (Arraylike) – Depth vector in km.

  • RF2 (Arraylike) – Vector containing the depth migrated RF.

  • delta (Arraylike) – Vector containing Euclidian distance of piercing point from the station at depth z.

pyglimer.rf.moveout.ppoint(q_a: float, q_b: float, dz: float, p: float, phase: str)[source]#

Calculate spherical distance between piercing point and station. INPUT has to be in Cartesian/ flat Earth:

Parameters:

  • q_a (float) – Vertical P-wave slowness in s/km.

  • q_b (float) – Vertical S-wave slowness in s/km.

  • dz (float) – Vertical distance between two consecutive piercing points [km].

  • p (float) – Horizontal slowness in s/km.

  • phase (string) – Either “S” or “P”.

Returns:

x_delta – Euclidean distance between station and piercing point.

Return type:

float