seispy package¶
seispy.ccpprofile module¶
- class seispy.ccpprofile.CCPProfile(cfg_file=None, log: ~seispy.setuplog.SetupLog = <seispy.setuplog.SetupLog object>)[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
- 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, includingbin_lat
,bin_lon
,profile_dis
,depth
,amp
,ci_low
,ci_high
andcount
.bin_lat
andbin_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
andci_high
mean confidence interval with bootstrap method;count
represents stacking number of each bin.
- Parameters:
format (str) – Format for stacked data
- seispy.ccpprofile.create_center_bin_profile(stations, val=5, method='linear')[source]¶
Create bins along a profile with given stations
- Parameters:
stations (
seispy.rf2depth_makedata.Station
) – Stations along a profileval (float) – The interval between two points in km
method (str) – Method for interpolation
- Returns:
The location of bins (bin_loca), and length between each bin and the start point (profile_range)
- Return type:
(numpy.array, numpy.array)
- seispy.ccpprofile.init_profile(lat1, lon1, lat2, lon2, val)[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
Class for 3-D CCP stacking, Usually used to study mantle transition zone structure.
- Parameters:
cfg_file (str, optional) – Path to configure file. If not defined a instance of CCP3D.cpara will be initialed, defaults to None
log (seispy.sutuplog.logger , optional) – A logger instance. If not defined, seispy.sutuplog.logger will be initialed, defaults to None
Methods
Initial grid points and search stations within a distance.
load_para
(cfg_file)Load parameters from configure file.
Read RFdepth data from npz file.
read_stack_data
(stack_data_path[, cfg_file, ...])Read stacked data from local.
save_good_410_660
(fname)Save good 410 and 660 km discontinuities to local.
save_stack_data
(fname)Save stacked data and parameters to local as a npz file.
search_good_410_660
([peak_410_min, ...])Search good 410 and 660 km discontinuities from stacked data.
stack
()Search conversion points falling within a bin and stack them with bootstrap method.
get_depth_err
- load_para(cfg_file)[source]¶
Load parameters from configure file.
- Parameters:
cfg_file (str) – Path to configure file.
- classmethod read_stack_data(stack_data_path, cfg_file=None, good_depth_path=None, ismtz=False)[source]¶
Read stacked data from local.
- Parameters:
stack_data_path (str) – Path to stacked data.
cfg_file (str, optional) – Path to configure file, defaults to None
good_depth_path (str, optional) – Path to good depth file, defaults to None
ismtz (bool, optional) – Whether the good depth file is in mtz format, defaults to False
- Returns:
A instance of CCP3D.
- Return type:
- save_good_410_660(fname)[source]¶
Save good 410 and 660 km discontinuities to local.
- Parameters:
fname (str) – file name of good 410 and 660 km discontinuities
- save_stack_data(fname)[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
- search_good_410_660(peak_410_min=380, peak_410_max=440, peak_660_min=630, peak_660_max=690)[source]¶
Search good 410 and 660 km discontinuities from stacked data.
- Parameters:
peak_410_min (float, optional) – Minimum depth of 410 km discontinuity, defaults to 380
peak_410_max (float, optional) – Maximum depth of 410 km discontinuity, defaults to 440
peak_660_min (float, optional) – Minimum depth of 660 km discontinuity, defaults to 630
peak_660_max (float, optional) – Maximum depth of 660 km discontinuity, defaults to 690
- seispy.ccp3d.bin_shape(cpara)[source]¶
Compute the radius of bins in degree.
- Parameters:
cpara (CCPPara) – Parameters of CCP stacking.
- Returns:
Radius of bins in degree.
- Return type:
1-D ndarray of floats with shape (n, ), where n is the number of bins.
- seispy.ccp3d.boot_bin_stack(data_bin, n_samples=3000)[source]¶
Stack data with bootstrap method.
- Parameters:
data_bin (1-D ndarray of floats) – Data falling within a bin.
n_samples (int, optional) – Number of bootstrap samples, defaults to 3000
- Returns:
Mean, confidence interval and number of data falling within a bin.
- Return type:
tuple of floats and int
- seispy.ccp3d.gen_center_bin(center_lat, center_lon, len_lat, len_lon, val)[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
Class for receiver function trace.
- Attributes:
id
Return a SEED compatible identifier of the trace.
- meta
Methods
attach_response
(inventories)Search for and attach channel response to the trace as
obspy.core.trace.Trace
.stats.response.copy
()Returns a deepcopy of the trace.
count
()Return number of data samples of the current trace.
decimate
(factor[, no_filter, strict_length])Downsample trace data by an integer factor.
deconvolute
(utr, wtr[, method])Deconvolute receiver function from waveforms.
detrend
([type])Remove a trend from the trace.
differentiate
([method])Differentiate the trace with respect to time.
filter
(type, **options)Filter the data of the current trace.
get_id
()Return a SEED compatible identifier of the trace.
integrate
([method])Integrate the trace with respect to time.
interpolate
(sampling_rate[, method, ...])Interpolate the data using various interpolation techniques.
max
()Returns the value of the absolute maximum amplitude in the trace.
normalize
([norm])Normalize the trace to its absolute maximum.
plot
(**kwargs)Create a simple graph of the current trace.
remove_response
([inventory, output, ...])Deconvolve instrument response.
remove_sensitivity
([inventory])Remove instrument sensitivity.
resample
(sampling_rate[, window, no_filter, ...])Resample trace data using Fourier method.
simulate
([paz_remove, paz_simulate, ...])Correct for instrument response / Simulate new instrument response.
slice
([starttime, endtime, nearest_sample])Return a new Trace object with data going from start to end time.
slide
(window_length, step[, offset, ...])Generator yielding equal length sliding windows of the Trace.
spectrogram
(**kwargs)Create a spectrogram plot of the trace.
split
()Split Trace object containing gaps using a NumPy masked array into several traces.
std
()Method to get the standard deviation of amplitudes in the trace.
taper
(max_percentage[, type, max_length, side])Taper the trace.
times
([type, reftime])For convenient plotting compute a NumPy array with timing information of all samples in the Trace.
trigger
(type, **options)Run a triggering algorithm on the data of the current trace.
trim
([starttime, endtime, pad, ...])Cut current trace to given start and end time.
verify
()Verify current trace object against available meta data.
write
(filename[, format])Save current trace into a file.
- classmethod deconvolute(utr, wtr, method='iter', **kwargs)[source]¶
Deconvolute receiver function from waveforms.
- Parameters:
utr (obspy.Trace) – R or Q component for the response function
wtr (obspy.Trace) – Z or L component for the source function
method (str, optional) – Method for deconvolution, defaults to ‘iter’
kwargs (dict) – Parameters for deconvolution
- Returns:
RFTrace object
- Return type:
- seispy.decon.correl(R, W, nfft)[source]¶
Correlation in frequency domain.
- Parameters:
R (np.ndarray) – numerator
W (np.ndarray) – denominator
- Returns:
Correlation in frequency domain
- Return type:
np.ndarray
- seispy.decon.deconit(uin, win, dt, nt=None, tshift=10, f0=2.0, itmax=400, minderr=0.001, phase='P')[source]¶
Iterative deconvolution using Ligorria & Ammon method. @author: Mijian Xu @ NJU Created on Wed Sep 10 14:21:38 2014
- Parameters:
uin (np.ndarray) – R or Q component for the response function
win (np.ndarray) – Z or L component for the source function
dt (float) – sample interval in second
nt (int, optional) – number of samples, defaults to None
tshift (float, optional) – Time shift before P arrival, defaults to 10.
f0 (float, optional) – Gauss factor, defaults to 2.0
itmax (int, optional) – Max iterations, defaults to 400
minderr (float, optional) – Min change in error required for stopping iterations, defaults to 0.001
phase (str, optional) – Phase of the RF, defaults to ‘P’
- Returns:
(RFI, rms, it) RF, rms and number of iterations.
- Return type:
(np.ndarray, np.ndarray, int)
- seispy.decon.deconvolute(uin, win, dt, method='iter', **kwargs)[source]¶
Deconvolute receiver function from waveforms. :param uin: R or Q component for the response function :type uin: np.ndarray :param win: Z or L component for the source function :type win: np.ndarray :param dt: sample interval in second :type dt: float :param method: Method for deconvolution, defaults to ‘iter’:type method: str, optional :param kwargs: Parameters for deconvolution :type kwargs: dict :return: RF, rms, [iter] :rtype: np.ndarray, float, int
- seispy.decon.deconwater(uin, win, dt, tshift=10.0, wlevel=0.05, f0=2.0, normalize=False, phase='P')[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.decon.gaussFilter(dt, nft, f0)[source]¶
Gaussian filter in frequency domain.
- Parameters:
dt (float) – sample interval in second
nft (int) – number of samples
f0 (float) – Gauss factor
- Returns:
Gaussian filter in frequency domain
- Return type:
np.ndarray
- seispy.decon.gfilter(x, nfft, gauss, dt)[source]¶
Apply Gaussian filter on time series.
- Parameters:
x (np.ndarray) – input trace
nfft (int) – number of samples
gauss (np.ndarray) – Gaussian filter in frequency domain
dt (float) – sample interval in second
- Returns:
Filtered data in time domain
- Return type:
np.ndarray
- seispy.decon.phaseshift(x, nfft, dt, tshift)[source]¶
Phase shift in frequency domain.
- Parameters:
x (np.ndarray) – input trace
nfft (int) – number of samples
dt (float) – sample interval in second
tshift (float) – Time shift before P arrival
- Returns:
Phase shifted data in time domain
- Return type:
np.ndarray
seispy.distaz module¶
- class seispy.distaz.distaz(lat1, lon1, lat2, lon2)[source]¶
Bases:
object
- Subroutine to calculate the Great Circle Arc distance
between two sets of geographic coordinates
Equations take from Bullen, pages 154, 155
- Owens, September 19, 1991
Sept. 25 – fixed az and baz calculations
P. Crotwell, Setember 27, 1995 Converted to c to fix annoying problem of fortran giving wrong answers if the input doesn’t contain a decimal point.
H. P. Crotwell, September 18, 1997 Java version for direct use in java programs. * * C. Groves, May 4, 2004 * Added enough convenience constructors to choke a horse and made public double * values use accessors so we can use this class as an immutable
H.P. Crotwell, May 31, 2006 Port to python, thus adding to the great list of languages to which distaz has been ported from the origin fortran: C, Tcl, Java and now python and I vaguely remember a perl port. Long live distaz!
Mijian Xu, Jan 01, 2016 Add np.ndarray to available input.
Methods
degreesToKilometers
getAz
getBaz
getDelta
- az¶
c c Make sure 0.0 is always 0.0, not 360. c
- baz¶
c c Bullen, Sec 10.2, eqn 7 / eqn 8 c c pt. 2 is unprimed, so this is technically the az c
- 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])Correct channel name for R, E, N components
decon_p
(tshift[, tcomp])Deconvolution for P wave
decon_s
(tshift, **kwargs)Deconvolution for S wave
deconvolute
(shift, time_after[, f0, method, ...])Deconvolution
detrend
()Detrend and demean
filter
([freqmin, freqmax, order])Bandpass filter
Fix channel name for R, E, N components
from_stream
(stream)Create EQ object from obspy stream
get_arrival
(model, evdp, dis[, phase])Get arrival time, ray parameter and incident angle from TauP model
get_time_offset
([event_time])Get time offset from SAC header
judge_rf
(gauss, shift, npts[, criterion, ...])Judge whether receiver function is valid
phase_trigger
(time_before, time_after[, ...])Trigger P or S phase
Read SAC files to stream
rotate
(baz[, inc, method, search_inc, baz_shift])Rotate to radial and transverse components
saverf
(path[, evtstr, shift, evla, evlo, ...])Save receiver function to SAC file
search_baz
(bazi[, time_b, time_e, offset])Search back azimuth for P wave
search_inc
(bazi)Search incident angle for S wave
set_comp
()Set component name
snr
([length])Calculate SNR
trim
(time_before, time_after[, isreturn])offset = sac.b - real o
write
(path, evt_datetime)Write raw stream to SAC files
- channel_correct(switchEN=False, reverseE=False, reverseN=False)[source]¶
Correct channel name for R, E, N components
- Parameters:
switchEN (bool, optional) – _description_, defaults to False
reverseE (bool, optional) – _description_, defaults to False
reverseN (bool, optional) – _description_, defaults to False
- decon_p(tshift, tcomp=False, **kwargs)[source]¶
Deconvolution for P wave
- Parameters:
tshift (float) – Time shift before P arrival
tcomp (bool, optional) – Whether calculate transverse component, defaults to False
- decon_s(tshift, **kwargs)[source]¶
Deconvolution for S wave
- Parameters:
tshift (float) – Time shift before P arrival
- deconvolute(shift, time_after, f0=2.0, method='iter', only_r=False, itmax=400, minderr=0.001, wlevel=0.05, target_dt=None)[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
- filter(freqmin=0.05, freqmax=1, order=4)[source]¶
Bandpass filter
- Parameters:
freqmin (float, optional) – minimum frequency, defaults to 0.05
freqmax (float, optional) – maximum frequency, defaults to 1
order (int, optional) – filter order, defaults to 4
- classmethod from_stream(stream)[source]¶
Create EQ object from obspy stream
- Parameters:
stream (obspy.Stream) – obspy stream
- Returns:
EQ object
- Return type:
- get_arrival(model, evdp, dis, phase='P')[source]¶
Get arrival time, ray parameter and incident angle from TauP model
- Parameters:
model (TauPyModel) – TauP model
evdp (float) – focal depth
dis (float) – epicentral distance
phase (str, optional) – phase name, defaults to ‘P’
- get_time_offset(event_time=None)[source]¶
Get time offset from SAC header
- Parameters:
event_time (obspy.core.utcdatetime.UTCDateTime, optional) – event time, defaults to None
- judge_rf(gauss, shift, npts, criterion='crust', rmsgate=None)[source]¶
Judge whether receiver function is valid
- Parameters:
gauss (float) – Gaussian factor
shift (float) – time shift before P arrival
npts (int) – number of points for RF
criterion (str, optional) – criterion for judging, defaults to ‘crust’
rmsgate (float, optional) – RMS gate, defaults to None
- Returns:
whether RF is valid
- Return type:
bool
- phase_trigger(time_before, time_after, prepick=True, stl=5, ltl=10)[source]¶
Trigger P or S phase
- Parameters:
time_before (float) – time before P or S arrival
time_after (float) – time after P or S arrival
prepick (bool, optional) – whether use prepick, defaults to True
stl (int, optional) – short time length for STA/LTA, defaults to 5
ltl (int, optional) – long time length for STA/LTA, defaults to 10
- rotate(baz, inc=None, method='NE->RT', search_inc=False, baz_shift=0)[source]¶
Rotate to radial and transverse components
- Parameters:
baz (float) – back azimuth
inc (float, optional) – incident angle, defaults to None
method (str, optional) – method for rotation, defaults to ‘NE->RT’
search_inc (bool, optional) – whether search incident angle, defaults to False
baz_shift (int, optional) – shift back azimuth, defaults to 0
- saverf(path, evtstr=None, shift=0, evla=-12345.0, evlo=-12345.0, evdp=-12345.0, mag=-12345.0, gauss=0, baz=-12345.0, gcarc=-12345.0, only_r=False, **kwargs)[source]¶
Save receiver function to SAC file
- Parameters:
path (str) – path to save SAC file
evtstr (str, optional) – event string, defaults to None
shift (int, optional) – time shift before P arrival, defaults to 0
evla (float, optional) – event latitude, defaults to -12345.
evlo (float, optional) – event longitude, defaults to -12345.
evdp (float, optional) – event depth, defaults to -12345.
mag (float, optional) – event magnitude, defaults to -12345.
gauss (float, optional) – Gaussian factor, defaults to 0
baz (float, optional) – back azimuth, defaults to -12345.
gcarc (float, optional) – epicentral distance, defaults to -12345.
only_r (bool, optional) – whether only save R component, defaults to False
- search_baz(bazi, time_b=10, time_e=20, offset=90)[source]¶
Search back azimuth for P wave
- Parameters:
bazi (float) – back azimuth
time_b (int, optional) – time before P arrival, defaults to 10
time_e (int, optional) – time after P arrival, defaults to 20
offset (int, optional) – offset for searching, defaults to 90
- Returns:
back azimuth and amplitude
- Return type:
(float, np.ndarray)
- search_inc(bazi)[source]¶
Search incident angle for S wave
- Parameters:
bazi (float) – back azimuth
- Returns:
incident angle
- Return type:
float
seispy.geo module¶
- seispy.geo.acosd(x)[source]¶
Inverse cosine function with degree as output
- Parameters:
x (float) – Cosine value
- Returns:
Degree
- Return type:
float
- seispy.geo.asind(x)[source]¶
Inverse sine function with degree as output
- Parameters:
x (float) – Sine value
- Returns:
Degree
- Return type:
float
- seispy.geo.atand(x)[source]¶
Inverse tangent function with degree as output
- Parameters:
x (float) – Tangent value
- Returns:
Degree
- Return type:
float
- seispy.geo.cosd(deg)[source]¶
Cosine function with degree as input
- Parameters:
deg (float) – Degree
- Returns:
Cosine value
- Return type:
float
- seispy.geo.cotd(deg)[source]¶
Cotangent function with degree as input
- Parameters:
deg (float) – Degree
- Returns:
Cotangent value
- Return type:
float
- seispy.geo.deg2km(deg)[source]¶
Convert degree to km
- Parameters:
deg (float) – Distance in degree
- Returns:
Distance in km
- Return type:
float
- 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.geoproject(lat_p, lon_p, lat1, lon1, lat2, lon2)[source]¶
Project a point to a line
- Parameters:
lat_p (float) – Latitude of the point
lon_p (float) – Longitude of the point
lat1 (float) – Latitude of the first point of the line
lon1 (float) – Longitude of the first point of the line
lat2 (float) – Latitude of the second point of the line
lon2 (float) – Longitude of the second point of the line
- Returns:
Latitude and longitude of the projected point
- Return type:
tuple
- seispy.geo.km2deg(km)[source]¶
Convert km to degree
- Parameters:
km (float) – Distance in km
- Returns:
Distance in degree
- Return type:
float
- seispy.geo.latlon_from(lat0, lon0, azimuth, gcarc_dist, ellps='WGS84')[source]¶
Determine position with given position of initial point, azimuth and distance
Accepted numeric scalar or array: -
int
-float
-numpy.floating
-numpy.integer
-list
-tuple
-array.array
-numpy.ndarray
-xarray.DataArray
-pandas.Series
- Parameters:
lat0 (float or array) – Latitude of original point
lon0 (float or array) – Longitude of original point
azimuth (float or array) – Azimuth(s) in degree
gcarc_dist (float or array) – Distance(s) between initial and terminus point(s) in degree
ellps (
str
, optional) – Ellipsoids supported bypyproj
, defaults to “WGS84”
- Returns:
- scalar or array:
Latitude(s) of terminus point(s)
- scalar or array:
Longitude(s) of terminus point(s)
- seispy.geo.rad2deg(rad)[source]¶
Convert radian to degree
- Parameters:
rad (float) – Radian
- Returns:
Degree
- Return type:
float
- 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
- Return type:
np.ndarray
- seispy.geo.rotateSeisENtoTR(e, n, baz)[source]¶
Rotate EN to TR
- Parameters:
e (np.ndarray) – East component
n (np.ndarray) – North component
baz (float) – Back-azimuth
- Returns:
T and R components
- Return type:
tuple
- seispy.geo.rotateSeisZENtoLQT(z, e, n, bazi, inc)[source]¶
Rotate ZEN to LQT
- Parameters:
z (np.ndarray) – Vertical component
e (np.ndarray) – East component
n (np.ndarray) – North component
bazi (float) – Back-azimuth
inc (float) – Incidence angle
- Returns:
L, Q and T components
- Return type:
tuple
- seispy.geo.rssq(x)[source]¶
Root sum square :param x: Input array :type x: np.ndarray :return: Root sum square :rtype: float
- seispy.geo.sdeg2skm(sdeg)[source]¶
Convert s/degree to s/km
- Parameters:
sdeg (float) – s/degree
- Returns:
s/km
- Return type:
float
- seispy.geo.sind(deg)[source]¶
Sine function with degree as input
- Parameters:
deg (float) – Degree
- Returns:
Sine value
- Return type:
float
- seispy.geo.skm2sdeg(skm)[source]¶
Convert s/km to s/degree
- Parameters:
skm (float) – s/km
- Returns:
s/degree
- Return type:
float
- seispy.geo.skm2srad(skm)[source]¶
Convert s/km to s/rad
- Parameters:
skm (float) – s/km
- Returns:
s/rad
- Return type:
float
- seispy.geo.snr(x, y)[source]¶
Signal-to-noise ratio
- Parameters:
x (np.ndarray) – Signal
y (np.ndarray) – Noise
- Returns:
SNR
- Return type:
float
- seispy.geo.sph2geo(r, theta, phi)[source]¶
Convert spherical coordinates to geographic coordinates.
- Parameters:
r (float) – radial distance from coordinate system origin {Units: km, Range: [0, inf)}
theta (float) – polar angle {Units: radians, Range: [0, π]}
phi (float) – azimuthal angle {Units: radians, Range: [-π, π]}
- Returns:
geographic coordinate conversion (lat, lon, depth) of spherical coordinates
- Return type:
(float, float, float)
- seispy.geo.spherical2cartesian(lon, lat, dep)[source]¶
Convert spherical coordinates to cartesian coordinates
- Parameters:
lon (float) – Longitude
lat (float) – Latitude
dep (float) – Depth
- Returns:
Cartesian coordinates
- Return type:
tuple
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
- Parameters:
allstack – stacked HK matrix
h – 1-D array of H
kappa – 1-D array of kappa
ev_num – event number
- Returns:
seispy.hkpara module¶
seispy.io module¶
- class seispy.io.Query(server='IRIS', **kwargs)[source]¶
Bases:
object
Methods
get_events
([starttime, endtime])Get events from IRIS
get_stations
- get_events(starttime=None, endtime=UTCDateTime(2024, 8, 13, 7, 38, 28, 336560), **kwargs)[source]¶
Get events from IRIS
- Parameters:
starttime (
obspy.UTCDateTime
, optional) – Start time of events, defaults to Noneendtime (
obspy.UTCDateTime
, optional) – End time of events, defaults to UTCDateTime.now()
- Returns:
Events
- Return type:
obspy.Catalog
seispy.mccc module¶
seispy.para module¶
- class seispy.para.RFPara[source]¶
Bases:
object
- Attributes:
- catalogpath
- criterion
- datapath
- decon_method
- rfpath
- server
Methods
read_para
(cfg_file)Read parameters from configure file :param cfg_file: Path to configure file :type cfg_file: str
get_para
- property catalogpath¶
- property criterion¶
- property datapath¶
- property decon_method¶
- classmethod read_para(cfg_file)[source]¶
Read parameters from configure file :param cfg_file: Path to configure file :type cfg_file: str
- property rfpath¶
- property server¶
- class seispy.para.StaInfo[source]¶
Bases:
object
Methods
Get station information
get_station_from_ws
(**kwargs)Get station information from web-service with given network and station or other optional condition.
link_server
(server[, user, password])_summary_
load_stainfo
(pathname, ref_comp, suffix)Load station information from SAC file
- get_station_from_ws(**kwargs)[source]¶
Get station information from web-service with given network and station or other optional condition.
seispy.plotR module¶
- class seispy.plotR.PlotR(stadata, enf=12, xlim=[2, 80], key='bazi')[source]¶
Bases:
object
Methods
Initialize figure
plot
([out_path, outformat, show])Plot receiver function in R component :param out_path: output path :type out_path: str :param outformat: output format :type outformat: str :param show: show figure :type show: bool
Plot PRFs with R component
set_fig
()Set figure
- seispy.plotR.plotr(rfsta, out_path='./', xlim=[-2, 80], key='bazi', enf=6, outformat='g', show=False)[source]¶
Plot receiver function in R component :param rfsta: Path to PRFs :type rfsta: seispy.rfcorrect.RFStation :param enf: The enlarge factor, defaults to 6 :type enf: int, optional :param out_path: The output path, defaults to current directory :type out_path: str, optional :param key: The key to sort PRFs, available for
event
,evla
,evlo
,evdp
,dis
,bazi
,rayp
,mag
,f0
, defaults tobazi
- Parameters:
outformat (str, optional) – File format of the image file, g as ‘png’, f as ‘pdf’, defaults to ‘g’
xlim (list, optional) – xlim, defaults to [-2, 80]
show (bool) – show figure
- Returns:
PlotR object
- Return type:
seispy.plotRT module¶
- class seispy.plotRT.PlotRT(stadata, enf=3, xlim=[2, 80], key='bazi')[source]¶
Bases:
object
Methods
plot
([out_path, outformat, show])Plot PRFs with R and T components
Plot PRFs with R and T components
set_fig
()Set figure
init_figure
- plot(out_path=None, outformat='g', show=False)[source]¶
Plot PRFs with R and T components
- Parameters:
rfsta (seispy.rfcorrect.RFStation) – Path to PRFs
enf (int, optional) – The enlarge factor, defaults to 3
out_path (str, optional) – The output path, defaults to current directory
key (str, optional) – The key to sort PRFs, avialible for
event
,evla
,evlo
,evdp
,dis
,bazi
,rayp
,mag
,f0
, defaults tobazi
outformat (str, optional) – File format of the image file, g as ‘png’, f as ‘pdf’, defaults to ‘g’
- seispy.plotRT.plotrt(rfsta, out_path='./', xlim=[-2, 30], key='bazi', enf=6, outformat='g', show=False)[source]¶
Plot PRFs with R and T components
- Parameters:
rfsta (seispy.rfcorrect.RFStation) – Path to PRFs
enf (int, optional) – The enlarge factor, defaults to 6
out_path (str, optional) – The output path, defaults to current directory
key (str, optional) – The key to sort PRFs, avialible for
event
,evla
,evlo
,evdp
,dis
,bazi
,rayp
,mag
,f0
, defaults tobazi
outformat (str, optional) – File format of the image file, g as ‘png’, f as ‘pdf’, defaults to ‘g’
xlim (list, optional) – xlim, defaults to [-2, 30]
show (bool) – show figure
seispy.psrayp module¶
seispy.rf module¶
- class seispy.rf.RF(cfg_file=None, log=None)[source]¶
Bases:
object
- Attributes:
- date_begin
- date_end
Methods
baz_correct
([time_b, time_e, offset, ...])Correct back-azimuth for all data
Calculate arrivals and ray parameters for all data
Correct channel components
deconv
()Deconvolution for all data
detrend
()Detrend all data
drop_eq_snr
([length, z_only])Drop earthquakes with low SNR
filter
([freqmin, freqmax, order])Filter all data
load_stainfo
([use_date_range])Load station information from local file or remote web-service
loadpjt
(path)Load project from local disk
match_eq
()Assosiate earthquakes with local data file or online data.
pick
([prepick, stl, ltl])Pick phase arrival
rotate
([search_inc])Rotate all data to ZNE or RTZ
Save raw data to local disk
savepjt
()Save project to local disk
saverf
([gauss])Save receiver functions to local disk
search_eq
([local, catalog])Search earthquakes from local or online data server
trim
()Trim waveforms from start to end
- baz_correct(time_b=10, time_e=20, offset=90, correct_angle=None)[source]¶
Correct back-azimuth for all data
- Parameters:
time_b (int, optional) – Begin time of searching, defaults to 10
time_e (int, optional) – End time of searching, defaults to 20
offset (int, optional) – Offset of searching, defaults to 90
correct_angle (float, optional) – Correct angle, defaults to None
- property date_begin¶
- property date_end¶
- drop_eq_snr(length=None, z_only=False)[source]¶
Drop earthquakes with low SNR
- Parameters:
length (int, optional) – Length of data, defaults to None
z_only (bool, optional) – Use Z component only, defaults to False
- filter(freqmin=None, freqmax=None, order=4)[source]¶
Filter all data
- Parameters:
freqmin (float, optional) – Minimum frequency, defaults to self.para.freqmin
freqmax (float, optional) – Maximum frequency, defaults to self.para.freqmax
order (int, optional) – Order of filter, defaults to 4
- load_stainfo(use_date_range=True)[source]¶
Load station information from local file or remote web-service
- Parameters:
use_date_range (bool, optional) – Use date range to search station information, defaults to True
- classmethod loadpjt(path)[source]¶
Load project from local disk
- Parameters:
path (str) – Path to project file
- Returns:
Project object
- Return type:
- pick(prepick=True, stl=5, ltl=10)[source]¶
Pick phase arrival
- Parameters:
prepick (bool, optional) – Use STA/LTA method, defaults to True
stl (int, optional) – Short time length, defaults to 5
ltl (int, optional) – Long time length, defaults to 10
- rotate(search_inc=False)[source]¶
Rotate all data to ZNE or RTZ
- Parameters:
search_inc (bool, optional) – Search incidence angle, defaults to False
- saverf(gauss=None)[source]¶
Save receiver functions to local disk
- Parameters:
gauss (float, optional) – Gaussian width, defaults to self.para.gauss
- seispy.rf.fetch_waveform(eq_lst, para, model, logger)[source]¶
Fetch waveforms from remote data server
- Parameters:
eq_lst (pandas.DataFrame) – Earthquake list
para (seispy.para.RFPara) – RFPara object
model (obspy.taup.TauPyModel) – TauPyModel object
logger (seispy.setuplog.SetupLog) – Logger
- Returns:
Earthquake list with fetched waveforms
- Return type:
pandas.DataFrame
- seispy.rf.match_eq(eq_lst, pathname, stla, stlo, logger, ref_comp='Z', suffix='SAC', offset=None, tolerance=210, dateformat='%Y.%j.%H.%M.%S')[source]¶
Match earthquakes with local SAC files
- Parameters:
eq_lst (pandas.DataFrame) – Earthquake list
pathname (str) – Path to SAC files
stla (float) – Station latitude
stlo (float) – Station longitude
logger (seispy.setuplog.SetupLog) – Logger
ref_comp (str, optional) – Reference component, defaults to ‘Z’
suffix (str, optional) – Suffix of SAC files, defaults to ‘SAC’
offset (float, optional) – Time offset between SAC files and earthquakes, defaults to None
tolerance (int, optional) – Tolerance of time offset, defaults to 210
dateformat (str, optional) – Date format of SAC files, defaults to ‘%Y.%j.%H.%M.%S’
- Returns:
Earthquake list with matched SAC files
- Return type:
pandas.DataFrame
- seispy.rf.read_catalog(logpath: str, b_time, e_time, stla: float, stlo: float, magmin=5.5, magmax=10.0, dismin=30.0, dismax=90.0, depthmin=0, depthmax=800)[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¶
class RF2depth : process cal RF2depth class sta_part, sta_full, _RFInd
- class seispy.rf2depth_makedata.RFDepth(cpara: ~seispy.ccppara.CCPPara, log: ~logging.Logger = <Logger RF2depth (INFO)>, raytracing3d=False, velmod3d=None, modfolder1d=None)[source]¶
Bases:
object
Convert receiver function to depth axis
Methods
makedata
([psphase])Convert receiver function to depth axis
- seispy.rf2depth_makedata.rf2depth()[source]¶
CLI for Convert receiver function to depth axis There’s 4 branch provided to do RF 2 depth conversion
only -d :do moveout correction
only -r : do raytracing but no moveout correction
-d and -r : do moveout correction and raytracing
-m : use {staname}.vel file for RF2depth conversion
- class seispy.rf2depth_makedata.sta_full(station, stla, stlo, stel)¶
Bases:
tuple
Methods
count
(value, /)Return number of occurrences of value.
index
(value[, start, stop])Return first index of value.
- station¶
Alias for field number 0
- stel¶
Alias for field number 3
- stla¶
Alias for field number 1
- stlo¶
Alias for field number 2
seispy.rfani module¶
- class seispy.rfani.RFAni(sacdatar, tb, te, tlen=3, val=10, rayp=0.06, model='iasp91')[source]¶
Bases:
object
Methods
baz_stack
([val])Stack RF data with back-azimuth.
joint_ani
([weight])Joint method for crustal anisotropy estimation.
plot_correct
([fvd, dt, enf, outpath])Plot the RF data with back-azimuth and time after P corrected.
plot_polar
([cmap, show, outpath])Polar map of crustal anisotropy.
plot_stack_baz
([enf, outpath])Plot the stack of RF data with back-azimuth.
Calculate the energy of radial component.
Rotate RF data to fast and slow direction.
search_peak
(energy[, opt])Search the peak of energy.
search_peak_list
(energy[, opt])Search the peak of energy.
xyz2grd
(energy)Interpolate energy to grid.
- baz_stack(val=10)[source]¶
Stack RF data with back-azimuth.
- Parameters:
- valint, optional
Interval of back-azimuth, by default 10
- joint_ani(weight=[0.5, 0.3, 0.2])[source]¶
Joint method for crustal anisotropy estimation.
- Parameters:
- weightlist, optional
Weight for three energy matrix, by default [0.5, 0.3, 0.2]
- plot_correct(fvd=0, dt=0.44, enf=80, outpath=None)[source]¶
Plot the RF data with back-azimuth and time after P corrected.
- Parameters:
- fvdint, optional
Fast velocity direction, by default 0
- dtfloat, optional
Time delay for correction, by default 0.44
- enfint, optional
Enlarge factor for back-azimuth, by default 80
- outpathstr, optional
Output path for saving the figure, by default None
- plot_polar(cmap=<matplotlib.colors.ListedColormap object>, show=False, outpath='./')[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.
- plot_stack_baz(enf=60, outpath='./')[source]¶
Plot the stack of RF data with back-azimuth.
- Parameters:
- enfint, optional
Enlarge factor for back-azimuth, by default 60
- outpathstr, optional
Output path for saving the figure, by default ‘./’
- search_peak(energy, opt='max')[source]¶
Search the peak of energy.
- Parameters:
- energynp.ndarray
Energy matrix
- optstr, optional
Option for searching peak, by default ‘max’
seispy.rfcorrect module¶
- class seispy.rfcorrect.RFStation(data_path, only_r=False, prime_comp='R')[source]¶
Bases:
object
- Attributes:
- stel
Methods
bin_stack
([key, lim, val])Stack RFs by bins of
key
with interval ofval
harmonic
([tb, te, is_stack])Harmonic decomposition for extracting anisotropic and isotropic features from the radial and transverse RFs
jointani
(tb, te[, tlen, stack_baz_val, ...])Eastimate crustal anisotropy with a joint method.
moveoutcorrect
([ref_rayp, dep_range, ...])Moveout correction with specified reference ray-parameter and depth
normalize
([method])Normalize amplitude of each RFs.
plotr
([outformat])Plot radial RFs
plotrt
([outformat])Plot radial and transverse RFs
psrf2depth
([dep_range])Time-to-depth conversion with specified depth series.
psrf_1D_raytracing
([dep_range])1D back ray tracing to obtained Ps conversion points at discret depths
psrf_3D_moveoutcorrect
(mod3dpath, **kwargs)3D moveout correction with 3D velocity model
psrf_3D_raytracing
(mod3dpath[, dep_range, srayp])3D back ray tracing to obtained Ps conversion points at discret depths
psrf_3D_timecorrect
(mod3dpath[, dep_range, ...])3D time-to-depth conversion with 3D velocity model
read_stream
(stream, rayp, baz[, prime_comp, ...])Create RFStation instance from
obspy.Stream
resample
(dt)Resample RFs with specified dt
slantstack
([ref_dis, rayp_range, tau_range])Slant stack for receiver function
sort
([key])Sort RFs by keys in given
event
,evla
,evlo
,evdp
,dis
,bazi
,rayp
,mag
,f0
- bin_stack(key='bazi', lim=[0, 360], val=10)[source]¶
Stack RFs by bins of
key
with interval ofval
- Parameters:
key (str, optional) – Key to stack, valid in [‘bazi’, ‘rayp’ (in s/km)], defaults to ‘bazi’
val (int, optional) – Interval of bins, defaults to 10 degree for bazi
- Returns:
Stacked RFs
- Return type:
dict
with keys ofdata_prime
,datat
andcount
- harmonic(tb=-5, te=10, is_stack=True)[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 10is_stack (bool, optional) – Wether stack the result, defaults to True
- Returns:
Harmonic components and unmodel components
- Return type:
numpy.ndarray
,numpy.ndarray
- jointani(tb, te, tlen=3.0, stack_baz_val=10, rayp=0.06, velmodel='iasp91', weight=[0.4, 0.4, 0.2])[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.
- Parameters:
method (str, optional) – 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.
- plotr(outformat=None, **kwargs)[source]¶
Plot radial RFs
- Parameters:
outformat (str, optional) – Output format for plot, defaults to None
- plotrt(outformat=None, **kwargs)[source]¶
Plot radial and transverse RFs
- Parameters:
outformat (str, optional) – Output format for plot, defaults to None
- psrf2depth(dep_range=array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149]), **kwargs)[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_moveoutcorrect(mod3dpath, **kwargs)[source]¶
3D moveout correction with 3D velocity model
- Parameters:
mod3dpath (str) – Path to 3D velocity model
dep_range (numpy.ndarray, optional) – Discret conversion depth, defaults to np.arange(0, 150)
srayp (numpy.lib.npyio.NpzFile, optional) – Ray-parameter lib for Ps phases, If set up to None the rayp of direct is used, defaults to None
- Returns:
2D array of RFs in depth
- Return type:
numpy.ndarray
- psrf_3D_raytracing(mod3dpath, dep_range=array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149]), srayp=None)[source]¶
3D back ray tracing to obtained Ps conversion points at discret depths
- Parameters:
mod3dpath (str) – Path to 3D velocity model
dep_range (numpy.ndarray, optional) – Discret conversion depth, defaults to np.arange(0, 150)
srayp (numpy.lib.npyio.NpzFile, optional) – Ray-parameter lib for Ps phases, If set up to None the rayp of direct is used, defaults to None
- Return pplat_s:
Latitude of conversion points
- Return pplon_s:
Longitude of conversion points
- Return tps:
Time difference of Ps at each depth
- Return type:
list
- psrf_3D_timecorrect(mod3dpath, dep_range=array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149]), normalize='single', **kwargs)[source]¶
3D time-to-depth conversion with 3D velocity model
- Parameters:
mod3dpath (str) – Path to 3D velocity model
dep_range (numpy.ndarray, optional) – Discret conversion depth, defaults to np.arange(0, 150)
normalize (str, optional) – Normalization method, defaults to ‘single’, see RFStation.normalize for detail
- Returns:
2D array of RFs in depth
- Return type:
numpy.ndarray
- classmethod read_stream(stream, rayp, baz, prime_comp='R', stream_t=None)[source]¶
Create RFStation instance from
obspy.Stream
- Parameters:
stream (
obspy.Stream()
) – Stream of RFsrayp (float or
np.ndarray()
) – Ray-parameter of RFsbaz (float or
np.ndarray()
) – Back-azimuth of RFsprime_comp (str, optional) – Prime component in RF filename.
R
orQ
for PRF andL
orZ
for SRF, defaults to ‘R’stream_t (
obspy.Stream()
, optional) – Stream of transverse RFs, defaults to None
- Returns:
RFStation instance
- Return type:
- slantstack(ref_dis=None, rayp_range=None, tau_range=None)[source]¶
Slant stack for receiver function
- Parameters:
ref_dis (int or float, optional) – reference distance, by default None
rayp_range (numpy.ndarray, optional) – range of ray parameter, by default None
tau_range (numpy.ndarray, optional) – range of tau, by default None
- sort(key='bazi')[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
bin_stack
([key, lim, val])Stack RFs by bins of
key
with interval ofval
harmonic
([tb, te, is_stack])Harmonic decomposition for extracting anisotropic and isotropic features from the radial and transverse RFs
jointani
(tb, te[, tlen, stack_baz_val, ...])Eastimate crustal anisotropy with a joint method.
moveoutcorrect
([ref_rayp, dep_range, ...])Moveout correction with specified reference ray-parameter and depth
normalize
([method])Normalize amplitude of each RFs.
plotr
([outformat])Plot radial RFs
plotrt
([outformat])Plot radial and transverse RFs
psrf2depth
([dep_range])Time-to-depth conversion with specified depth series.
psrf_1D_raytracing
([dep_range])1D back ray tracing to obtained Ps conversion points at discret depths
psrf_3D_moveoutcorrect
(mod3dpath, **kwargs)3D moveout correction with 3D velocity model
psrf_3D_raytracing
(mod3dpath[, dep_range, srayp])3D back ray tracing to obtained Ps conversion points at discret depths
psrf_3D_timecorrect
(mod3dpath[, dep_range, ...])3D time-to-depth conversion with 3D velocity model
read_stream
(stream, rayp, baz[, prime_comp, ...])Create RFStation instance from
obspy.Stream
resample
(dt)Resample RFs with specified dt
slantstack
([ref_dis, rayp_range, tau_range])Slant stack for receiver function
sort
([key])Sort RFs by keys in given
event
,evla
,evlo
,evdp
,dis
,bazi
,rayp
,mag
,f0
- seispy.rfcorrect.moveoutcorrect_ref(stadatar, raypref, YAxisRange, chan='', velmod='iasp91', sphere=True, phase=1)[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: DepModel, 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¶
- class seispy.setuplog.SetupLog(filename='/home/runner/.RF.log')[source]¶
Bases:
object
- default_level = 20¶
- default_logs = {'Batlog': ('Bat', 20, 'stream_handler', 'file_handler'), 'CCPlog': ('CCP', 20, 'stream_handler'), 'ModCreatorlog': ('ModCreator', 20, 'stream_handler'), 'PickDepthlog': ('PickDepth', 20, 'stream_handler'), 'RF2depthlog': ('RF2depth', 20, 'stream_handler'), 'RFlog': ('RF', 20, 'stream_handler')}¶
seispy.signal module¶
- seispy.signal.smooth(x, half_len=5, window='flat')[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