seispy package

seispy.bootstrap module

exception seispy.bootstrap.InstabilityWarning

Bases: UserWarning

Issued when results may be unstable.

seispy.bootstrap.bootstrap_indexes(data, n_samples=10000)

Given data points data, where axis 0 is considered to delineate points, return an generator for sets of bootstrap indexes. This can be used as a list of bootstrap indexes (with list(bootstrap_indexes(data))) as well.

seispy.bootstrap.bootstrap_indexes_array(data, n_samples=10000)

seispy.bootstrap.bootstrap_indexes_moving_block(data, n_samples=10000, block_length=3, wrap=False)

Generate moving-block bootstrap samples.

Given data points data, where axis 0 is considered to delineate points, return a generator for sets of bootstrap indexes. This can be used as a list of bootstrap indexes (with list(bootstrap_indexes_moving_block(data))) as well.

n_samples [default 10000]: the number of subsamples to generate.

block_length [default 3]: the length of block.

wrap [default False]: if false, choose only blocks within the data, making the last block for data of length L start at L-block_length. If true, choose blocks starting anywhere, and if they extend past the end of the data, wrap around to the beginning of the data again.

seispy.bootstrap.ci(data, statfunction=<function average>, alpha=0.05, n_samples=10000, method='bca', output='lowhigh', epsilon=0.001, multi=None, _iter=True)

Given a set of data data, and a statistics function statfunction that applies to that data, computes the bootstrap confidence interval for statfunction on that data. Data points are assumed to be delineated by axis 0.

data: array_like, shape (N, …) OR tuple of array_like all with shape (N, …)

Input data. Data points are assumed to be delineated by axis 0. Beyond this, the shape doesn’t matter, so long as statfunction can be applied to the array. If a tuple of array_likes is passed, then samples from each array (along axis 0) are passed in order as separate parameters to the statfunction. The type of data (single array or tuple of arrays) can be explicitly specified by the multi parameter.

statfunction: function (data, weights=(weights, optional)) -> value

This function should accept samples of data from data. It is applied to these samples individually.

If using the ABC method, the function _must_ accept a named weights parameter which will be an array_like with weights for each sample, and must return a _weighted_ result. Otherwise this parameter is not used or required. Note that numpy’s np.average accepts this. (default=np.average)

alpha: float or iterable, optional

The percentiles to use for the confidence interval (default=0.05). If this is a float, the returned values are (alpha/2, 1-alpha/2) percentile confidence intervals. If it is an iterable, alpha is assumed to be an iterable of each desired percentile.

n_samples: float, optional

The number of bootstrap samples to use (default=10000)

method: string, optional

The method to use: one of ‘pi’, ‘bca’, or ‘abc’ (default=’bca’)

output: string, optional

The format of the output. ‘lowhigh’ gives low and high confidence interval values. ‘errorbar’ gives transposed abs(value-confidence interval value) values that are suitable for use with matplotlib’s errorbar function. (default=’lowhigh’)

epsilon: float, optional (only for ABC method)

The step size for finite difference calculations in the ABC method. Ignored for all other methods. (default=0.001)

multi: boolean, optional

If False, assume data is a single array. If True, assume data is a tuple/other iterable of arrays of the same length that should be sampled together. If None, decide based on whether the data is an actual tuple. (default=None)

confidences: tuple of floats

The confidence percentiles specified by alpha

‘pi’: Percentile Interval (Efron 13.3)

The percentile interval method simply returns the 100*alphath bootstrap sample’s values for the statistic. This is an extremely simple method of confidence interval calculation. However, it has several disadvantages compared to the bias-corrected accelerated method, which is the default.

‘bca’: Bias-Corrected Accelerated (BCa) Non-Parametric (Efron 14.3) (default)

This method is much more complex to explain. However, it gives considerably better results, and is generally recommended for normal situations. Note that in cases where the statistic is smooth, and can be expressed with weights, the ABC method will give approximated results much, much faster. Note that in a case where the statfunction results in equal output for every bootstrap sample, the BCa confidence interval is technically undefined, as the acceleration value is undefined. To match the percentile interval method and give reasonable output, the implementation of this method returns a confidence interval of zero width using the 0th bootstrap sample in this case, and warns the user.

‘abc’: Approximate Bootstrap Confidence (Efron 14.4, 22.6)

This method provides approximated bootstrap confidence intervals without actually taking bootstrap samples. This requires that the statistic be smooth, and allow for weighting of individual points with a weights= parameter (note that np.average allows this). This is _much_ faster than all other methods for situations where it can be used.

To calculate the confidence intervals for the mean of some numbers:

>> boot.ci( np.randn(100), np.average )

Given some data points in arrays x and y calculate the confidence intervals for all linear regression coefficients simultaneously:

>> boot.ci( (x,y), scipy.stats.linregress )

Efron, An Introduction to the Bootstrap. Chapman & Hall 1993

seispy.bootstrap.jackknife_indexes(data)

Given data points data, where axis 0 is considered to delineate points, return a list of arrays where each array is a set of jackknife indexes.

For a given set of data Y, the jackknife sample J[i] is defined as the data set Y with the ith data point deleted.

seispy.bootstrap.randint(low, high=None, size=None, dtype=int)

Return random integers from low (inclusive) to high (exclusive).

Return random integers from the “discrete uniform” distribution of the specified dtype in the “half-open” interval [low, high). If high is None (the default), then results are from [0, low).

Note

New code should use the integers method of a default_rng() instance instead; see random-quick-start.

lowint or array-like of ints

Lowest (signed) integers to be drawn from the distribution (unless high=None, in which case this parameter is one above the highest such integer).

highint or array-like of ints, optional

If provided, one above the largest (signed) integer to be drawn from the distribution (see above for behavior if high=None). If array-like, must contain integer values

sizeint or tuple of ints, optional

Output shape. If the given shape is, e.g., (m, n, k), then m * n * k samples are drawn. Default is None, in which case a single value is returned.

dtypedtype, optional

Desired dtype of the result. Byteorder must be native. The default value is int.

New in version 1.11.0.

outint or ndarray of ints

size-shaped array of random integers from the appropriate distribution, or a single such random int if size not provided.

random_integerssimilar to randint, only for the closed

interval [low, high], and 1 is the lowest value if high is omitted.

Generator.integers: which should be used for new code.

>>> np.random.randint(2, size=10)
array([1, 0, 0, 0, 1, 1, 0, 0, 1, 0]) # random
>>> np.random.randint(1, size=10)
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0])

Generate a 2 x 4 array of ints between 0 and 4, inclusive:

>>> np.random.randint(5, size=(2, 4))
array([[4, 0, 2, 1], # random
       [3, 2, 2, 0]])

Generate a 1 x 3 array with 3 different upper bounds

>>> np.random.randint(1, [3, 5, 10])
array([2, 2, 9]) # random

Generate a 1 by 3 array with 3 different lower bounds

>>> np.random.randint([1, 5, 7], 10)
array([9, 8, 7]) # random

Generate a 2 by 4 array using broadcasting with dtype of uint8

>>> np.random.randint([1, 3, 5, 7], [[10], [20]], dtype=np.uint8)
array([[ 8,  6,  9,  7], # random
       [ 1, 16,  9, 12]], dtype=uint8)

seispy.bootstrap.subsample_indexes(data, n_samples=1000, size=0.5)

Given data points data, where axis 0 is considered to delineate points, return a list of arrays where each array is indexes a subsample of the data of size size. If size is >= 1, then it will be taken to be an absolute size. If size < 1, it will be taken to be a fraction of the data size. If size == -1, it will be taken to mean subsamples the same size as the sample (ie, permuted samples)

seispy.ccp module

seispy.ccp.ccp_profile()

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

seispy.ccp.get_sta(rfdep, stalist, line_loca, dep_axis, log)

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

seispy.ccp.line_proj(lat1, lon1, lat2, lon2)

seispy.ccp.search_pierce(rfdep, bin_loca, profile_range, stack_range, dep_axis, log, bin_radius=None, isci=False, domperiod=5)

seispy.ccp.select_sta(rfdep, stalist, line_loca, width, dep_axis, log)

seispy.ccp.stack(rfdep, cpara, log=<seispy.setuplog.setuplog object>)

Parameters

• rfdep – RFdepth struct

• cpara – parameters for CCP stacking

• log – class for seispy.setuplog.setuplog

Returns

seispy.ccp.writedat(dat_path, stack_data, stack_range, isci=False)

seispy.ccppara module

class seispy.ccppara.CCPPara

Bases: object

property bin_radius

property shape

seispy.ccppara.ccppara(cfg_file)

seispy.decov module

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

seispy.decov.decovit(uin, win, dt, nt=None, tshift=10, f0=2.0, itmax=400, minderr=0.001, info=False)

Created on Wed Sep 10 14:21:38 2014 [RFI, rms, it]=makeRFitdecon(uin,win,dt,nt,tshift,f0,itmax,minderr)

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.decov.gaussFilter(dt, nft, f0)

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

seispy.decov.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!

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, switchEN=False, reverseE=False, reverseN=False)

Bases: object

arr_correct(write_to_sac=True)

offset = sac.b - real o

deconvolute(shift, time_after, f0, phase='P', only_r=False, itmax=400, minderr=0.001, target_dt=None)

detrend()

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

fix_channel_name()

get_arrival(model, evdp, dis)

get_raypara(model, evdp, dis)

get_time_offset(event_time)

judge_rf(shift, npts, criterion='crust')

rotate(bazi, inc=None, method='NE->RT', phase='P')

saverf(path, evtstr=None, phase='P', 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)

snr(length=50, phase='P')

trim(time_before, time_after, phase='P', isreturn=False)

offset = sac.b - real o

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

Parameters

• bazi –

• inc –

Returns

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)];

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.slantstack(seis, timeaxis, rayp_range, tau_range, ref_dis, dis)

seispy.geo.snr(x, y)

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

seispy.geo.srad2skm(srad)

seispy.geo.tand(deg)

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)

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

seispy.hk.hktest()

seispy.hk.load_cyan_map()

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

property hrange

property krange

seispy.hkpara.hkpara(cfg_file)

seispy.io module

seispy.io.lsnc()

seispy.io.nc2npz(ncdata, minlat=- 90, maxlat=90, minlon=- 180, maxlon=180, mindep=0, maxdep=6371, key='dvs')

seispy.io.wsfetch(server, starttime=None, endtime=None, minlatitude=None, maxlatitude=None, minlongitude=None, maxlongitude=None, latitude=None, longitude=None, minradius=None, maxradius=None, mindepth=None, maxdepth=None, minmagnitude=None, maxmagnitude=None, magnitudetype=None, includeallorigins=None, includeallmagnitudes=None, includearrivals=None, eventid=None, limit=None, offset=None, orderby='time-asc', catalog=None, contributor=None)

seispy.mccc module

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

seispy.para module

class seispy.para.para

Bases: object

property catalogpath

property criterion

property datapath

get_para()

property imagepath

property rfpath

seispy.pickfigure module

class seispy.pickfigure.RFFigure(rfpath, width=21, height=11, dpi=100)

Bases: matplotlib.figure.Figure

butnext()

butprevious()

enlarge()

finish()

init_figure(width=21, height=11, dpi=100)

init_variables()

onclick(event)

plot()

plotbaz()

plotwave()

read_sac(dt=0.2)

reduce()

set_figure()

Set the .Figure instance the artist belongs to.

fig : .Figure

set_page()

set_ylabels()

sort_baz_()

class seispy.pickfigure.StaData(filenames, rrf, trf, baz, goodrf)

Bases: object

seispy.pickfigure.indexpags(evt_num, maxidx=20)

seispy.pickui module

seispy.plotR module

seispy.plotR.init_figure()

seispy.plotR.main()

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

seispy.plotR.plotr(station, cfg_file, enf=6)

seispy.plotR.read_process_data(lst)

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

seispy.plotRT module

seispy.plotRT.init_figure()

seispy.plotRT.main()

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

seispy.plotRT.plotrt(rfpath, enf=3, out_path='./', outformat='g')

seispy.plotRT.read_process_data(lst, resamp_dt=0.1)

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

seispy.psrayp module

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

Bases: object

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)

seispy.rf.CfgParser(cfg_file)

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

Bases: object

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

cal_phase()

property date_begin

property date_end

deconv(itmax=None, minderr=None)

detrend()

drop_eq_snr(length=None)

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

load_stainfo()

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

rotate(method='NE->RT')

saverf()

search_eq(local=False, server=None, catalog='GCMT')

trim()

seispy.rf.datestr2regex(datestr)

seispy.rf.load_station_info(pathname, ref_comp, suffix)

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

seispy.rf.prf()

seispy.rf.read_catalog(logpath, b_time, e_time, stla, stlo, magmin=5.5, magmax=10, dismin=30, dismax=90)

seispy.rf.setpar()

class seispy.rf.stainfo

Bases: object

get_stainfo()

load_stainfo(pathname, ref_comp, suffix)

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, log=<seispy.setuplog.setuplog object>)

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

seispy.rf2depth_makedata.rf2depth()

seispy.rfcorrect module

class seispy.rfcorrect.DepModel(YAxisRange, velmod='iasp91')

Bases: object

class seispy.rfcorrect.Mod3DPerturbation(modpath, YAxisRange, velmod='iasp91')

Bases: object

interpdvp(points)

interpdvs(points)

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

Bases: object

resample(dt)

seispy.rfcorrect.from_file(mode_name)

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

seispy.rfcorrect.moveoutcorrect_ref(stadatar, raypref, YAxisRange, sampling, shift, velmod='iasp91')

Parameters

• stadatar – data class of SACStation

• raypref – referred ray parameter in rad

• YAxisRange – Depth range in nd.array type

• sampling – dt

• shift – time before P

• velmod – Path to velocity model

Returns

Newdatar, EndIndex, x_s, x_p

seispy.rfcorrect.psrf2depth(stadatar, YAxisRange, sampling, shift, velmod='iasp91', velmod_3d=None, srayp=None)

Parameters

• stadatar –

• YAxisRange –

• sampling –

• shift –

• velmod –

Returns

seispy.rfcorrect.psrf_1D_raytracing(stadatar, YAxisRange, velmod='iasp91', srayp=None)

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

seispy.rfcorrect.psrf_3D_raytracing(stadatar, YAxisRange, mod3d, srayp=None)

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

stla: float

The latitude of the station

stlo: float

The longitude of the station

stadatar: object SACStation

The data class including PRFs and more parameters

YAxisRange: array_like

The depth array with the same intervals

mod3d: ‘Mod3DPerturbation’ object

The 3D velocity model with fields of dep, lat, lon, vp and vs.

seispy.rfcorrect.time2depth(stadatar, YAxisRange, Tpds)

seispy.setuplog module

class seispy.setuplog.setuplog(filename='/home/travis/.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

seispy.updatecatalog module

seispy.updatecatalog.convertinfo(info)

seispy.updatecatalog.fetch_cata(inlog='/home/travis/virtualenv/python3.6.7/lib/python3.6/site-packages/seispy/data/EventCMT.dat', outlog='')

seispy.updatecatalog.main()

seispy.updatecatalog.ndk2dat()

seispy.updatecatalog.ndkparse(ndk_str)

Module contents