_images/PyRaysum_logo.png

Synthetic seismograms are computed from three classes: Model, Geometry, and Control using prs.run(), and are stored in an instance of class Result.

Classes

pyraysum defines the following classes:

Model

class pyraysum.prs.Model(thickn, rho, vp, vs=None, flag=1, ani=None, trend=None, plunge=None, strike=None, dip=None, vpvs=1.73, maxlay=15)[source]

Model of the subsurface seismic velocity structure.

Note

This object holds the infromation of the .mod file in classic Raysum

Parameters:

  • thickn (array_like) – Thickness of layers (m)

  • rho (float or array_like) – Density (kg/m^3)

  • vp (float or array_like) – P-wave velocity (m/s)

  • vs (float or array_like, optional) – S-wave velocity (m/s) If None, computed from vpvs

  • flag (array_like of str, optional) – 1 for isotropic, 0 for anisotropic

  • ani (float or array_like, optional) – Anisotropy (percent)

  • trend (float or array_like, optional) – Trend of symmetry axis (degree)

  • plunge (float or array_like, optional) – Plunge of symmetry axis (degree)

  • strike (float or array_like, optional) – azimuth of interface in RHR (degree clockwise from north)

  • dip (float or array_like, optional) – dip of interface in RHR (degree down from horizontal)

  • vpvs (float or array_like, optional) – P-to-S velocity ratio. Defaults to 1.73. Ignored if vs is set.

  • maxlay (int) – Maximum number of layers defined in params.h

Warning

When setting vpvs, vs is adjusted to satisfy vs = vp / vpvs.

The following attributes are set upon initialization, when setting a layer property via Model.__setitem__() (i.e. model[0] or model[0, "thickn"]), and when execuding Model.update() after setting a single element of the above model attributes. f prefixes indicate attributes used for interaction with fraysum.run_bare() and fraysum.run_full(). Set these directly for best performance.

nlay

Number of layers

parameters

Convenience attribute that collects the f-attributes in the order expected by fraysum.run_bare() and fraysum.run_full()

fthickn

Thickness of layers (m)

frho

Density (kg/m^3)

fvp

P-wave velocity (m/s)

fvs

S-wave velocity (m/s)

fflag

Flag indicating isotropy (1) or anisotropy (0) of layer material

fani

Anisotropy (percent)

ftrend

Trend of symmetry axis (radians)

fplunge

Plunge of symmetry axis (radians)

fstrike

azimuth of interface in RHR (radians)

fdip

dip of interface in RHR (radians)

Example

>>> from pyraysum import Model
>>> model = Model([10000, 0], [3000, 4500], [6000, 8000], [3500, 4600])
>>> print(model)
#  thickn     rho      vp      vs  flag aniso   trend plunge strike   dip
  10000.0  3000.0  6000.0  3500.0    1    0.0     0.0    0.0    0.0   0.0
      0.0  4500.0  8000.0  4600.0    1    0.0     0.0    0.0    0.0   0.0
>>> model[0]["thickn"]
10000.0
>>> model.thickn[0]
10000.0
>>> model.thickn[0] = 15000  # model.fthickn has not yet changed
>>> model.update()  # Now everything works as expected
>>> model.fthickn[0]
15000.0
>>> model[1, "vp"]
8000.0
>>> model[1, "vp"] = 7400.  # Does not require model.update()
>>> model[1, "vp"]
7400.0
>>> model[1] = {"vs": 4200., "rho": 4000.}
>>> model[1]["vs"]
4200.0
>>> model[1, "rho"]
4000.0
>>> model[1] = {"thickn": 5000, "vp": 8000., "vpvs": 2.}
>>> model[1]["vs"]
4000.0
>>> model += [5000, 3600, 8000, 4000]  # thickn, rho, vp, vs
>>> print(model)
#  thickn     rho      vp      vs  flag aniso   trend plunge strike   dip
  15000.0  3000.0  6000.0  3500.0    1    0.0     0.0    0.0    0.0   0.0
   5000.0  4000.0  8000.0  4000.0    1    0.0     0.0    0.0    0.0   0.0
   5000.0  3600.0  8000.0  4000.0    1    0.0     0.0    0.0    0.0   0.0
>>> model += {"thickn": 0, "rho": 3800, "vp": 8500., "dip": 20, "strike": 90}
>>> print(model)
#  thickn     rho      vp      vs  flag aniso   trend plunge strike   dip
  15000.0  3000.0  6000.0  3500.0    1    0.0     0.0    0.0    0.0   0.0
   5000.0  4000.0  8000.0  4000.0    1    0.0     0.0    0.0    0.0   0.0
   5000.0  3600.0  8000.0  4000.0    1    0.0     0.0    0.0    0.0   0.0
      0.0  3800.0  8500.0  4913.3    1    0.0     0.0    0.0   90.0  20.0
update(change='vpvs')[source]

Update all attributes after one of them was changed.

Parameters:

change (str) –

Decide which of vp, vs or vpvs should depend on the other two:

  • ’vpvs’: Calculate vpvs from vpvs = vp / vs

  • ’vp’: Calculate vp from vp = vs * vpvs

  • ’vs’: Calculate vs from vs = vp / vpvs

Returns:

None

copy()[source]

Return a copy of the model

change(command, verbose=True)[source]

Change layer properties using a command string.

Parameters:

  • command (str) – An arbitrary number of command substrings separated by ‘;’

  • verbose (bool) – Print changed parameters to screen

Returns:

List of changes applied of the form: (attribute, layer, new value)

Return type:

list of 3*tuple

Note

In the command argument, each substring has the form:

KEY LAYER SIGN VAL;

where

KEY [t|vp|vs|psp|pss|s|d|a|tr|pl] is the attribute to change

  • t: thickness (km)

  • vp: P wave velocity (km/s)

  • vs: S wave velocity (km/s)

  • psp: P to S wave velocity ratio with fixed vs (changing vp)

  • pss: P to S wave velocity ratio with fixed vp (changing vs)

  • s: strike (degree)

  • d: dip (degree)

  • a: anisotropy (%)

  • tr: trend of the anisotropy axis (degree)

  • pl: plunge ot the anisotropy axis (degree)

LAYER (int) is the index of the layer

SIGN [=|+|-] is to set / increase / decrease the attribute

VAL (float) is the value to set / increase / decrease

Hint

For example, Model.change('t0+10;psp0-0.2;d1+5;s1=45') does:

  1. Increase the thickness of the first layer by 10 km

  2. Decrease Vp/Vs of the of the first layer by 0.2, holding Vs fixed

  3. Increase the dip of the second layer by 5 degree

  4. Set the strike of the second layer to 45 degree

split_layer(n)[source]

Split layer n into two with half the thickness each, but otherwise identical parameters.

Parameters:

n – (int) Index of the layer to split

remove_layer(n)[source]

Remove layer n

Parameters:

n (int) – Index of the layer to remove

average_layers(top, bottom)[source]

Combine layers between top and bottom indices into one with summed thicknesses and averaged vp, vs, and rho.

Parameters:

  • top (int) – Index before top-most layer to include in combination

  • bottom (int) – Index after bottom-most layer to include in combination

Raises:

  • IndexError – if bottom is less or equal to top

  • ValueError – if any layer is anisotropic

  • ValueError – if any layer has a difference in strike or dip

save(fname='sample.mod', comment='', hint=False, version='prs')[source]

Alias for write()

write(fname='sample.mod', comment='', hint=False, version='prs')[source]

Write seismic velocity model to disk as Raysum ASCII model file

Parameters:

  • fname (str) – Name of the output file (including extension)

  • comment (str) – String to write into file header

  • hint (bool) – Include usage comment to model file

  • version ("prs" or "raysum") – Use PyRaysum or Raysum file format

plot(zmax=75.0)[source]

Plot model as stair case, layers and labeled interfaces, and show it

Parameters:

zmax (float) – Maximum depth of model to plot (km)

plot_profile(zmax=75.0, ax=None)[source]

Plot model as stair case and show it

Parameters:

  • zmax (float) – Maximum depth of model to plot (km)

  • ax (plt.axis) – Axis handle for plotting. If None, show the plot

Returns:

Axis handle for plotting

Return type:

plt.axis

plot_layers(zmax=75.0, ax=None)[source]

Plot model as horizontal layers and show it

Parameters:

  • zmax (float) – Maximum depth of model to plot (km)

  • ax (plt.axis) – Axis handle for plotting. If None, show the plot

Returns:

Axis handle for plotting

Return type:

plt.axis

plot_interfaces(zmax=75, info=['vp', 'vs', 'vpvs', 'rho'], ax=None)[source]

Plot model as labeled interfaces with possibly dipping layers

Parameters:

  • zmax (float) – Maximum depth of model to plot (km)

  • info (list of str) – Which Model attributes to list for each layer

  • ax (plt.axis) – Axis handle for plotting. If None, show the plot

Returns:

Axis handle for plotting

Return type:

plt.axis

Raises:

ValueError – If element in info is not recognized

Geometry

class pyraysum.prs.Geometry(baz, slow, dn=0, de=0, maxtr=500)[source]

Ray geometry of seismic events and station configuration.

One set of synthetic traces will be computed for each array element.

Note

This object holds the infromation of the .geom file in classic Raysum

Parameters:

  • baz (float or array_like) – Ray back-azimuths (deg)

  • slow (float or array_like) – Ray slownesses (s/km)

  • dn (float) – North-offset of the seismic station (m)

  • de (float) – East-offset of the seismic station (m)

  • maxtr (int) – Maximum number of traces defined in params.h

Hint

If baz and slow do not have a different length, one ray will be computed for each combination of baz and slow (See examples)

The following attributes are set upon initialization. f prefixes indicate attributes used for interaction with fraysum.run_bare() and fraysum.run_full().

ntr (int):

Number of traces

parameters (list):

Convenience attribute that collects the f-attributes in the order expected by fraysum.run_bare() and fraysum.run_full()

fbaz (np.ndarray):

Ray back-azimuth (radians)

fslow (np.ndarray):

Ray slowness (m/s)

fdn (np.ndarray):

North-offset of the seismic station (m)

fde (np.ndarray):

East-offset of the seismic station (m)

Example

>>> from pyraysum import Geometry
>>> geom = Geometry(60, 0.06)
>>> print(geom)
#Back-azimuth, Slowness, N-offset, E-offset
        60.00    0.0600      0.00      0.00
>>> geom += [[90, 120], 0.04]
>>> print(geom)
#Back-azimuth, Slowness, N-offset, E-offset
        60.00    0.0600      0.00      0.00
        90.00    0.0400      0.00      0.00
       120.00    0.0400      0.00      0.00
>>> geom = Geometry(range(0, 360, 90), [0.04, 0.08], 10, 35)
>>> print(geom)
#Back-azimuth, Slowness, N-offset, E-offset
         0.00    0.0400     10.00     35.00
        90.00    0.0400     10.00     35.00
       180.00    0.0400     10.00     35.00
       270.00    0.0400     10.00     35.00
         0.00    0.0800     10.00     35.00
        90.00    0.0800     10.00     35.00
       180.00    0.0800     10.00     35.00
       270.00    0.0800     10.00     35.00
copy()[source]

Return a copy of the geometry

save(fname='sample.geom')[source]

Alias for write()

write(fname='sample.geom')[source]

Write ray geometry to disk as ascii file

Parameters:

fname (str) – Name of file

plot(show=True)[source]

Plot ray geometry in polar coordinates

Returns:

Axis handle for plotting

Return type:

plt.axis

Control

class pyraysum.prs.Control(verbose=0, wvtype='P', mults=0, npts=300, dt=0.025, align=1, shift=None, rot=1, maxseg=45, maxtr=500, maxph=40000, maxsamp=100000)[source]

Run Control parameters for prs.run().

Parameters:

  • wvtype (str) – Wave type of incoming wavefield ('P', 'SV', or 'SH')

  • mults (int) –

    ID for calculating free surface multiples

    • 0: no multiple

    • 1: First interface multiples only

    • 2: all first-order multiples (once reflected from the surface)

    • 3: supply phases to be computed via Control.set_phaselist()

  • npts (int) – Number of samples in time series

  • dt (float) – Sampling interval in seconds

  • align (int) –

    ID for time alignment of seismograms

    • 0 or “none”: do not align

    • 1 or “P”: align at P

    • 2 or “SV”: align at SV

    • 3 or “SH”: align at SH

  • shift (float or None) – Time shift in seconds (positive shift moves seismograms to greater lags). If None, set to dt, to include direct wave when align > 0.

  • rot (int or str) –

    ID for rotation:

    • 0 or “ZNE”: up, north, east [left-handed]

    • 1 or “RTZ”: radial, transverse, down [positive towards source]

    • 2 or “PVH”: P-, SH-, SV-polarization [positive in ray direction]

  • verbose (int or bool) –

    Verbosity:

    • 0 or False: silent

    • 1 or True: verbose

  • maxseg (int) – Maximum number of segments per ray, as defined in param.h

  • maxtr (int) – Maximum number of traces, as defined in param.h

  • maxph (int) – Maximum number of phases per trace, as defined in param.h

  • maxsamp (int) – Maximum number of samples per trace, as defined in param.h

parameters

Convenience attribute that collects the attributes in the order expected by fraysum.run_bare() and fraysum.run_full()

Type:

list

Warning

The option rot=0 returns the left-handed ZNE coordinates (Z positive up), conventionally used for seismometers. The internal Fortran rountines (everything imported through fraysum) return right-handed NEZ coordinates (Z positive down).

set_phaselist(descriptors, equivalent=False, kinematic=False)[source]

Explicitly set phases to be calculated.

Parameters:

  • descriptors (list of str) –

    List of phase descriptors, where each descriptor is a string of consecutive PHASE SEGMENT pairs, where PHASE is:

    • P upgoing P-wave

    • S upgoing fast S-wave

    • T upgoing slow S-wave

    • p downgoing P-wave

    • s downgoing fast S-wave

    • t downgoing slow S-wave

    SEGMENT is the index of the model layer.

  • equivalent (bool) – Augment phaselist by equivalent phases (e.g., ‘1P0P0s0P’; ‘1P0S0p0P’)

  • kinematic (bool) – Limit equivalent phases to those with same polarity as input phases

Raises:

IndexError – If resulting phaselist is longer than maxph.

Caution

At every interface, S-wave energy split up into a slow and a fast S-wave, S and T, regardless of layer anisotropy. Make sure to include all contributions when setting a explicit phaselist. It is best to inspect the output of Result.descriptors() of a previous run with mults=2 for the ray combinations relevant for your problem.

Caution

Using equivalent=False may produce incorrect amplitudes, because contributions of equivalent phases may be missing. At the same time, amplitudes are oftentimes correct within 10%. Using equivalent=True may cause performance issues.

Hint

In a two layer model (index 0 is the topmost subsurface layer, index 1 is the underlying half-space) the descriptors compute the phases:

  • ['1P0P']: direct P-wave

  • ['1P0S']: P-to-S1 converted wave

  • ['1P0T']: P-to-S2 converted wave

  • ['1P0P', '1P0S', '1P0T]: direct P-wave and all P-to-S converted waves

  • ['1P0P0s0S']: P reflected s at the surface, reflected S at the top of the half-space

See also

Result.descriptors() returns all phases computed in a previous run.

See also

equivalent_phases() outputs equivalent phases explicitly

null_phaselist(mults=2)[source]

Do not use phaselist, but compute using mults keyword.

Parameters:

mults (int) – Set Control.mults to this value.

save(fname='raysum-param', version='prs')[source]

Alias for write()

write(fname='raysum-param', version='prs')[source]

Write parameter file to disk

Parameters:

  • fname – (str) Name of file

  • version – (“prs” or “raysum”) Pyraysum or Raysum raysum compatible

Result

class pyraysum.prs.Result(model=None, geometry=None, rc=[], streams=[])[source]

Result of a PyRaysum wavefield simulation. 3-component synthetic seismograms are strored as a list of streams in the stream attribute. Includes methods to calculate receiver functions, filter and plot the results. See run() for examples.

Parameters:

  • model (Model) – Subsurface velocity model

  • geometry (Geometry) – Recording geometry

  • rc (Control) – Run-control parameters

  • streams (list) – List of Stream objects.

If created with mode='full' in run(), the stats attribute of each Trace in each Stream in Result.streams holds the additional attributes:

phase_times

Arrival times of seismic phases

phase_amplitudes

Amplitudes of seismic phases

phase_descriptors

Descriptors of seismic phases. Index indicates layer through which phase propagates. (e.g. “1P0S”, see also descriptors())

phase_names

Short names of seismic phases (e.g. “PS”)

conversion_names

Conversion names of seismic phases. Index indicates top of layer at which conversion occurs. (e.g. “P1S”)

calculate_rfs()[source]

Generate receiver functions from spectral division of displacement traces. Will be stored in rflist.

Raises:

ValueError – In case Control parameters a unsuitable.

descriptors()[source]
Returns:

list

Unique list of all phase descriptors present

Raises:

AttributeError – If no descriptors are present

Example

>>> from pyraysum import Model, Control, Geometry, run
>>> model = Model([30000., 0], [2800., 3300.], [6000., 8000.], [3600., 4500.])
>>> geom = Geometry(0., 0.06)
>>> rc = Control(mults=0)
>>> seismogram = run(model, geom, rc)
>>> seismogram.descriptors()
['1P0P', '1P0S']
write(fname='sample.tr')[source]

Write streams to file, using legacy Raysum trace format

Parameters:

fname (str) – Name of the output file (including extension)

plot(typ, **kwargs)[source]

Plot the displacement seismograms and/or receiver functions stored in Result streams.

Parameters:

  • typ (str) – Type of plot to show. Options are 'streams', 'rfs', or 'all' for the displacement seismograms, receiver functions, or both

  • **kwargs – Are passed to the underlying stream_wiggles() or rf_wiggles()

Hint

If Result contains only one event (i.e., single baz and slow values), it is more appropriate to plot the waveforms using the default Stream.plot() method on Result.streams[0].

filter(typ, ftype, **kwargs)[source]

Filters the displacement seismograms and/or receiver functions stored in Result streams.

Parameters:

  • typ (str) – Type of plot to show. Options are 'streams', 'rfs', or 'all' for the displacement seismograms, receiver functions, or both

  • ftype (str) – Type of filter to use in obspy.Trace.filter

  • **kwargs –

    Keyword arguments passed to obspy.Trace.filter

Modules

Module prs

pyraysum.prs.run(model, geometry, rc, mode='full', rf=False)[source]

Run a wave-field simulation. This function calls the compiled fraysum binaries.

Parameters:

  • model (Model) – Subsurface velocity model

  • geometry (Geometry) – Recording geometry

  • rc (Control) – Run-control parameters

  • mode (str) –

    • 'full': Compute seismograms, phase arrivals and descriptors (slower)

    • 'bare': Only compute seismograms (faster)

  • rf (bool) – Whether or not to calculate receiver functions

Returns:

Synthetic seimograms

Return type:

Result

Hint

Best performance is achived by calling fraysum.run_bare()

Example

>>> from pyraysum import prs, Model, Geometry, Control
>>> # Define two-layer model with isotropic crust over isotropic half-space
>>> model = Model([30000., 0], [2800., 3300.], [6000., 8000.], [3600., 4500.])
>>> geom = Geometry(0., 0.06) # baz = 0 deg; slow = 0.06 s/km
>>> rc = Control(npts=1500, dt=0.025)
>>> result = prs.run(model, geom, rc, rf=True)
>>> type(result.streams[0])
<class 'obspy.core.stream.Stream'>
>>> type(result.rfs[0])
<class 'obspy.core.stream.Stream'>
>>> seis, rf = result[0]
>>> print(seis)
3 Trace(s) in Stream:
...SYR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:37.475000Z | 40.0 Hz, 1500 samples
...SYT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:37.475000Z | 40.0 Hz, 1500 samples
...SYZ | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:37.475000Z | 40.0 Hz, 1500 samples
>>> print(rf)
2 Trace(s) in Stream:
...RFR | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:37.475000Z | 40.0 Hz, 1500 samples
...RFT | 1970-01-01T00:00:00.000000Z - 1970-01-01T00:00:37.475000Z | 40.0 Hz, 1500 samples
pyraysum.prs.read_model(modfile, encoding=None, version='prs')[source]

Reads model parameters from file

Parameters:

version ("prs" or "raysum") – File has PyRaysum or Raysum format

Returns:

Seismic velocity model

Return type:

Model

pyraysum.prs.read_geometry(geomfile, encoding=None)[source]

Reads geometry parameters from file

Returns:

Ray geometry

Return type:

Geometry

pyraysum.prs.read_control(paramfile, version='prs')[source]

Read Raysum run control parameters from file.

Parameters:

version – (“prs” or “raysum”) File has PyRaysum or Raysum format

Returns:

Run control parameters

Return type:

Control

pyraysum.prs.equivalent_phases(descriptors, kinematic=False)[source]

Return descriptors of equivalent phases, i.e. those arriving at the same time as the input phases

Parameters:

  • descriptors (list of strings) – Phase descriptors as in Control.set_phaselist()

  • kinematic (boolean) – If True, restrict to kinematically equivalent phases, i.e. those that have the same polarization as input phases

Returns:

List of unique descriptors of equivalent phases

Return type:

list of strings

Example

>>> from pyraysum import prs
>>> prs.equivalent_phases(['1P0P0p0S'])
['1P0P0s0P', '1P0S0p0P']

Module frs

pyraysum.frs.read_traces(traces, rc, geometry, arrivals=None)[source]

Create a Result from the array produced by fraysum.run_bare() and fraysum.run_full().

Parameters:

  • rc (Control) – Run-control parameters

  • geometry (Geometry) – Ray parameters

  • arrivals (list) – Output of read_arrivals(). List of arrival times, amplitudes, and names

Returns:

List of obspy.Stream objects

Return type:

list

pyraysum.frs.read_arrivals(ttimes, amplitudes, phaselist, geometry)[source]

Convert the phaselist, amplitude and traveltime output of fraysum.run_full() to lists of phase arrival times, amplitudes, long phase descriptors, short phase names, and conversion names.

Parameters:

  • ttimes (np.array) – Travel time array …

  • amplitudes (np.array) – Amplitude array …

  • phaselist (np.array) – Phase identifier array returned by fraysum.run_full()

  • geometry (prs.Geometry) – Ray geometry

Returns:

3-component phase arrival lists, where indices translate to:

  • 0: phase arrival times

  • 1: phase amplitudes

  • 2: (long) phase descriptors, e.g. “1P0S”

  • 3: (short) phase names, e.g. “PS”

  • 4: (intermediate) conversion names, e.g. “P1S”

Return type:

list

pyraysum.frs.make_array(geometry, rc)[source]

Initialize array for NumPy-based post processing

Parameters:

  • geometry (prs.Geometry) – Recording geometry

  • rc (prs.Control) – Run-control parameters

Returns:

Array in shape to be used by filtered_rf_array() and filtered_array().

Return type:

numpy.zeros((geometry.ntr, 2, rc.npts))

pyraysum.frs.filtered_rf_array(traces, rfarray, ntr, npts, dt, fmin, fmax)[source]

Fast, NumPy-based, receiver function computation and filtering of fraysum.run_bare() output

Roughly equivalent to subsequent calls to read_traces(), Result.calculate_rfs() and Result.filter(), stripped down for computational efficiency, for use in inversion/probabilistic approaches.

Parameters:

  • traces (np.ndarray) – Output of fraysum.run_bare()

  • rfarray (np.ndarray) – Initialized array of shape (ntr, 2, npts) to store output (See: make_array())

  • ntr (int) – Number of traces (Geometry.ntr)

  • npts (int) – Number of points per trace (Control.npts)

  • dt (float) – Sampling interval (Control.dt)

  • fmin (float) – Lower bandpass frequency corner (Hz)

  • fmax (float) – Upper bandpass frequency corner (Hz)

Returns:

None. Output is written to rfarray

Return type:

None

Warning

Assumes PVH alignment (ray-polarization), i.e. Control.rot="PVH".

pyraysum.frs.filtered_array(traces, rfarray, ntr, npts, dt, fmin, fmax)[source]

Fast, NumPy-based filtering of fraysum.run_bare() output

Roughly equivalent to subsequent calls to read_traces(), and Result.filter(), stripped down for computational efficiency, for use in inversion/probabilistic approaches.

Parameters:

  • traces (np.ndarray) – Output of fraysum.run_bare()

  • rfarray (np.ndarray) – Initialized array of shape (ntr, 2, npts) to store output (See: make_array())

  • ntr (int) – Number of traces (Geometry.ntr)

  • npts (int) – Number of points per trace (Control.npts)

  • dt (float) – Sampling intervall (Control.dt)

  • fmin (float) – Lower bandpass frequency corner (Hz)

  • fmax (float) – Upper bandpass frequency corner (Hz)

Returns:

None. Output is written to rfarray

Return type:

None

Warning

Assumes PVH alignment (ray-polarization), i.e. Control.rot="PVH".

Module plot

Functions to plot single station displacement or receiver function seismograms.

pyraysum.plot.stack_all(stream, pws=False)[source]

Stack all traces in Stream objects.

Parameters:

  • stream (Stream) – Contains traces to stack

  • pws (bool) – Enables Phase-Weighted Stacking

Returns:

Stacked trace

Return type:

Trace

pyraysum.plot.rf_wiggles(rflist, btyp='baz', wvtype='P', pws=False, tmin=-5.0, tmax=20, scale=None, pcolor='red', ncolor='blue', save=False, axs=None, ftitle='Figure_rf_wiggle', fmt='png', plot_kwargs={'color': 'black', 'linewidth': 0.1}, figure_kwargs={})[source]

Plot receiver function seismograms sorted by back-azimuth or slowness.

Parameters:

  • rflist (list or prs.Result) – list of obspy.core.Stream objects containing receiver functions

  • btyp (str) – Type of sorting for panel

  • wvtype (str) – Wavet type (‘P’, ‘SV’, or ‘SH’)

  • pws (bool) – Enables Phase-Weighted Stacking

  • tmin (float) – Lower bound of time axis (s)

  • tmax (float) – Upper bound of time axis (s)

  • scale (float) – Amplitude scaling factor

  • save (bool) – Whether or not to save the figure

  • axs (4*tuple) –

    matplotlib axes to place the subplots.

    • axs[0]: radial stack

    • axs[1]: radial section

    • axs[2]: transverse stack

    • axs[3]: transverse section

  • pcolor (str) – Color to fill positive wiggles

  • ncolor (str) – Color to fill negative wiggles

  • ftitle (str) – Filename of figure to be saved (without format suffix, see fmt)

  • fmt (str) – Output format (‘png’, ‘jpg’, or ‘eps’, etc.)

  • plot_kwargs (dict) – Keyword arguments passed to matplotlib.pyplot.plot()

  • figure (dict) – Keyword arguments passed to matplotlib.pyplot.figure()

Returns:

Figure handle

Return type:

matplotlib.pyplot.figure

pyraysum.plot.stream_wiggles(streamlist, btyp='baz', wvtype='P', tmin=-5.0, tmax=20.0, scale=None, pcolor='red', ncolor='blue', save=False, axs=None, ftitle='Figure_pw_wiggles', fmt='png', plot_kwargs={'color': 'black', 'linewidth': 0.1}, figure_kwargs={})[source]

Plot displacement seismograms sorted by back-azimuth or slowness.

Parameters:

  • streamlist (list or prs.Result) –

    list of obspy.core.Stream objects containing displacement seismograms

  • btyp (str) – Type of sorting for panel

  • wvtype (str) – Wave type (‘P’, ‘SV’, or ‘SH’)

  • tmin (float) – Lower bound of time axis (s)

  • tmax (float) – Upper bound of time axis (s)

  • scale (float) – Scaling factor

  • pcolor (str) – Color to fill positive wiggles

  • ncolor (str) – Color to fill negative wiggles

  • save (bool) – Whether or not to save the figure

  • axs (3*tuple) –

    matplotlib axes to place the subplots.

    • axs[0]: P or R section

    • axs[1]: V or T section

    • axs[2]: H or Z section

  • ftitle (str) – Filename of figure to be saved (without format suffix, see fmt)

  • fmt (str) – Output format (‘png’, ‘jpg’, or ‘eps’, etc.)

  • plot_kwargs (dict) –

    Keyword arguments passed to matplotlib.pyplot.plot()

  • figure (dict) –

    Keyword arguments passed to matplotlib.pyplot.figure()

Returns:

Figure handle

Return type:

matplotlib.pyplot.figure

pyraysum.plot.seis_wiggles(stream, tmin=-5.0, tmax=20.0, save=False, ftitle='Figure_pw_wiggles_3c', fmt='png', figure_kwargs={})[source]

Plots 3-component wiggles.

Parameters:

  • stream (Stream) –

    obspy.core.Stream containing 3 traces

  • tmin (float) – Lower bound of time axis (s)

  • tmax (float, optional) – Upper bound of time axis (s)

  • save (bool) – Whether or not to save the figure

  • ftitle (str) – Filename of figure to be saved (without format suffix, see fmt)

  • fmt (str) – Output format (‘png’, ‘jpg’, or ‘eps’, etc.)

Returns:

Figure handle

Return type:

matplotlib.pyplot.figure

Module fraysum

These are the the access points to the low-level fortran routines.

Tip

The corresponding functions can be more succinctly called with the unpacked lists of parameters, like this: fraysum.run_bare(*model.parameters, *geometry.parameters, *rc.parameters).

tr_ph = fraysum.run_bare(thick, rho, alpha, beta, isoflag, pct, trend,plunge, strike, dip, nlay, baz, slow, sta_dx, sta_dy, ntr, iphase, mults, nsamp,dt, align, shift, out_rot, verb, nsegin, numphin, phaselistin)
Parameters:

  • thick – Model.fthickn

  • rho – Model.frho

  • alpha – Model.fvp

  • beta – Model.fvs

  • isoflag – Model.fflag

  • pct – Model.fani

  • trend – Model.ftrend

  • plunge – Model.fplunge

  • strike – Model.fstrike

  • dip – Model.fdip

  • nlay – Model.nlay

  • baz – Geometry.fbaz

  • slow – Geometry.fslow

  • sta_dx – Geometry.fdn

  • sta_dy – Geometry.fde

  • ntr – Geometry.ntr

  • iphase – Control.wvtype

  • mults – Control.mults

  • nsamp – Control.nsamp

  • dt – Control.dt

  • align – Control.align

  • shift – Control.shift

  • out_rot – Control.rot

  • verb – Control.verbose

  • nsegin – Control._nseg

  • numphin – Control._numph

  • phaselistin – Control.phaselist

Returns:

tr_ph

Return type:

rank-3 array(‘f’) with bounds (3, maxsamp, maxtr)

tr_ph, travel_time, amplitudeout, phaselist = fraysum.run_full(thick,rho, alpha, beta, isoflag, pct, trend, plunge, strike, dip, nlay, baz, slow,sta_dx, sta_dy, ntr, iphase, mults, nsamp, dt, align, shift, out_rot, verb, nsegin,numphin, phaselistin)
Parameters:

  • thick – Model.fthickn

  • rho – Model.frho

  • alpha – Model.fvp

  • beta – Model.fvs

  • isoflag – Model.fflag

  • pct – Model.fani

  • trend – Model.ftrend

  • plunge – Model.fplunge

  • strike – Model.fstrike

  • dip – Model.fdip

  • nlay – Model.nlay

  • baz – Geometry.fbaz

  • slow – Geometry.fslow

  • sta_dx – Geometry.fdn

  • sta_dy – Geometry.fde

  • ntr – Geometry.ntr

  • iphase – Control.wvtype

  • mults – Control.mults

  • nsamp – Control.nsamp

  • dt – Control.dt

  • align – Control.align

  • shift – Control.shift

  • out_rot – Control.rot

  • verb – Control.verbose

  • nsegin – Control._nseg

  • numphin – Control._numph

  • phaselistin – Control.phaselist

Returns:

tr_ph

Return type:

rank-3 array(‘f’) with bounds (3, maxsamp, maxtr):

Returns:

travel_time

Return type:

rank-2 array(‘f’) with bounds (maxph, maxtr)

Returns:

amplitudeout

Return type:

rank-3 array(‘f’) with bounds (3, maxph,`maxtr`)

Returns:

phaselist

Return type:

rank-3 array(‘i’) with bounds (maxseg, 2, maxph)