seispy package#
seispy.ccpprofile module#
- class seispy.ccpprofile.CCPProfile(cfg_file=None, log=None)[source]#
Bases:
object
Methods
Initialize bins of profile
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
- save_stack_data(format='npz')[source]#
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
- seispy.ccpprofile.init_profile(lat1, lon1, lat2, lon2, val)[source]#
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)[source]#
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
- classmethod read_stack_data(stack_data_path, cfg_file=None, good_depth_path=None, ismtz=False)[source]#
- save_stack_data(fname)[source]#
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
- seispy.ccp3d.gen_center_bin(center_lat, center_lon, len_lat, len_lon, val)[source]#
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#
seispy.decon module#
- class seispy.decon.RFTrace(data=Ellipsis, header=None)[source]#
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
- seispy.decon.deconit(uin, win, dt, nt=None, tshift=10, f0=2.0, itmax=400, minderr=0.001, phase='P')[source]#
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.deconwater(uin, win, dt, tshift=10.0, wlevel=0.05, f0=2.0, normalize=False, phase='P')[source]#
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.distaz module#
- class seispy.distaz.distaz(lat1, lon1, lat2, lon2)[source]#
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
- 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
seispy.eq module#
- class seispy.eq.EQ(pathname, datestr, suffix='SAC')[source]#
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
- channel_correct(switchEN=False, reverseE=False, reverseN=False)[source]#
_summary_
- Parameters:
switchEN (bool, optional) – _description_, defaults to False
reverseE (bool, optional) – _description_, defaults to False
reverseN (bool, optional) – _description_, defaults to False
- deconvolute(shift, time_after, f0=2.0, method='iter', only_r=False, itmax=400, minderr=0.001, wlevel=0.05, target_dt=None)[source]#
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
orwater
, by defaultiter
- 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
seispy.geo module#
- seispy.geo.geo2sph(dep, lat, lon)[source]#
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.rot3D(bazi, inc)[source]#
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.sph2geo(r, theta, phi)[source]#
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.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)[source]#
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)[source]#
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')[source]#
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)[source]#
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)[source]#
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.hkpara module#
seispy.io module#
seispy.mccc module#
seispy.para module#
- class seispy.para.RFPara[source]#
Bases:
object
- Attributes:
- catalogpath
- criterion
- datapath
- decon_method
- rfpath
- server
Methods
get_para
read_para
- property catalogpath#
- property criterion#
- property datapath#
- property decon_method#
- property rfpath#
- property server#
- class seispy.para.StaInfo[source]#
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
seispy.plotR module#
seispy.plotRT module#
- seispy.plotRT.plotrt(rfsta, enf=3, out_path='./', key='bazi', outformat='g', xlim=[-2, 30], show=False)[source]#
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.psrayp module#
seispy.rf module#
- class seispy.rf.RF(cfg_file=None, log=None)[source]#
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
- property date_begin#
- property date_end#
- 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')[source]#
- 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)[source]#
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.rf2depth_makedata module#
- seispy.rf2depth_makedata.makedata(cpara, velmod3d=None, modfolder1d=None, log=<seispy.setuplog.setuplog object>)[source]#
seispy.rfani module#
- class seispy.rfani.RFAni(sacdatar, tb, te, tlen=3, val=10, rayp=0.06, model='iasp91')[source]#
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
- plot_polar(cmap=<matplotlib.colors.ListedColormap object>, show=False, outpath='./')[source]#
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.
seispy.rfcorrect module#
- class seispy.rfcorrect.RFStation(data_path, only_r=False, prime_comp='R')[source]#
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
andaverage
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)[source]#
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 -5te (
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.
- jointani(tb, te, tlen=3.0, stack_baz_val=10, rayp=0.06, velmodel='iasp91', weight=[0.4, 0.4, 0.2])[source]#
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)[source]#
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
orNone
Corrected RFs in transverse component. If
only_r
isTrue
, this variable isNone
- rf_corr:
- normalize(method='single')[source]#
Normalize amplitude of each RFs. :param method: Method of normalization with
single
andaverage
avaliable.single
for normalization with max amplitude of current RF.average
for normalization with average amplitude of current station.
- 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)[source]#
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)[source]#
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_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)[source]#
- 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)[source]#
- classmethod read_stream(stream, rayp, baz, prime_comp='R')[source]#
Create RFStation instance from
obspy.Stream
- Parameters:
- stream
obspy.Stream
_description_
- rayp
numpy.ndarray
orfloat
Ray-parameters in s/km.
- baz
numpy.ndarray
orfloat
Back-azimuth
- prime_compstr, optional
Prime component of RF, by default ‘R’
- stream
- resample(dt)[source]#
Resample RFs with specified dt
- Parameters:
- dt
float
Target sampling interval in sec
- dt
- sort(key='bazi')[source]#
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)[source]#
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
andaverage
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)[source]#
Interpolate Vp and Vs from 3D velocity with a specified depth range.
- Parameters:
- mod3d
np.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_dep
np.ndarray()
1D array of depths in km
- mod3d
- Returns:
- Vp
np.ndarray()
Vp in
new_dep
- Vs
np.ndarray()
Vs in
new_dep
- Vp
- seispy.rfcorrect.moveoutcorrect_ref(stadatar, raypref, YAxisRange, chan='r', velmod='iasp91', sphere=True, phase=1)[source]#
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)[source]#
Time-to-depth conversion with S-wave backprojection.
- Parameters:
stadatar (
RFStation()
) – Data class of RFStationYAxisRange (
numpy.ndarray
) – Depth range for conversionvelmod (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. Seeseispy.psrayp()
in detail, defaults to Nonenormalize (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)[source]#
- seispy.rfcorrect.psrf_3D_migration(pplat_s, pplon_s, pplat_p, pplon_p, raylength_s, raylength_p, Tpds, dep_range, mod3d)[source]#
- 3D time difference correction with specified ray path and 3D velocity model.
The input parameters can be generated with
psrf_1D_raytracing()
.
- Parameters:
- pplat_s
np.ndarray()
2D array of latitude of S-wave in
dep_range
, (RFStation.ev_num()
,dep_range.size
)- pplon_s
np.ndarray()
2D array of longitude of S-wave in
dep_range
, (RFStation.ev_num()
,dep_range.size
)- pplat_p
np.ndarray()
2D array of latitude of P-wave in
dep_range
, (RFStation.ev_num()
,dep_range.size
)- pplon_p
np.ndarray()
2D array of longitude of P-wave in
dep_range
, (RFStation.ev_num()
,dep_range.size
)- raylength_s
np.ndarray()
2D array of ray path length of S-wave in
dep_range
, (RFStation.ev_num()
,dep_range.size
)- raylength_p
np.ndarray()
2D array of ray path length of P-wave in
dep_range
, (RFStation.ev_num()
,dep_range.size
)- Tpds
np.ndarray()
1D array of time difference in
dep_range
(dep_range.size
)- dep_range
np.ndarray()
1D array of depths in km, (
dep_range.size
)- mod3d
np.lib.npyio.NpzFile()
3D velocity loaded from a
.npz
file
- pplat_s
- Returns:
np.ndarray()
Corrected time difference in dep_range
- seispy.rfcorrect.psrf_3D_raytracing(stadatar, YAxisRange, mod3d, srayp=None, elevation=0, sphere=True)[source]#
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
andvs
.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')[source]#
Interpolate RF amplitude with specified time difference and depth range
- Parameters:
- stadatar
RFStation()
Data class of
RFStation()
- dep_range
np.ndarray()
1D array of depths in km, (
dep_range.size
)- Tpds
np.ndarray()
1D array of time difference in
dep_range
(dep_range.size
)- normalizestr, optional
Normlization option,
'sinlge'
and'average'
are available , by default ‘single’SeeRFStation.normalize()
in detail.
- stadatar
- 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()
)
- PS_RFdepth:
- seispy.rfcorrect.xps_tps_map(dep_mod, srayp, prayp, is_raylen=False, sphere=True, phase=1)[source]#
Calculate horizontal distance and time difference at depths
- Parameters:
dep_mod (
seispy.util.DepModel()
) – 1D velocity model classsrayp (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 forPpPs
, 3 forPsPs+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
- If
seispy.setuplog module#
seispy.signal module#
- seispy.signal.smooth(x, half_len=5, window='flat')[source]#
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)[source]#
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 whitenNfft (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