seispy package

seispy.ccpprofile module

class seispy.ccpprofile.CCPProfile(cfg_file=None, log=None)

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

load_para

initial_profile()

Initialize bins of profile

load_para(cfg_file)

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. where 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.bin_shape(cpara)

seispy.ccpprofile.fix_filename(filename, typ='dat')

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.ccpprofile.line_proj(lat1, lon1, lat2, lon2)

seispy.cc3d module

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

Bases: object

Methods

save_stack_data(fname)	Save stacked data and parameters to local as a npz file.
stack()	Search conversion points falling within a bin and stack them with bootstrap method.

get_depth_err
initial_grid
load_para
read_rfdep
read_stack_data
save_good_410_660
search_good_410_660

get_depth_err(type='std')

initial_grid()

load_para(cfg_file)

read_rfdep()

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

save_good_410_660(fname)

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)

stack()

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

seispy.ccp3d.bin_shape(cpara)

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

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)

seispy.decon module

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

Bases: 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.
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.

deconvolute

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

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

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

Created on Wed Sep 10 14:21:38 2014

In: uin = numerator (radial for PdS) win = denominator (vertical component for PdS) dt = sample interval (s) nt = number of samples tshift = Time until beginning of receiver function (s) f0 = width of gaussian filter itmax = max # iterations minderr = Min change in error required for stopping iterations

Out: RFI = receiver function rms = Root mean square error for predicting numerator after each iteration

@author: Mijian Xu @ NJU

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

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)

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

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

seispy.distaz module

seispy.distaz.acosd(x)

seispy.distaz.asind(x)

seispy.distaz.atand(x)

seispy.distaz.cosd(deg)

seispy.distaz.cotd(deg)

seispy.distaz.deg2km(degree)

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

Bases: object

c Subroutine to calculate the Great Circle Arc distance c between two sets of geographic coordinates c c Equations take from Bullen, pages 154, 155 c c T. Owens, September 19, 1991 c Sept. 25 – fixed az and baz calculations c 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 availible input varible.

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.distaz.km2deg(kilometers)

seispy.distaz.sind(deg)

seispy.distaz.tand(deg)

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])	_summary_
deconvolute(shift, time_after[, f0, method, ...])	Deconvolution
trim(time_before, time_after[, isreturn])	offset = sac.b - real o

cleanstream
decon_p
decon_s
detrend
filter
fix_channel_name
from_stream
get_arrival
get_time_offset
judge_rf
phase_trigger
readstream
rotate
s_condition
saverf
search_baz
search_inc
set_comp
snr
write

arr_correct(write_to_sac=True)

offset = sac.b - real o

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

_summary_

Parameters:

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

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

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

cleanstream()

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

decon_s(tshift, **kwargs)

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()

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

fix_channel_name()

classmethod from_stream(stream)

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

get_time_offset(event_time=None)

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

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

readstream()

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

s_condition(trrf, shift)

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)

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

search_inc(bazi)

set_comp()

snr(length=50)

trim(time_before, time_after, isreturn=False)

offset = sac.b - real o

write(path, evt_datetime)

exception seispy.eq.NotEnoughComponent(matchkey)

Bases: Exception

exception seispy.eq.TooMoreComponents(matchkey)

Bases: Exception

seispy.eq.rotateZNE(st)

seispy.geo module

seispy.geo.acosd(x)

seispy.geo.asind(x)

seispy.geo.atand(x)

seispy.geo.cosd(deg)

seispy.geo.cotd(deg)

seispy.geo.deg2km(deg)

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)

seispy.geo.km2deg(km)

seispy.geo.latlon_from(lat1, lon1, azimuth, gcarc_dist)

seispy.geo.rad2deg(rad)

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

seispy.geo.rotateSeisENtoTR(E, N, BAZ)

seispy.geo.rotateSeisZENtoLQT(Z, E, N, bazi, inc)

seispy.geo.rssq(x)

seispy.geo.sdeg2skm(sdeg)

seispy.geo.sind(deg)

seispy.geo.skm2sdeg(skm)

seispy.geo.skm2srad(skm)

seispy.geo.snr(x, y)

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

Convert spherical coordinates to geographic coordinates. :param float r: radial distance from coordinate system origin

{Units: km, Range: [0, inf)}

Parameters:

• 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)

seispy.geo.srad2skm(srad)

seispy.geo.tand(deg)

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 :param allstack: stacked HK matrix :param h: 1-D array of H :param kappa: 1-D array of kappa :param ev_num: event number :return:

seispy.hk.hk()

seispy.hk.hk_sta_test()

seispy.hk.hksta(hpara, 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.hktest()

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')

Bases: object

Methods

get_events
get_stations

get_events(starttime=None, endtime=UTCDateTime(2022, 12, 13, 14, 56, 44, 404506), **kwargs)

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

get_para
read_para

property catalogpath

property criterion

property datapath

property decon_method

get_para()

classmethod read_para(cfg_file)

property rfpath

property server

class seispy.para.StaInfo

Bases: object

Methods

get_station_from_ws([server])

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

get_stainfo
load_stainfo

get_stainfo()

get_station_from_ws(server='IRIS', **kwargs)

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

Parameters:

serverstr, optional

Sever name of FDSN web-service, by default ‘IRIS’

load_stainfo(pathname, ref_comp, suffix)

seispy.plotR module

seispy.plotR.init_figure()

seispy.plotR.plot_waves(axr, axb, axs, stadata, enf=12)

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

seispy.plotR.set_fig(axr, axb, axs, stadata, xmin=-2, xmax=80)

seispy.plotRT module

seispy.plotRT.init_figure()

seispy.plotRT.plot_waves(axr, axt, axb, axr_sum, axt_sum, stadata, enf=3)

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

Plot PRFs with R and T components

Parameters:

• rfpath (str) – Path to PRFs

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

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

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

seispy.plotRT.set_fig(axr, axt, axb, axr_sum, axt_sum, stadata, station, xmin=-2, xmax=30, comp='R')

seispy.psrayp module

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

Bases: object

Methods

taup_rayp([this_dis, this_dep, taup])

param taup: Path to taup_time

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
cal_phase
channel_correct
deconv
detrend
drop_eq_snr
filter
load_stainfo
loadpjt
match_eq
pick
rotate
save_raw_data
savepjt
saverf
search_eq
trim

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

cal_phase()

channel_correct()

property date_begin

property date_end

deconv()

detrend()

drop_eq_snr(length=None)

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

load_stainfo()

classmethod loadpjt(path)

match_eq()

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

rotate(search_inc=False)

save_raw_data()

savepjt()

saverf(gauss=None)

search_eq(local=False, catalog=None)

trim()

exception seispy.rf.SACFileNotFoundError(matchkey)

Bases: Exception

seispy.rf.datestr2regex(datestr)

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

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')

seispy.rf.pickphase(eqs, para, 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)

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()

seispy.rf2depth_makedata module

class seispy.rf2depth_makedata.Station(sta_lst)

Bases: object

seispy.rf2depth_makedata.init_mat(sta_num)

seispy.rf2depth_makedata.makedata(cpara, velmod3d=None, modfolder1d=None, log=<seispy.setuplog.setuplog object>)

seispy.rf2depth_makedata.makedata3d(cpara, velmod3d, log=<seispy.setuplog.setuplog object>, raytracing3d=True)

seispy.rf2depth_makedata.rf2depth()

seispy.rfani module

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

Bases: object

Methods

plot_polar([cmap, show, outpath])

Polar map of crustal anisotropy.

baz_stack
cut_energy_waveform
init_ani_para
joint_ani
plot_correct
plot_stack_baz
radial_energy_max
rotate_to_fast_slow
search_peak
search_peak_amp
search_peak_list
xyz2grd

baz_stack(val=10)

cut_energy_waveform(idx, nb, ne)

init_ani_para()

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

plot_correct(fvd=0, dt=0.44, enf=80, outpath=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='./')

radial_energy_max()

rotate_to_fast_slow()

search_peak(energy, opt='max')

search_peak_amp()

search_peak_list(energy, opt='max')

xyz2grd(energy)

seispy.rfani.average_delay(fd, td)

seispy.rfani.joint_stack(energy_r, energy_cc, energy_tc, weight=[0.4, 0.3, 0.3])

seispy.rfcorrect module

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

Bases: object

Attributes:

stel

Methods

harmonic([tb, te])	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. :param method: 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. :type method: str, optional.
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
read_stream(stream, rayp, baz[, prime_comp])	Create RFStation instance from obspy.Stream
resample(dt)	Resample RFs with specified dt
sort([key])	Sort RFs by keys in given event, evla, evlo, evdp, dis, bazi, rayp, mag, f0

init_property
plotr
plotrt
psrf_3D_moveoutcorrect
psrf_3D_raytracing
psrf_3D_timecorrect
read_sample
slantstack

harmonic(tb=-5, te=10)

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

Returns:

harmonic_trans: numpy.ndarray, float

Harmonic components with shape of (5, nsamp), nsamp = (te-tb)/RFStation.sampling

unmodel_trans: numpy.ndarray, float

Unmodel components with shape same as harmonic_trans.

init_property(ev_num)

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. :param method: 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)

plotrt(outformat=None, **kwargs)

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)

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)

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)

read_sample(data_path)

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

Create RFStation instance from obspy.Stream

Parameters:

streamobspy.Stream

_description_

raypnumpy.ndarray or float

Ray-parameters in s/km.

baznumpy.ndarray or float

Back-azimuth

prime_compstr, optional

Prime component of RF, by default ‘R’

resample(dt)

Resample RFs with specified dt

Parameters:

dtfloat

Target sampling interval in sec

slantstack(ref_dis=None, rayp_range=None, tau_range=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

harmonic([tb, te])	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. :param method: 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. :type method: str, optional.
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
read_stream(stream, rayp, baz[, prime_comp])	Create RFStation instance from obspy.Stream
resample(dt)	Resample RFs with specified dt
sort([key])	Sort RFs by keys in given event, evla, evlo, evdp, dis, bazi, rayp, mag, f0

init_property
plotr
plotrt
psrf_3D_moveoutcorrect
psrf_3D_raytracing
psrf_3D_timecorrect
read_sample
slantstack

seispy.rfcorrect.interp_depth_model(model, lat, lon, new_dep)

Interpolate Vp and Vs from 3D velocity with a specified depth range.

Parameters:

mod3dnp.lib.npyio.NpzFile()

3D velocity loaded from a .npz file

latfloat

Latitude of position in 3D velocity model

lonfloat

Longitude of position in 3D velocity model

new_depnp.ndarray()

1D array of depths in km

Returns:

Vpnp.ndarray()

Vp in new_dep

Vsnp.ndarray()

Vs in new_dep

seispy.rfcorrect.moveoutcorrect_ref(stadatar, raypref, YAxisRange, chan='r', 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, 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='/Users/runner/.RF.log')

Bases: object

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