seispy package

seispy.ccpprofile module

class seispy.ccpprofile.CCPProfile(cfg_file=None, log: ~seispy.setuplog.SetupLog = <seispy.setuplog.SetupLog object>)

Bases: object

Methods

initial_profile()	Initialize bins of profile
read_rfdep()	Read RFdepth file
save_stack_data([format])	If format is npz, saving stacked data and parameters to local as a npz file.
stack()	Stack RFs in bins

initial_profile()

Initialize bins of profile

read_rfdep()

Read RFdepth file

Raises:

FileNotFoundError – Not Found RFdepth file

save_stack_data(format='npz')

If format is npz, saving stacked data and parameters to local as a npz file. To load the file, please use data = np.load(fname, allow_pickle=True). data['cpara'] is the parameters when CCP stacking. data['stack_data'] is the result of stacked data.

If format is dat the stacked data will be save into a txt file with 8 columns, including bin_lat, bin_lon, profile_dis, depth, amp, ci_low, ci_high and count.

• bin_lat and bin_lon represent the position of each bin;

• profile_dis represents the distance in km between each bin and the start point of the profile;

• depth represents depth of each bin;

• amp means the stacked amplitude;

• ci_low and ci_high mean confidence interval with bootstrap method;

• count represents stacking number of each bin.

Parameters:

format (str) – Format for stacked data

stack()

Stack RFs in bins

seispy.ccpprofile.create_center_bin_profile(stations, val=5, method='linear')

Create bins along a profile with given stations

Parameters:

• stations (seispy.rf2depth_makedata.Station) – Stations along a profile

• val (float) – The interval between two points in km

• method (str) – Method for interpolation

Returns:

The location of bins (bin_loca), and length between each bin and the start point (profile_range)

Return type:

(numpy.array, numpy.array)

seispy.ccpprofile.init_profile(lat1, lon1, lat2, lon2, val)

Initial bins along a profile with given position of two points.

Parameters:

• lat1 (float) – The latitude of the start point

• lon1 (float) – The lontitude of the start point

• lat2 (float) – The latitude of the end point

• lon2 (float) – The lontitude of the end point

• val (float) – The interval between two points in km

Returns:

The location of bins (bin_loca), and length between each bin and the start point (profile_range)

The bin_loca is positions of bins with a numpy.array with two column. The profile_range is distance between bin center and the start point with an 1D numpy.array. :rtype: (numpy.array, numpy.array)

seispy.cc3d module

class seispy.ccp3d.CCP3D(cfg_file=None, log=None)

Bases: object

Class for 3-D CCP stacking, Usually used to study mantle transition zone structure.

Parameters:

• cfg_file (str, optional) – Path to configure file. If not defined a instance of CCP3D.cpara will be initialed, defaults to None

• log (seispy.sutuplog.logger , optional) – A logger instance. If not defined, seispy.sutuplog.logger will be initialed, defaults to None

Methods

initial_grid()	Initial grid points and search stations within a distance.
load_para(cfg_file)	Load parameters from configure file.
read_rfdep()	Read RFdepth data from npz file.
read_stack_data(stack_data_path[, cfg_file, ...])	Read stacked data from local.
save_good_410_660(fname)	Save good 410 and 660 km discontinuities to local.
save_stack_data(fname)	Save stacked data and parameters to local as a npz file.
search_good_410_660([peak_410_min, ...])	Search good 410 and 660 km discontinuities from stacked data.
stack()	Search conversion points falling within a bin and stack them with bootstrap method.

get_depth_err

get_depth_err(type='std')

initial_grid()

Initial grid points and search stations within a distance.

load_para(cfg_file)

Load parameters from configure file.

Parameters:

cfg_file (str) – Path to configure file.

read_rfdep()

Read RFdepth data from npz file.

classmethod read_stack_data(stack_data_path, cfg_file=None, good_depth_path=None, ismtz=False)

Read stacked data from local.

Parameters:

• stack_data_path (str) – Path to stacked data.

• cfg_file (str, optional) – Path to configure file, defaults to None

• good_depth_path (str, optional) – Path to good depth file, defaults to None

• ismtz (bool, optional) – Whether the good depth file is in mtz format, defaults to False

Returns:

A instance of CCP3D.

Return type:

CCP3D

save_good_410_660(fname)

Save good 410 and 660 km discontinuities to local.

Parameters:

fname (str) – file name of good 410 and 660 km discontinuities

save_stack_data(fname)

Save stacked data and parameters to local as a npz file. To load the file, please use data = np.load(fname, allow_pickle=True). data[‘cpara’] is the parameters when CCP stacking. data[‘stack_data’] is the result of stacked data.

Parameters:

fname (str) – file name of stacked data

search_good_410_660(peak_410_min=380, peak_410_max=440, peak_660_min=630, peak_660_max=690)

Search good 410 and 660 km discontinuities from stacked data.

Parameters:

• peak_410_min (float, optional) – Minimum depth of 410 km discontinuity, defaults to 380

• peak_410_max (float, optional) – Maximum depth of 410 km discontinuity, defaults to 440

• peak_660_min (float, optional) – Minimum depth of 660 km discontinuity, defaults to 630

• peak_660_max (float, optional) – Maximum depth of 660 km discontinuity, defaults to 690

stack()

Search conversion points falling within a bin and stack them with bootstrap method.

seispy.ccp3d.bin_shape(cpara)

Compute the radius of bins in degree.

Parameters:

cpara (CCPPara) – Parameters of CCP stacking.

Returns:

Radius of bins in degree.

Return type:

1-D ndarray of floats with shape (n, ), where n is the number of bins.

seispy.ccp3d.boot_bin_stack(data_bin, n_samples=3000)

Stack data with bootstrap method.

Parameters:

• data_bin (1-D ndarray of floats) – Data falling within a bin.

• n_samples (int, optional) – Number of bootstrap samples, defaults to 3000

Returns:

Mean, confidence interval and number of data falling within a bin.

Return type:

tuple of floats and int

seispy.ccp3d.gen_center_bin(center_lat, center_lon, len_lat, len_lon, val)

Create spaced grid point with coordinates of the center point in the area in spherical coordinates.

Parameters:

• center_lat (float) – Latitude of the center point.

• center_lon (float) – Longitude of the center point.

• len_lat (float) – Half length in degree along latitude axis.

• len_lon (float) – Half length in degree along longitude axis.

• val (float) – Interval in degree between adjacent grid point.

Returns:

Coordinates of Grid points.

Return type:

2-D ndarray of floats with shape (n, 2), where n is the number of grid points.

seispy.ccppara module

class seispy.ccppara.CCPPara

Bases: object

Attributes:

bin_radius

shape

property bin_radius

property shape

seispy.ccppara.ccppara(cfg_file)

Read configure file for CCP stacking :param cfg_file: Path to configure file :type cfg_file: str :return: CCPPara object :rtype: CCPPara

seispy.decon module

class seispy.decon.RFTrace(data=Ellipsis, header=None)

Bases: Trace

Class for receiver function trace.

Attributes:

id

Return a SEED compatible identifier of the trace.

meta

Methods

attach_response(inventories)	Search for and attach channel response to the trace as obspy.core.trace.Trace.stats.response.
copy()	Returns a deepcopy of the trace.
count()	Return number of data samples of the current trace.
decimate(factor[, no_filter, strict_length])	Downsample trace data by an integer factor.
deconvolute(utr, wtr[, method])	Deconvolute receiver function from waveforms.
detrend([type])	Remove a trend from the trace.
differentiate([method])	Differentiate the trace with respect to time.
filter(type, **options)	Filter the data of the current trace.
get_id()	Return a SEED compatible identifier of the trace.
integrate([method])	Integrate the trace with respect to time.
interpolate(sampling_rate[, method, ...])	Interpolate the data using various interpolation techniques.
max()	Returns the value of the absolute maximum amplitude in the trace.
normalize([norm])	Normalize the trace to its absolute maximum.
plot(**kwargs)	Create a simple graph of the current trace.
remove_response([inventory, output, ...])	Deconvolve instrument response.
remove_sensitivity([inventory])	Remove instrument sensitivity.
resample(sampling_rate[, window, no_filter, ...])	Resample trace data using Fourier method.
simulate([paz_remove, paz_simulate, ...])	Correct for instrument response / Simulate new instrument response.
slice([starttime, endtime, nearest_sample])	Return a new Trace object with data going from start to end time.
slide(window_length, step[, offset, ...])	Generator yielding equal length sliding windows of the Trace.
spectrogram(**kwargs)	Create a spectrogram plot of the trace.
split()	Split Trace object containing gaps using a NumPy masked array into several traces.
std()	Method to get the standard deviation of amplitudes in the trace.
taper(max_percentage[, type, max_length, side])	Taper the trace.
times([type, reftime])	For convenient plotting compute a NumPy array with timing information of all samples in the Trace.
trigger(type, **options)	Run a triggering algorithm on the data of the current trace.
trim([starttime, endtime, pad, ...])	Cut current trace to given start and end time.
verify()	Verify current trace object against available meta data.
write(filename[, format])	Save current trace into a file.

classmethod deconvolute(utr, wtr, method='iter', **kwargs)

Deconvolute receiver function from waveforms.

Parameters:

• utr (obspy.Trace) – R or Q component for the response function

• wtr (obspy.Trace) – Z or L component for the source function

• method (str, optional) – Method for deconvolution, defaults to ‘iter’

• kwargs (dict) – Parameters for deconvolution

Returns:

RFTrace object

Return type:

RFTrace

seispy.decon.correl(R, W, nfft)

Correlation in frequency domain.

Parameters:

• R (np.ndarray) – numerator

• W (np.ndarray) – denominator

Returns:

Correlation in frequency domain

Return type:

np.ndarray

seispy.decon.deconit(uin, win, dt, nt=None, tshift=10, f0=2.0, itmax=400, minderr=0.001, phase='P')

Iterative deconvolution using Ligorria & Ammon method. @author: Mijian Xu @ NJU Created on Wed Sep 10 14:21:38 2014

Parameters:

• uin (np.ndarray) – R or Q component for the response function

• win (np.ndarray) – Z or L component for the source function

• dt (float) – sample interval in second

• nt (int, optional) – number of samples, defaults to None

• tshift (float, optional) – Time shift before P arrival, defaults to 10.

• f0 (float, optional) – Gauss factor, defaults to 2.0

• itmax (int, optional) – Max iterations, defaults to 400

• minderr (float, optional) – Min change in error required for stopping iterations, defaults to 0.001

• phase (str, optional) – Phase of the RF, defaults to ‘P’

Returns:

(RFI, rms, it) RF, rms and number of iterations.

Return type:

(np.ndarray, np.ndarray, int)

seispy.decon.deconvolute(uin, win, dt, method='iter', **kwargs)

Deconvolute receiver function from waveforms. :param uin: R or Q component for the response function :type uin: np.ndarray :param win: Z or L component for the source function :type win: np.ndarray :param dt: sample interval in second :type dt: float :param method: Method for deconvolution, defaults to ‘iter’:type method: str, optional :param kwargs: Parameters for deconvolution :type kwargs: dict :return: RF, rms, [iter] :rtype: np.ndarray, float, int

seispy.decon.deconwater(uin, win, dt, tshift=10.0, wlevel=0.05, f0=2.0, normalize=False, phase='P')

Frequency-domain deconvolution using waterlevel method.

Parameters:

• uin (np.ndarray) – R or Q component for the response function

• win (np.ndarray) – Z or L component for the source function

• dt (float) – sample interval in second

• tshift (float, optional) – Time shift before P arrival, defaults to 10.

• wlevel (float, optional) – Waterlevel to stabilize the deconvolution, defaults to 0.05

• f0 (float, optional) – Gauss factor, defaults to 2.0

• normalize (bool, optional) – If normalize the amplitude of the RF, defaults to False

Returns:

(rf, rms) RF and final rms.

Return type:

(np.ndarray, float)

seispy.decon.gaussFilter(dt, nft, f0)

Gaussian filter in frequency domain.

Parameters:

• dt (float) – sample interval in second

• nft (int) – number of samples

• f0 (float) – Gauss factor

Returns:

Gaussian filter in frequency domain

Return type:

np.ndarray

seispy.decon.gfilter(x, nfft, gauss, dt)

Apply Gaussian filter on time series.

Parameters:

• x (np.ndarray) – input trace

• nfft (int) – number of samples

• gauss (np.ndarray) – Gaussian filter in frequency domain

• dt (float) – sample interval in second

Returns:

Filtered data in time domain

Return type:

np.ndarray

seispy.decon.phaseshift(x, nfft, dt, tshift)

Phase shift in frequency domain.

Parameters:

• x (np.ndarray) – input trace

• nfft (int) – number of samples

• dt (float) – sample interval in second

• tshift (float) – Time shift before P arrival

Returns:

Phase shifted data in time domain

Return type:

np.ndarray

seispy.distaz module

class seispy.distaz.distaz(lat1, lon1, lat2, lon2)

Bases: object

Subroutine to calculate the Great Circle Arc distance

between two sets of geographic coordinates

Equations take from Bullen, pages 154, 155

20. Owens, September 19, 1991

    Sept. 25 – fixed az and baz calculations

P. Crotwell, Setember 27, 1995 Converted to c to fix annoying problem of fortran giving wrong answers if the input doesn’t contain a decimal point.

H. P. Crotwell, September 18, 1997 Java version for direct use in java programs. * * C. Groves, May 4, 2004 * Added enough convenience constructors to choke a horse and made public double * values use accessors so we can use this class as an immutable

H.P. Crotwell, May 31, 2006 Port to python, thus adding to the great list of languages to which distaz has been ported from the origin fortran: C, Tcl, Java and now python and I vaguely remember a perl port. Long live distaz!

Mijian Xu, Jan 01, 2016 Add np.ndarray to available input.

Methods

degreesToKilometers
getAz
getBaz
getDelta

az

c c Make sure 0.0 is always 0.0, not 360. c

baz

c c Bullen, Sec 10.2, eqn 7 / eqn 8 c c pt. 2 is unprimed, so this is technically the az c

degreesToKilometers()

delta

c c Bullen, Sec 10.2, eqn 7 / eqn 8 c c pt. 1 is unprimed, so this is technically the baz c c Calculate baz this way to avoid quadrant problems c

evtlon

if (lat1 == lat2) and (lon1 == lon2):

self.delta = 0.0 self.az = 0.0 self.baz = 0.0 return

getAz()

getBaz()

getDelta()

seispy.eq module

class seispy.eq.EQ(pathname, datestr, suffix='SAC')

Bases: object

Methods

arr_correct([write_to_sac])	offset = sac.b - real o
channel_correct([switchEN, reverseE, reverseN])	Correct channel name for R, E, N components
decon_p(tshift[, tcomp])	Deconvolution for P wave
decon_s(tshift, **kwargs)	Deconvolution for S wave
deconvolute(shift, time_after[, f0, method, ...])	Deconvolution
detrend()	Detrend and demean
filter([freqmin, freqmax, order])	Bandpass filter
fix_channel_name()	Fix channel name for R, E, N components
from_stream(stream)	Create EQ object from obspy stream
get_arrival(model, evdp, dis[, phase])	Get arrival time, ray parameter and incident angle from TauP model
get_time_offset([event_time])	Get time offset from SAC header
judge_rf(gauss, shift, npts[, criterion, ...])	Judge whether receiver function is valid
phase_trigger(time_before, time_after[, ...])	Trigger P or S phase
readstream()	Read SAC files to stream
rotate(baz[, inc, method, search_inc, baz_shift])	Rotate to radial and transverse components
saverf(path[, evtstr, shift, evla, evlo, ...])	Save receiver function to SAC file
search_baz(bazi[, time_b, time_e, offset])	Search back azimuth for P wave
search_inc(bazi)	Search incident angle for S wave
set_comp()	Set component name
snr([length])	Calculate SNR
trim(time_before, time_after[, isreturn])	offset = sac.b - real o
write(path, evt_datetime)	Write raw stream to SAC files

arr_correct(write_to_sac=True)

offset = sac.b - real o

channel_correct(switchEN=False, reverseE=False, reverseN=False)

Correct channel name for R, E, N components

Parameters:

• switchEN (bool, optional) – _description_, defaults to False

• reverseE (bool, optional) – _description_, defaults to False

• reverseN (bool, optional) – _description_, defaults to False

decon_p(tshift, tcomp=False, **kwargs)

Deconvolution for P wave

Parameters:

• tshift (float) – Time shift before P arrival

• tcomp (bool, optional) – Whether calculate transverse component, defaults to False

decon_s(tshift, **kwargs)

Deconvolution for S wave

Parameters:

tshift (float) – Time shift before P arrival

deconvolute(shift, time_after, f0=2.0, method='iter', only_r=False, itmax=400, minderr=0.001, wlevel=0.05, target_dt=None)

Deconvolution

Parameters:

shiftfloat

Time shift before P arrival

time_afterfloat

Time length after P arrival

f0float or list, optional

Gaussian factors, by default 2.0

methodstr, optional

method for deconvolution in iter or water, by default iter

only_rbool, optional

Whether only calculate RF in prime component, by default False

itmaxint, optional

Maximum iterative number, valid for method of iter, by default 400

minderrfloat, optional

Minium residual error, valid for method of iter, by default 0.001

wlevelfloat, optional

Water level, valid for method of water, by default 0.05

target_dtNone or float, optional

Time delta for resampling, by default None

detrend()

Detrend and demean

filter(freqmin=0.05, freqmax=1, order=4)

Bandpass filter

Parameters:

• freqmin (float, optional) – minimum frequency, defaults to 0.05

• freqmax (float, optional) – maximum frequency, defaults to 1

• order (int, optional) – filter order, defaults to 4

fix_channel_name()

Fix channel name for R, E, N components

classmethod from_stream(stream)

Create EQ object from obspy stream

Parameters:

stream (obspy.Stream) – obspy stream

Returns:

EQ object

Return type:

EQ

get_arrival(model, evdp, dis, phase='P')

Get arrival time, ray parameter and incident angle from TauP model

Parameters:

• model (TauPyModel) – TauP model

• evdp (float) – focal depth

• dis (float) – epicentral distance

• phase (str, optional) – phase name, defaults to ‘P’

get_time_offset(event_time=None)

Get time offset from SAC header

Parameters:

event_time (obspy.core.utcdatetime.UTCDateTime, optional) – event time, defaults to None

judge_rf(gauss, shift, npts, criterion='crust', rmsgate=None)

Judge whether receiver function is valid

Parameters:

• gauss (float) – Gaussian factor

• shift (float) – time shift before P arrival

• npts (int) – number of points for RF

• criterion (str, optional) – criterion for judging, defaults to ‘crust’

• rmsgate (float, optional) – RMS gate, defaults to None

Returns:

whether RF is valid

Return type:

bool

phase_trigger(time_before, time_after, prepick=True, stl=5, ltl=10)

Trigger P or S phase

Parameters:

• time_before (float) – time before P or S arrival

• time_after (float) – time after P or S arrival

• prepick (bool, optional) – whether use prepick, defaults to True

• stl (int, optional) – short time length for STA/LTA, defaults to 5

• ltl (int, optional) – long time length for STA/LTA, defaults to 10

readstream()

Read SAC files to stream

rotate(baz, inc=None, method='NE->RT', search_inc=False, baz_shift=0)

Rotate to radial and transverse components

Parameters:

• baz (float) – back azimuth

• inc (float, optional) – incident angle, defaults to None

• method (str, optional) – method for rotation, defaults to ‘NE->RT’

• search_inc (bool, optional) – whether search incident angle, defaults to False

• baz_shift (int, optional) – shift back azimuth, defaults to 0

saverf(path, evtstr=None, shift=0, evla=-12345.0, evlo=-12345.0, evdp=-12345.0, mag=-12345.0, gauss=0, baz=-12345.0, gcarc=-12345.0, only_r=False, **kwargs)

Save receiver function to SAC file

Parameters:

• path (str) – path to save SAC file

• evtstr (str, optional) – event string, defaults to None

• shift (int, optional) – time shift before P arrival, defaults to 0

• evla (float, optional) – event latitude, defaults to -12345.

• evlo (float, optional) – event longitude, defaults to -12345.

• evdp (float, optional) – event depth, defaults to -12345.

• mag (float, optional) – event magnitude, defaults to -12345.

• gauss (float, optional) – Gaussian factor, defaults to 0

• baz (float, optional) – back azimuth, defaults to -12345.

• gcarc (float, optional) – epicentral distance, defaults to -12345.

• only_r (bool, optional) – whether only save R component, defaults to False

search_baz(bazi, time_b=10, time_e=20, offset=90)

Search back azimuth for P wave

Parameters:

• bazi (float) – back azimuth

• time_b (int, optional) – time before P arrival, defaults to 10

• time_e (int, optional) – time after P arrival, defaults to 20

• offset (int, optional) – offset for searching, defaults to 90

Returns:

back azimuth and amplitude

Return type:

(float, np.ndarray)

search_inc(bazi)

Search incident angle for S wave

Parameters:

bazi (float) – back azimuth

Returns:

incident angle

Return type:

float

set_comp()

Set component name

snr(length=50)

Calculate SNR

Parameters:

length (int, optional) – length for noise, defaults to 50

Returns:

SNR of E, N, Z components

Return type:

(float, float, float)

trim(time_before, time_after, isreturn=False)

offset = sac.b - real o

write(path, evt_datetime)

Write raw stream to SAC files

Parameters:

• path (str) – path to save SAC files

• evt_datetime (obspy.core.utcdatetime.UTCDateTime) – event datetime

exception seispy.eq.NotEnoughComponent(matchkey)

Bases: Exception

exception seispy.eq.TooMoreComponents(matchkey)

Bases: Exception

seispy.eq.rotateZNE(st)

Rotate Z, N, E components to Z, N, E components

seispy.geo module

seispy.geo.acosd(x)

Inverse cosine function with degree as output

Parameters:

x (float) – Cosine value

Returns:

Degree

Return type:

float

seispy.geo.asind(x)

Inverse sine function with degree as output

Parameters:

x (float) – Sine value

Returns:

Degree

Return type:

float

seispy.geo.atand(x)

Inverse tangent function with degree as output

Parameters:

x (float) – Tangent value

Returns:

Degree

Return type:

float

seispy.geo.cosd(deg)

Cosine function with degree as input

Parameters:

deg (float) – Degree

Returns:

Cosine value

Return type:

float

seispy.geo.cotd(deg)

Cotangent function with degree as input

Parameters:

deg (float) – Degree

Returns:

Cotangent value

Return type:

float

seispy.geo.deg2km(deg)

Convert degree to km

Parameters:

deg (float) – Distance in degree

Returns:

Distance in km

Return type:

float

seispy.geo.extrema(x, opt='max')

seispy.geo.geo2sph(dep, lat, lon)

Convert geographic coordinates to spherical coordinates.

Parameters:

• dep (float or np.ndarray) – Depth from surface

• lat (float or np.ndarray) – Latitude

• lon (float or np.ndarray) – Longitude

Returns:

spherical coordinate r, theta, phi

Return type:

list

seispy.geo.geoproject(lat_p, lon_p, lat1, lon1, lat2, lon2)

Project a point to a line

Parameters:

• lat_p (float) – Latitude of the point

• lon_p (float) – Longitude of the point

• lat1 (float) – Latitude of the first point of the line

• lon1 (float) – Longitude of the first point of the line

• lat2 (float) – Latitude of the second point of the line

• lon2 (float) – Longitude of the second point of the line

Returns:

Latitude and longitude of the projected point

Return type:

tuple

seispy.geo.km2deg(km)

Convert km to degree

Parameters:

km (float) – Distance in km

Returns:

Distance in degree

Return type:

float

seispy.geo.latlon_from(lat0, lon0, azimuth, gcarc_dist, ellps='WGS84')

Determine position with given position of initial point, azimuth and distance

Accepted numeric scalar or array: - int - float - numpy.floating - numpy.integer - list - tuple - array.array - numpy.ndarray - xarray.DataArray - pandas.Series

Parameters:

• lat0 (float or array) – Latitude of original point

• lon0 (float or array) – Longitude of original point

• azimuth (float or array) – Azimuth(s) in degree

• gcarc_dist (float or array) – Distance(s) between initial and terminus point(s) in degree

• ellps (str, optional) – Ellipsoids supported by pyproj, defaults to “WGS84”

Returns:

scalar or array:

Latitude(s) of terminus point(s)

scalar or array:

Longitude(s) of terminus point(s)

seispy.geo.rad2deg(rad)

Convert radian to degree

Parameters:

rad (float) – Radian

Returns:

Degree

Return type:

float

seispy.geo.rot3D(bazi, inc)

Rotate components from ZRT to LQT

M = [[cos(inc) -sin(inc)*sin(bazi) -sin(inc)*cos(bazi)],

[sin(inc) cos(inc)*sin(bazi) cos(inc)*cos(bazi)], [0 -cos(bazi) sin(bazi)]]

Parameters:

• bazi – back-azimuth of station-event pair

• inc – Incidence angle

Returns:

Coefficient matrix m

Return type:

np.ndarray

seispy.geo.rotateSeisENtoTR(e, n, baz)

Rotate EN to TR

Parameters:

• e (np.ndarray) – East component

• n (np.ndarray) – North component

• baz (float) – Back-azimuth

Returns:

T and R components

Return type:

tuple

seispy.geo.rotateSeisZENtoLQT(z, e, n, bazi, inc)

Rotate ZEN to LQT

Parameters:

• z (np.ndarray) – Vertical component

• e (np.ndarray) – East component

• n (np.ndarray) – North component

• bazi (float) – Back-azimuth

• inc (float) – Incidence angle

Returns:

L, Q and T components

Return type:

tuple

seispy.geo.rssq(x)

Root sum square :param x: Input array :type x: np.ndarray :return: Root sum square :rtype: float

seispy.geo.sdeg2skm(sdeg)

Convert s/degree to s/km

Parameters:

sdeg (float) – s/degree

Returns:

s/km

Return type:

float

seispy.geo.sind(deg)

Sine function with degree as input

Parameters:

deg (float) – Degree

Returns:

Sine value

Return type:

float

seispy.geo.skm2sdeg(skm)

Convert s/km to s/degree

Parameters:

skm (float) – s/km

Returns:

s/degree

Return type:

float

seispy.geo.skm2srad(skm)

Convert s/km to s/rad

Parameters:

skm (float) – s/km

Returns:

s/rad

Return type:

float

seispy.geo.snr(x, y)

Signal-to-noise ratio

Parameters:

• x (np.ndarray) – Signal

• y (np.ndarray) – Noise

Returns:

SNR

Return type:

float

seispy.geo.sph2geo(r, theta, phi)

Convert spherical coordinates to geographic coordinates.

Parameters:

• r (float) – radial distance from coordinate system origin {Units: km, Range: [0, inf)}

• theta (float) – polar angle {Units: radians, Range: [0, π]}

• phi (float) – azimuthal angle {Units: radians, Range: [-π, π]}

Returns:

geographic coordinate conversion (lat, lon, depth) of spherical coordinates

Return type:

(float, float, float)

seispy.geo.spherical2cartesian(lon, lat, dep)

Convert spherical coordinates to cartesian coordinates

Parameters:

• lon (float) – Longitude

• lat (float) – Latitude

• dep (float) – Depth

Returns:

Cartesian coordinates

Return type:

tuple

seispy.geo.srad2skm(srad)

Convert s/rad to s/km

Parameters:

srad (float) – s/rad

Returns:

s/km

Return type:

float

seispy.geo.tand(deg)

Tangent function with degree as input

Parameters:

deg (float) – Degree

Returns:

Tangent value

Return type:

float

seispy.get_cpt module

Created: 2020.03 Copyright: (c) Dimitrios Bouziotas (bouziot) Licence: GGNU General Public License v3 (GPL-3) -You may freely copy, distribute and modify the software, in accordance with the provided license. -You may not sublicense or hold the original author liable. This software comes with no warranty at all. -Active contributions, forks and redevelopments are welcome. -If you would like to include this software in your work, please reference it using the zenodo or github link. Please also reference the original pycpt package (https://github.com/j08lue/pycpt) ——————————————————————————–

seispy.get_cpt.get_cmap(cpt_path, name=None, method='cdict', ret_cmap_type='LinearSegmented', N=None)

Get the cpt colormap as a LinearSegmented colormap. Utilizes the gmtColormap_openfile parser. Parameters ———- cpt_path : str, with one of the following options:

• the full dir path to a .cpt file

• the filename of a .cpt file that is in the local repo (check get_cpt.basedir)

• a url.

namestr, optional

colormap name if not provided, the name will be derived automatically using _getname()

method: str, optional

Choose between ‘cdict’ and ‘list’. The first one fetches all info from the .cpt file. The latter allows you to control the number of colors to sample from. Check gmtColormap_openfile for more info.

N: int, optional

the number of colors in the colormap to be returned. Can define the granularity of the returned colormap. Only useful when method=’list’

seispy.get_cpt.get_listed_cmap(cpt_path, name=None, N=None)

Get the cpt colormap as a ListedColormap. Utilizes the gmtColormap_openfile parser. Parameters ———- cpt_path : str, with one of the following options:

• the full dir path to a .cpt file

• the filename of a .cpt file that is in the local repo (check get_cpt.basedir)

• a url

namestr, optional

colormap name if not provided, the name will be derived automatically using _getname()

N: int, optional

the number of colors in the colormap to be returned. Leave None to derive the colors from the .cpt file. If you use a number less than the colors included in that file, a subset of colors will be returned.

seispy.get_cpt.gmtColormap_openfile(cptf, name=None, method='cdict', N=None, ret_cmap_type='LinearSegmented')

Read a GMT color map from an OPEN cpt file Edited by: bouziot, 2020.03

Parameters:

cptfstr, open file or url handle

path to .cpt file

namestr, optional

name for color map if not provided, the file name will be used

methodstr, suggests the method to use.

If method = ‘cdict’, generates the LinearSegmentedColormap using a color dictionary (cdict), disregarding any value in N.

If method = ‘list’, generates the LinearSegmentedColor using the .from_list() method, passing a list of (value, (r,g,b)) tuples obtained from the cpt file. This allows the selection of colormap resolution by the user, using the N parameter

Nint, the number of colors in the colormap. Only useful when method=’list’.

ret_cmap_type: str, the type of matplotlib cmap object to be returned. Accepts either ‘LinearSegmented’, which returns a matplotlib.colors.LinearSegmentedColormap, or ‘Listed’, which returns a ListedColormap

In case ‘Listed’ is selected, the method argument from the user is ignored and method is set to ‘list’ (‘Linear’ doesn’t work with ‘cdict’).

N is then passed to matplotlib.colors.ListedColormap().

- If N is set to None: all colors of the cpt file will be returned as a list.

- In case of a user-defined N, colors will be truncated or extended by repetition (see matplotlib.colors.ListedColormap).

Returns:

a matplotlib colormap object (matplotlib.colors.LinearSegmentedColormap or matplotlib.colors.ListedColormap)

seispy.get_cpt.plot_cmaps(cmap_list, width=6, cmap_height=0.5, axes_off=False)

Plot a colormap or list of colormaps with their names. Parameters ——- cmap_list (str, cmap object or list of cmap objects anr strings): a list of colormaps to plot, either as cmap objects OR as preinstalled matplotlib colormap strings width (float): width of plot cmap_height (float): height of each colormap in plot axes_off (bool): boolean to erase axes

Returns:

a matplotlib figure object (matplotlib.figure.Figure)

seispy.hk module

seispy.hk.ci(allstack, h, kappa, ev_num)

Search best H and kappa from stacked matrix. Calculate error for H and kappa

Parameters:

• allstack – stacked HK matrix

• h – 1-D array of H

• kappa – 1-D array of kappa

• ev_num – event number

Returns:

seispy.hk.hk()

seispy.hk.hksta(hpara: HKPara, isplot=False, isdisplay=False)

seispy.hk.hkstack(seis, t0, dt, p, h, kappa, vp=6.3, weight=(0.7, 0.2, 0.1))

seispy.hk.plot(stack, allstack, h, kappa, besth, bestk, cvalue, cmap=<matplotlib.colors.ListedColormap object>, title=None, path=None)

seispy.hk.print_result(besth, bestk, maxhsig, maxksig, print_comment=True)

seispy.hk.time2idx(times, ti0, dt)

seispy.hk.tppps(depth, eta_p, eta_s)

seispy.hk.tps(depth, eta_p, eta_s)

seispy.hk.tpsps(depth, eta_s)

seispy.hk.transarray(array, axis=0)

seispy.hk.vslow(v, rayp)

seispy.hkpara module

class seispy.hkpara.HKPara

Bases: object

Attributes:

hrange

krange

property hrange

property krange

seispy.hkpara.hkpara(cfg_file)

seispy.io module

class seispy.io.Query(server='IRIS', **kwargs)

Bases: object

Methods

get_events([starttime, endtime])

Get events from IRIS

get_stations

get_events(starttime=None, endtime=UTCDateTime(2024, 8, 13, 7, 38, 28, 336560), **kwargs)

Get events from IRIS

Parameters:

• starttime (obspy.UTCDateTime, optional) – Start time of events, defaults to None

• endtime (obspy.UTCDateTime, optional) – End time of events, defaults to UTCDateTime.now()

Returns:

Events

Return type:

obspy.Catalog

get_stations(includerestricted=False, **kwargs)

seispy.mccc module

seispy.mccc.mccc(data, dt, twin=0)

seispy.para module

class seispy.para.RFPara

Bases: object

Attributes:

catalogpath

criterion

datapath

decon_method

rfpath

server

Methods

read_para(cfg_file)

Read parameters from configure file :param cfg_file: Path to configure file :type cfg_file: str

get_para

property catalogpath

property criterion

property datapath

property decon_method

get_para()

classmethod read_para(cfg_file)

Read parameters from configure file :param cfg_file: Path to configure file :type cfg_file: str

property rfpath

property server

class seispy.para.StaInfo

Bases: object

Methods

get_stainfo()	Get station information
get_station_from_ws(**kwargs)	Get station information from web-service with given network and station or other optional condition.
link_server(server[, user, password])	_summary_
load_stainfo(pathname, ref_comp, suffix)	Load station information from SAC file

get_stainfo()

Get station information

Returns:

Station information

Return type:

dict

get_station_from_ws(**kwargs)

Get station information from web-service with given network and station or other optional condition.

link_server(server, user=None, password=None)

_summary_

Parameters:

• server (str) – Sever name of FDSN web-service, by default ‘IRIS’

• user (str, optional) – username, defaults to None

• password (str, optional) – password, defaults to None

load_stainfo(pathname, ref_comp, suffix)

Load station information from SAC file

Parameters:

• pathname (str) – Path to SAC file

• ref_comp (str) – Reference component

• suffix (str) – Suffix of SAC file

seispy.plotR module

class seispy.plotR.PlotR(stadata, enf=12, xlim=[2, 80], key='bazi')

Bases: object

Methods

init_figure()	Initialize figure
plot([out_path, outformat, show])	Plot receiver function in R component :param out_path: output path :type out_path: str :param outformat: output format :type outformat: str :param show: show figure :type show: bool
plot_waves()	Plot PRFs with R component
set_fig()	Set figure

init_figure()

Initialize figure

plot(out_path=None, outformat='g', show=False)

Plot receiver function in R component :param out_path: output path :type out_path: str :param outformat: output format :type outformat: str :param show: show figure :type show: bool

plot_waves()

Plot PRFs with R component

set_fig()

Set figure

seispy.plotR.plotr(rfsta, out_path='./', xlim=[-2, 80], key='bazi', enf=6, outformat='g', show=False)

Plot receiver function in R component :param rfsta: Path to PRFs :type rfsta: seispy.rfcorrect.RFStation :param enf: The enlarge factor, defaults to 6 :type enf: int, optional :param out_path: The output path, defaults to current directory :type out_path: str, optional :param key: The key to sort PRFs, available for event, evla, evlo, evdp,

dis, bazi, rayp, mag, f0, defaults to bazi

Parameters:

• outformat (str, optional) – File format of the image file, g as ‘png’, f as ‘pdf’, defaults to ‘g’

• xlim (list, optional) – xlim, defaults to [-2, 80]

• show (bool) – show figure

Returns:

PlotR object

Return type:

PlotR

seispy.plotRT module

class seispy.plotRT.PlotRT(stadata, enf=3, xlim=[2, 80], key='bazi')

Bases: object

Methods

plot([out_path, outformat, show])	Plot PRFs with R and T components
plot_waves()	Plot PRFs with R and T components
set_fig()	Set figure

init_figure

init_figure()

plot(out_path=None, outformat='g', show=False)

Plot PRFs with R and T components

Parameters:

• rfsta (seispy.rfcorrect.RFStation) – Path to PRFs

• enf (int, optional) – The enlarge factor, defaults to 3

• out_path (str, optional) – The output path, defaults to current directory

• key (str, optional) – The key to sort PRFs, avialible for event, evla, evlo, evdp, dis, bazi, rayp, mag, f0, defaults to bazi

• outformat (str, optional) – File format of the image file, g as ‘png’, f as ‘pdf’, defaults to ‘g’

plot_waves()

Plot PRFs with R and T components

set_fig()

Set figure

seispy.plotRT.plotrt(rfsta, out_path='./', xlim=[-2, 30], key='bazi', enf=6, outformat='g', show=False)

Plot PRFs with R and T components

Parameters:

• rfsta (seispy.rfcorrect.RFStation) – Path to PRFs

• enf (int, optional) – The enlarge factor, defaults to 6

• out_path (str, optional) – The output path, defaults to current directory

• key (str, optional) – The key to sort PRFs, avialible for event, evla, evlo, evdp, dis, bazi, rayp, mag, f0, defaults to bazi

• outformat (str, optional) – File format of the image file, g as ‘png’, f as ‘pdf’, defaults to ‘g’

• xlim (list, optional) – xlim, defaults to [-2, 30]

• show (bool) – show figure

seispy.psrayp module

class seispy.psrayp.PsRayp(dis, dep, laymin=0, laymax=800)

Bases: object

Methods

taup_rayp([this_dis, this_dep, taup])

get_rayp
make_phase_list
save

get_rayp()

make_phase_list()

save(path='Ps_rayp')

taup_rayp(this_dis=50, this_dep=10, taup='taup_time')

Parameters:

taup – Path to taup_time

Returns:

seispy.psrayp.gen_rayp_lib()

seispy.psrayp.get_psrayp(rayp_lib, dis, dep, layers)

seispy.rf module

seispy.rf.CfgModify(cfg_file, session, key, value)

class seispy.rf.RF(cfg_file=None, log=None)

Bases: object

Attributes:

date_begin

date_end

Methods

baz_correct([time_b, time_e, offset, ...])	Correct back-azimuth for all data
cal_phase()	Calculate arrivals and ray parameters for all data
channel_correct()	Correct channel components
deconv()	Deconvolution for all data
detrend()	Detrend all data
drop_eq_snr([length, z_only])	Drop earthquakes with low SNR
filter([freqmin, freqmax, order])	Filter all data
load_stainfo([use_date_range])	Load station information from local file or remote web-service
loadpjt(path)	Load project from local disk
match_eq()	Assosiate earthquakes with local data file or online data.
pick([prepick, stl, ltl])	Pick phase arrival
rotate([search_inc])	Rotate all data to ZNE or RTZ
save_raw_data()	Save raw data to local disk
savepjt()	Save project to local disk
saverf([gauss])	Save receiver functions to local disk
search_eq([local, catalog])	Search earthquakes from local or online data server
trim()	Trim waveforms from start to end

baz_correct(time_b=10, time_e=20, offset=90, correct_angle=None)

Correct back-azimuth for all data

Parameters:

• time_b (int, optional) – Begin time of searching, defaults to 10

• time_e (int, optional) – End time of searching, defaults to 20

• offset (int, optional) – Offset of searching, defaults to 90

• correct_angle (float, optional) – Correct angle, defaults to None

cal_phase()

Calculate arrivals and ray parameters for all data

channel_correct()

Correct channel components

property date_begin

property date_end

deconv()

Deconvolution for all data

detrend()

Detrend all data

drop_eq_snr(length=None, z_only=False)

Drop earthquakes with low SNR

Parameters:

• length (int, optional) – Length of data, defaults to None

• z_only (bool, optional) – Use Z component only, defaults to False

filter(freqmin=None, freqmax=None, order=4)

Filter all data

Parameters:

• freqmin (float, optional) – Minimum frequency, defaults to self.para.freqmin

• freqmax (float, optional) – Maximum frequency, defaults to self.para.freqmax

• order (int, optional) – Order of filter, defaults to 4

load_stainfo(use_date_range=True)

Load station information from local file or remote web-service

Parameters:

use_date_range (bool, optional) – Use date range to search station information, defaults to True

classmethod loadpjt(path)

Load project from local disk

Parameters:

path (str) – Path to project file

Returns:

Project object

Return type:

seispy.rf.RF

match_eq()

Assosiate earthquakes with local data file or online data.

pick(prepick=True, stl=5, ltl=10)

Pick phase arrival

Parameters:

• prepick (bool, optional) – Use STA/LTA method, defaults to True

• stl (int, optional) – Short time length, defaults to 5

• ltl (int, optional) – Long time length, defaults to 10

rotate(search_inc=False)

Rotate all data to ZNE or RTZ

Parameters:

search_inc (bool, optional) – Search incidence angle, defaults to False

save_raw_data()

Save raw data to local disk

savepjt()

Save project to local disk

saverf(gauss=None)

Save receiver functions to local disk

Parameters:

gauss (float, optional) – Gaussian width, defaults to self.para.gauss

search_eq(local=False, catalog=None)

Search earthquakes from local or online data server

Parameters:

• local (bool, optional) – Search from local data, defaults to False

• catalog (str, optional) – Catalog type, defaults to None

trim()

Trim waveforms from start to end

exception seispy.rf.SACFileNotFoundError(matchkey)

Bases: Exception

seispy.rf.datestr2regex(datestr)

seispy.rf.fetch_waveform(eq_lst, para, model, logger)

Fetch waveforms from remote data server

Parameters:

• eq_lst (pandas.DataFrame) – Earthquake list

• para (seispy.para.RFPara) – RFPara object

• model (obspy.taup.TauPyModel) – TauPyModel object

• logger (seispy.setuplog.SetupLog) – Logger

Returns:

Earthquake list with fetched waveforms

Return type:

pandas.DataFrame

seispy.rf.match_eq(eq_lst, pathname, stla, stlo, logger, ref_comp='Z', suffix='SAC', offset=None, tolerance=210, dateformat='%Y.%j.%H.%M.%S')

Match earthquakes with local SAC files

Parameters:

• eq_lst (pandas.DataFrame) – Earthquake list

• pathname (str) – Path to SAC files

• stla (float) – Station latitude

• stlo (float) – Station longitude

• logger (seispy.setuplog.SetupLog) – Logger

• ref_comp (str, optional) – Reference component, defaults to ‘Z’

• suffix (str, optional) – Suffix of SAC files, defaults to ‘SAC’

• offset (float, optional) – Time offset between SAC files and earthquakes, defaults to None

• tolerance (int, optional) – Tolerance of time offset, defaults to 210

• dateformat (str, optional) – Date format of SAC files, defaults to ‘%Y.%j.%H.%M.%S’

Returns:

Earthquake list with matched SAC files

Return type:

pandas.DataFrame

seispy.rf.pickphase(eqs, para, logger)

seispy.rf.process_row(i, size, row, para, model, query, tb, te, logger)

seispy.rf.read_catalog(logpath: str, b_time, e_time, stla: float, stlo: float, magmin=5.5, magmax=10.0, dismin=30.0, dismax=90.0, depthmin=0, depthmax=800)

Read local catalog with seispy or QUAKEML format

Parameters:

• logpath (str) – Path to catalogs

• b_time (obspy.UTCDateTime) – Start time

• e_time (obspy.UTCDateTime) – End time

• stla (float) – Station latitude

• stlo (float) – Station longitude

• magmin (float, optional) – Minimum magnitude, defaults to 5.5

• magmax (float, optional) – Maximum magnitude, defaults to 10

• dismin (float, optional) – Minimum distance, defaults to 30

• dismax (float, optional) – Maximum distance, defaults to 90

Returns:

list of earthquakes

Return type:

pandas.DataFrame

seispy.rf.setpar()

Set parameters to configure file

seispy.rf2depth_makedata module

class RF2depth : process cal RF2depth class sta_part, sta_full, _RFInd

class seispy.rf2depth_makedata.RFDepth(cpara: ~seispy.ccppara.CCPPara, log: ~logging.Logger = <Logger RF2depth (INFO)>, raytracing3d=False, velmod3d=None, modfolder1d=None)

Bases: object

Convert receiver function to depth axis

Methods

makedata([psphase])

Convert receiver function to depth axis

makedata(psphase=1)

Convert receiver function to depth axis

Parameters:

psphase (int) – 1 for Ps, 2 for PpPs, 3 for PpSs

class seispy.rf2depth_makedata.Station(sta_lst: str)

Bases: object

seispy.rf2depth_makedata.rf2depth()

CLI for Convert receiver function to depth axis There’s 4 branch provided to do RF 2 depth conversion

1. only -d :do moveout correction

2. only -r : do raytracing but no moveout correction

3. -d and -r : do moveout correction and raytracing

4. -m : use {staname}.vel file for RF2depth conversion

class seispy.rf2depth_makedata.sta_full(station, stla, stlo, stel)

Bases: tuple

Methods

count(value, /)	Return number of occurrences of value.
index(value[, start, stop])	Return first index of value.

station

Alias for field number 0

stel

Alias for field number 3

stla

Alias for field number 1

stlo

Alias for field number 2

class seispy.rf2depth_makedata.sta_part(station, stla, stlo)

Bases: tuple

Methods

count(value, /)	Return number of occurrences of value.
index(value[, start, stop])	Return first index of value.

station

Alias for field number 0

stla

Alias for field number 1

stlo

Alias for field number 2

seispy.rfani module

class seispy.rfani.RFAni(sacdatar, tb, te, tlen=3, val=10, rayp=0.06, model='iasp91')

Bases: object

Methods

baz_stack([val])	Stack RF data with back-azimuth.
joint_ani([weight])	Joint method for crustal anisotropy estimation.
plot_correct([fvd, dt, enf, outpath])	Plot the RF data with back-azimuth and time after P corrected.
plot_polar([cmap, show, outpath])	Polar map of crustal anisotropy.
plot_stack_baz([enf, outpath])	Plot the stack of RF data with back-azimuth.
radial_energy_max()	Calculate the energy of radial component.
rotate_to_fast_slow()	Rotate RF data to fast and slow direction.
search_peak(energy[, opt])	Search the peak of energy.
search_peak_list(energy[, opt])	Search the peak of energy.
xyz2grd(energy)	Interpolate energy to grid.

baz_stack(val=10)

Stack RF data with back-azimuth.

Parameters:

valint, optional

Interval of back-azimuth, by default 10

joint_ani(weight=[0.5, 0.3, 0.2])

Joint method for crustal anisotropy estimation.

Parameters:

weightlist, optional

Weight for three energy matrix, by default [0.5, 0.3, 0.2]

plot_correct(fvd=0, dt=0.44, enf=80, outpath=None)

Plot the RF data with back-azimuth and time after P corrected.

Parameters:

fvdint, optional

Fast velocity direction, by default 0

dtfloat, optional

Time delay for correction, by default 0.44

enfint, optional

Enlarge factor for back-azimuth, by default 80

outpathstr, optional

Output path for saving the figure, by default None

plot_polar(cmap=<matplotlib.colors.ListedColormap object>, show=False, outpath='./')

Polar map of crustal anisotropy. See Liu and Niu (2012, doi: 10.1111/j.1365-246X.2011.05249.x) in detail.

Parameters:

• cmap (optional) – Colormap of matplotlib, defaults to ‘rainbow’

• show (bool, optional) – If show the polar map in the Matplotlib window, defaults to True

• outpath (str, optional) – Output path to saving the figure. If show the figure in the Matplotlib window, this option will be invalid, defaults to current directory.

plot_stack_baz(enf=60, outpath='./')

Plot the stack of RF data with back-azimuth.

Parameters:

enfint, optional

Enlarge factor for back-azimuth, by default 60

outpathstr, optional

Output path for saving the figure, by default ‘./’

radial_energy_max()

Calculate the energy of radial component.

rotate_to_fast_slow()

Rotate RF data to fast and slow direction.

search_peak(energy, opt='max')

Search the peak of energy.

Parameters:

energynp.ndarray

Energy matrix

optstr, optional

Option for searching peak, by default ‘max’

search_peak_list(energy, opt='max')

Search the peak of energy.

Parameters:

energynp.ndarray

Energy matrix

optstr, optional

Option for searching peak, by default ‘max’

xyz2grd(energy)

Interpolate energy to grid.

seispy.rfcorrect module

class seispy.rfcorrect.RFStation(data_path, only_r=False, prime_comp='R')

Bases: object

Attributes:

stel

Methods

bin_stack([key, lim, val])	Stack RFs by bins of key with interval of val
harmonic([tb, te, is_stack])	Harmonic decomposition for extracting anisotropic and isotropic features from the radial and transverse RFs
jointani(tb, te[, tlen, stack_baz_val, ...])	Eastimate crustal anisotropy with a joint method.
moveoutcorrect([ref_rayp, dep_range, ...])	Moveout correction with specified reference ray-parameter and depth
normalize([method])	Normalize amplitude of each RFs.
plotr([outformat])	Plot radial RFs
plotrt([outformat])	Plot radial and transverse RFs
psrf2depth([dep_range])	Time-to-depth conversion with specified depth series.
psrf_1D_raytracing([dep_range])	1D back ray tracing to obtained Ps conversion points at discret depths
psrf_3D_moveoutcorrect(mod3dpath, **kwargs)	3D moveout correction with 3D velocity model
psrf_3D_raytracing(mod3dpath[, dep_range, srayp])	3D back ray tracing to obtained Ps conversion points at discret depths
psrf_3D_timecorrect(mod3dpath[, dep_range, ...])	3D time-to-depth conversion with 3D velocity model
read_stream(stream, rayp, baz[, prime_comp, ...])	Create RFStation instance from obspy.Stream
resample(dt)	Resample RFs with specified dt
slantstack([ref_dis, rayp_range, tau_range])	Slant stack for receiver function
sort([key])	Sort RFs by keys in given event, evla, evlo, evdp, dis, bazi, rayp, mag, f0

bin_stack(key='bazi', lim=[0, 360], val=10)

Stack RFs by bins of key with interval of val

Parameters:

• key (str, optional) – Key to stack, valid in [‘bazi’, ‘rayp’ (in s/km)], defaults to ‘bazi’

• val (int, optional) – Interval of bins, defaults to 10 degree for bazi

Returns:

Stacked RFs

Return type:

dict with keys of data_prime, datat and count

harmonic(tb=-5, te=10, is_stack=True)

Harmonic decomposition for extracting anisotropic and isotropic features from the radial and transverse RFs

Parameters:

• tb (float, optional) – Start time relative to P, defaults to -5

• te (float, optional) – End time relative to P, defaults to 10

• is_stack (bool, optional) – Wether stack the result, defaults to True

Returns:

Harmonic components and unmodel components

Return type:

numpy.ndarray, numpy.ndarray

jointani(tb, te, tlen=3.0, stack_baz_val=10, rayp=0.06, velmodel='iasp91', weight=[0.4, 0.4, 0.2])

Eastimate crustal anisotropy with a joint method. See Liu and Niu (2012, doi: 10.1111/j.1365-246X.2011.05249.x) in detail.

Parameters:

• tb (float) – Time before Pms for search Ps peak

• te (float) – Time after Pms for search Ps peak

• tlen (float, optional) – Half time length for cut out Ps phase, defaults to 3.0

• stack_baz_val (float, optional) – The interval for stacking binned by back-azimuth, defaults to 10

• rayp (float, optional) – Reference ray-parameter for moveout correction, defaults to 0.06

• velmodel (str, optional) – velocity model for moveout correction. ‘iasp91’, ‘prem’ and ‘ak135’ is valid for internal model. Specify path to velocity model for the customized model. The format is the same as in Taup, but the depth should be monotonically increasing, defaults to ‘iasp91’

• weight (list, optional) – Weight for three different method, defaults to [0.4, 0.4, 0.2]

Returns:

Dominant fast velocity direction and time delay

Return type:

List, List

moveoutcorrect(ref_rayp=0.06, dep_range=array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149]), velmod='iasp91', replace=False, **kwargs)

Moveout correction with specified reference ray-parameter and depth

Parameters:

• ref_rayp (float, optional) – reference ray-parameter in s/km, defaults to 0.06

• dep_range (np.ndarray, optional) – Depth range used for extracting velocity in velocity model, defaults to np.arange(0, 150)

• velmod (str, optional) – Velocity model for moveout correction. ‘iasp91’, ‘prem’ and ‘ak135’ is valid for internal model. Specify path to velocity model for the customized model. The format is the same as in Taup, but the depth should be monotonically increasing, defaults to ‘iasp91’

• replace (bool, optional) – whether replace original data, False to return new array, defaults to False

Returns:

rf_corr: np.ndarray

Corrected RFs with component of RFStation.comp

t_corr: np.ndarray or None

Corrected RFs in transverse component. If only_r is True, this variable is None

normalize(method='single')

Normalize amplitude of each RFs.

Parameters:

method (str, optional) – Method of normalization with single and average avaliable. - single for normalization with max amplitude of current RF. - average for normalization with average amplitude of current station.

plotr(outformat=None, **kwargs)

Plot radial RFs

Parameters:

outformat (str, optional) – Output format for plot, defaults to None

plotrt(outformat=None, **kwargs)

Plot radial and transverse RFs

Parameters:

outformat (str, optional) – Output format for plot, defaults to None

psrf2depth(dep_range=array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149]), **kwargs)

Time-to-depth conversion with specified depth series.

Parameters:

• dep_range (np.ndarray() , optional) – Discret conversion depth, defaults to np.arange(0, 150)

• velmod (str, optional) – Velocity model for time-to-depth conversion. ‘iasp91’, ‘prem’ and ‘ak135’ is valid for internal model. Specify path to velocity model for the customized model. The format is the same as in Taup, but the depth should be monotonically increasing, defaults to ‘iasp91’

• srayp (numpy.lib.npyio.NpzFile, optional) – Ray-parameter lib for Ps phases, If set up to None the rayp of direct is used, defaults to None

Returns:

2D array of RFs in depth

Return type:

np.ndarray()

psrf_1D_raytracing(dep_range=array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149]), **kwargs)

1D back ray tracing to obtained Ps conversion points at discret depths

Parameters:

• dep_range (numpy.ndarray, optional) – Discret conversion depth, defaults to np.arange(0, 150)

• velmod (str, optional) – Velocity model for time-to-depth conversion. ‘iasp91’, ‘prem’ and ‘ak135’ is valid for internal model. Specify path to velocity model for the customized model. The format is the same as in Taup, but the depth should be monotonically increasing, defaults to ‘iasp91’

• srayp (numpy.lib.npyio.NpzFile, optional) – Ray-parameter lib for Ps phases, If set up to None the rayp of direct is used, defaults to None

Return pplat_s:

Latitude of conversion points

Return pplon_s:

Longitude of conversion points

Return tps:

Time difference of Ps at each depth

Return type:

list

psrf_3D_moveoutcorrect(mod3dpath, **kwargs)

3D moveout correction with 3D velocity model

Parameters:

• mod3dpath (str) – Path to 3D velocity model

• dep_range (numpy.ndarray, optional) – Discret conversion depth, defaults to np.arange(0, 150)

• srayp (numpy.lib.npyio.NpzFile, optional) – Ray-parameter lib for Ps phases, If set up to None the rayp of direct is used, defaults to None

Returns:

2D array of RFs in depth

Return type:

numpy.ndarray

psrf_3D_raytracing(mod3dpath, dep_range=array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149]), srayp=None)

3D back ray tracing to obtained Ps conversion points at discret depths

Parameters:

• mod3dpath (str) – Path to 3D velocity model

• dep_range (numpy.ndarray, optional) – Discret conversion depth, defaults to np.arange(0, 150)

• srayp (numpy.lib.npyio.NpzFile, optional) – Ray-parameter lib for Ps phases, If set up to None the rayp of direct is used, defaults to None

Return pplat_s:

Latitude of conversion points

Return pplon_s:

Longitude of conversion points

Return tps:

Time difference of Ps at each depth

Return type:

list

psrf_3D_timecorrect(mod3dpath, dep_range=array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149]), normalize='single', **kwargs)

3D time-to-depth conversion with 3D velocity model

Parameters:

• mod3dpath (str) – Path to 3D velocity model

• dep_range (numpy.ndarray, optional) – Discret conversion depth, defaults to np.arange(0, 150)

• normalize (str, optional) – Normalization method, defaults to ‘single’, see RFStation.normalize for detail

Returns:

2D array of RFs in depth

Return type:

numpy.ndarray

classmethod read_stream(stream, rayp, baz, prime_comp='R', stream_t=None)

Create RFStation instance from obspy.Stream

Parameters:

• stream (obspy.Stream()) – Stream of RFs

• rayp (float or np.ndarray()) – Ray-parameter of RFs

• baz (float or np.ndarray()) – Back-azimuth of RFs

• prime_comp (str, optional) – Prime component in RF filename. R or Q for PRF and L or Z for SRF, defaults to ‘R’

• stream_t (obspy.Stream(), optional) – Stream of transverse RFs, defaults to None

Returns:

RFStation instance

Return type:

RFStation

resample(dt)

Resample RFs with specified dt

Parameters:

dt (float) – New sampling rate

slantstack(ref_dis=None, rayp_range=None, tau_range=None)

Slant stack for receiver function

Parameters:

• ref_dis (int or float, optional) – reference distance, by default None

• rayp_range (numpy.ndarray, optional) – range of ray parameter, by default None

• tau_range (numpy.ndarray, optional) – range of tau, by default None

sort(key='bazi')

Sort RFs by keys in given event, evla, evlo, evdp, dis, bazi, rayp, mag, f0

Parameters:

key (str, optional) – key to sort, defaults to bazi

property stel

class seispy.rfcorrect.SACStation(data_path, only_r=False)

Bases: RFStation

Attributes:

stel

Methods

bin_stack([key, lim, val])	Stack RFs by bins of key with interval of val
harmonic([tb, te, is_stack])	Harmonic decomposition for extracting anisotropic and isotropic features from the radial and transverse RFs
jointani(tb, te[, tlen, stack_baz_val, ...])	Eastimate crustal anisotropy with a joint method.
moveoutcorrect([ref_rayp, dep_range, ...])	Moveout correction with specified reference ray-parameter and depth
normalize([method])	Normalize amplitude of each RFs.
plotr([outformat])	Plot radial RFs
plotrt([outformat])	Plot radial and transverse RFs
psrf2depth([dep_range])	Time-to-depth conversion with specified depth series.
psrf_1D_raytracing([dep_range])	1D back ray tracing to obtained Ps conversion points at discret depths
psrf_3D_moveoutcorrect(mod3dpath, **kwargs)	3D moveout correction with 3D velocity model
psrf_3D_raytracing(mod3dpath[, dep_range, srayp])	3D back ray tracing to obtained Ps conversion points at discret depths
psrf_3D_timecorrect(mod3dpath[, dep_range, ...])	3D time-to-depth conversion with 3D velocity model
read_stream(stream, rayp, baz[, prime_comp, ...])	Create RFStation instance from obspy.Stream
resample(dt)	Resample RFs with specified dt
slantstack([ref_dis, rayp_range, tau_range])	Slant stack for receiver function
sort([key])	Sort RFs by keys in given event, evla, evlo, evdp, dis, bazi, rayp, mag, f0

seispy.rfcorrect.moveoutcorrect_ref(stadatar, raypref, YAxisRange, chan='', velmod='iasp91', sphere=True, phase=1)

Moveout correction refer to a specified ray-parameter

Parameters:

• stadatar – data class of RFStation

• raypref – referred ray parameter in rad

• YAxisRange – Depth range in nd.array type

• velmod – Path to velocity model

• chan – channel name for correction, ‘r’, ‘t’…

Returns:

Newdatar, EndIndex

seispy.rfcorrect.psrf2depth(stadatar, YAxisRange, velmod='iasp91', srayp=None, normalize='single', sphere=True, phase=1)

Time-to-depth conversion with S-wave backprojection.

Parameters:

• stadatar (RFStation()) – Data class of RFStation

• YAxisRange (numpy.ndarray) – Depth range for conversion

• velmod (str, optional) – Velocity for conversion, whcih can be a path to velocity file, defaults to ‘iasp91’

• srayp (str or seispy.psrayp.PsRayp(), optional) – ray-parameter library of conversion phases. See seispy.psrayp() in detail, defaults to None

• normalize (str, optional) – method of normalization, defaults to ‘single’. Please refer to RFStation.normalize()

• sphere (bool, optional) – Wether do earth-flattening transformation, defaults to True

Returns:

ps_rfdepth: 2-D numpy.ndarray, float

RFs in depth with shape of (stadatar.ev_num, YAxisRange.size), stadatar.ev_num is the number of RFs in current station. YAxisRange.size is the size of depth axis.

endindex: numpy.ndarray, int

End index of each RF in depth

x_s: 2-D numpy.ndarray, float

Horizontal distance between station and S-wave conversion points with shape of (stadatar.ev_num, YAxisRange.size)

x_p: 2-D numpy.ndarray, float

Horizontal distance between station and P-wave conversion points with shape of (stadatar.ev_num, YAxisRange.size)

seispy.rfcorrect.psrf_1D_raytracing(stadatar, YAxisRange, velmod='iasp91', srayp=None, sphere=True, phase=1)

seispy.rfcorrect.psrf_3D_migration(pplat_s, pplon_s, pplat_p, pplon_p, raylength_s, raylength_p, Tpds, dep_range, mod3d)

3D time difference correction with specified ray path and 3D velocity model.

The input parameters can be generated with psrf_1D_raytracing().

Parameters:

pplat_snp.ndarray()

2D array of latitude of S-wave in dep_range, (RFStation.ev_num(), dep_range.size)

pplon_snp.ndarray()

2D array of longitude of S-wave in dep_range, (RFStation.ev_num(), dep_range.size)

pplat_pnp.ndarray()

2D array of latitude of P-wave in dep_range, (RFStation.ev_num(), dep_range.size)

pplon_pnp.ndarray()

2D array of longitude of P-wave in dep_range, (RFStation.ev_num(), dep_range.size)

raylength_snp.ndarray()

2D array of ray path length of S-wave in dep_range, (RFStation.ev_num(), dep_range.size)

raylength_pnp.ndarray()

2D array of ray path length of P-wave in dep_range, (RFStation.ev_num(), dep_range.size)

Tpdsnp.ndarray()

1D array of time difference in dep_range (dep_range.size)

dep_rangenp.ndarray()

1D array of depths in km, (dep_range.size)

mod3dnp.lib.npyio.NpzFile()

3D velocity loaded from a .npz file

Returns:

np.ndarray()

Corrected time difference in dep_range

seispy.rfcorrect.psrf_3D_raytracing(stadatar, YAxisRange, mod3d, srayp=None, elevation=0, sphere=True)

Back ray trace the S wavs with a assumed ray parameter of P.

Parameters:

• stadatar (object RFStation) – The data class including PRFs and more parameters

• YAxisRange (numpy.ndarray) – The depth array with the same intervals

• mod3d ('Mod3DPerturbation' object) – The 3D velocity model with fields of dep, lat, lon, vp and vs.

• elevation (float) – Elevation of this station relative to sea level

Returns:

pplat_s, pplon_s, pplat_p, pplon_p, tps

Type:

numpy.ndarray * 5

seispy.rfcorrect.time2depth(stadatar, dep_range, Tpds, normalize='single')

Interpolate RF amplitude with specified time difference and depth range

Parameters:

stadatarRFStation()

Data class of RFStation()

dep_rangenp.ndarray()

1D array of depths in km, (dep_range.size)

Tpdsnp.ndarray()

1D array of time difference in dep_range (dep_range.size)

normalizestr, optional

Normlization option, 'sinlge' and 'average' are available , by default ‘single’See RFStation.normalize() in detail.

Returns:

PS_RFdepth: np.ndarray()

2D array of RFs in depth, (RFStation.ev_num(), dep_range.size)

end_index: np.ndarray()

1D array of the last digit of effective depth (RFStation.ev_num())

seispy.rfcorrect.xps_tps_map(dep_mod: DepModel, srayp, prayp, is_raylen=False, sphere=True, phase=1)

Calculate horizontal distance and time difference at depths

Parameters:

• dep_mod (seispy.util.DepModel()) – 1D velocity model class

• srayp (float) – conversion phase ray-parameters

• prayp (float) – S-wave ray-parameters

• is_raylen (bool, optional) – Wether calculate ray length at depths, defaults to False

• sphere (bool, optional) – Wether do earth-flattening transformation, defaults to True, defaults to True

• phase (int, optional) – Phases to calculate 1 for Ps, 2 for PpPs, 3 for PsPs+PpSs, defaults to 1

Returns:

If is_raylen = False

tps: 2-D numpy.ndarray, float

RFs in depth with shape of (stadatar.ev_num, YAxisRange.size), stadatar.ev_num is the number of RFs in current station. YAxisRange.size is the size of depth axis.

x_s: 2-D numpy.ndarray, float

Horizontal distance between station and S-wave conversion points with shape of (stadatar.ev_num, YAxisRange.size)

x_p: 2-D numpy.ndarray, float

Horizontal distance between station and P-wave conversion points with shape of (stadatar.ev_num, YAxisRange.size)

otherwise, two more variables will be returned

raylength_s: 2-D numpy.ndarray, float

raylength_p: 2-D numpy.ndarray, float

seispy.setuplog module

class seispy.setuplog.SetupLog(filename='/home/runner/.RF.log')

Bases: object

default_level = 20

default_logs = {'Batlog': ('Bat', 20, 'stream_handler', 'file_handler'), 'CCPlog': ('CCP', 20, 'stream_handler'), 'ModCreatorlog': ('ModCreator', 20, 'stream_handler'), 'PickDepthlog': ('PickDepth', 20, 'stream_handler'), 'RF2depthlog': ('RF2depth', 20, 'stream_handler'), 'RFlog': ('RF', 20, 'stream_handler')}

seispy.signal module

seispy.signal.smooth(x, half_len=5, window='flat')

smooth the data using a window with requested size.

This method is based on the convolution of a scaled window with the signal. The signal is prepared by introducing reflected copies of the signal (with the window size) in both ends so that transient parts are minimized in the begining and end part of the output signal.

input:

x: the input signal helf_len: the half dimension of the smoothing window; should be an odd integer window: the type of window from ‘flat’, ‘hanning’, ‘hamming’, ‘bartlett’, ‘blackman’

flat window will produce a moving average smoothing.

output:

the smoothed signal

example:

t=linspace(-2,2,0.1) x=sin(t)+randn(len(t))*0.1 y=smooth(x)

see also:

numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve scipy.signal.lfilter

TODO: the window parameter could be the window itself if an array instead of a string NOTE: length(output) != length(input), to correct this: return y[(window_len/2-1):-(window_len/2)] instead of just y.

seispy.signal.whiten(data, Nfft, delta, f1, f2, f3, f4)

This function takes 1-dimensional data timeseries array, goes to frequency domain using fft, whitens the amplitude of the spectrum in frequency domain between freqmin and freqmax and returns the whitened fft.

Parameters:

• data (numpy.ndarray) – Contains the 1D time series to whiten

• Nfft (int) – The number of points to compute the FFT

• delta (float) – The sampling frequency of the data

• freqmin (float) – The lower frequency bound

• freqmax (float) – The upper frequency bound

• plot (bool) – Whether to show a raw plot of the action (default: False)

Return type:

numpy.ndarray

Returns:

The FFT of the input trace, whitened between the frequency bounds

Module contents