API#

seismic.trace_data#

seismic.trace_data.waveform#

Waveform download and handling of raw data

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: GNU Lesser General Public License, Version 3 (https://www.gnu.org/copyleft/lesser.html)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Thursday, 18th February 2021 02:30:02 pm Last Modified: Friday, 27th September 2024 10:28:06 am

class seismic.trace_data.waveform.Store_Client(Client: Client, path: str, read_only: bool = False)[source]#

Bases: object

Client for request and local storage of waveform data

Request client that stores downloaded data for later local retrieval. When reading stored data the client reads from the local copy of the data. Inventory data is stored in the folder inventory and attached to data that is read.

Initialize the client

Parameters:

Client (obspy request client that supports the 'get_waveform' method) – obspy request client that supports the ‘get_waveform’ method
path (str) – path of the store’s sds root directory
read_only (Bool) – If True the local archive is not extended.

download_waveforms_mdl(starttime: UTCDateTime, endtime: UTCDateTime, clients: list = None, minlat: float = None, maxlat: float = None, minlon: float = None, maxlon: float = None, network: str = None, station: str = None, location: str = None, channel: str = None, channel_priorities: List[str] = None, location_priorities: List[str] = None)[source]#

Download data using the obspy’s MassDownloader. For all args except starttime and endtime, None is interpreted as a wildcard. Wildcards are generally allowed for all arguments except starttime and endtime. Channel and location priorities are given in the form of for example [‘BH[ZNE]’, ‘HH[ZNE]’]. Priority list are not used if the channel or location is specified.

Parameters:

starttime (UTCDateTime) – Starttime of the requested data
endtime (UTCDateTime) – Endtime of the requested data
clients (list or None, optional) – List of clients to use for downloading data
minlat (float or None, optional) – Minimum latitude of the requested data
maxlat (float or None, optional) – Maximum latitude of the requested data
minlon (float or None, optional) – Minimum longitude of the requested data
maxlon (float or None, optional) – Maximum longitude of the requested data
network (str or None, optional) – Network code of the requested data
station (str or None, optional) – Station code of the requested data
location (str or None, optional) – Location code of the requested data
channel (str or None, optional) – Channel code of the requested data
channel_priorities (List[str], optional) – List of channel priorities
location_priorities (List[str], optional) – List of location priorities

… seealso::: MassDownloader()

get_waveforms(network: str, station: str, location: str, channel: str, starttime: UTCDateTime, endtime: UTCDateTime, attach_response: bool = True, _check_times: bool = True) → Stream[source]#

get_available_stations(network: str = None) → List[Tuple[str, str]][source]#

Returns a list of stations for which raw data is available. Very similar to the get_all_stations() method with the difference that this one allows to filter for particular networks.

Parameters:: network (str or None, optional) – Only return stations from this network. network==None is same as a wildcard. Defaults to None.
Returns:: List of network and station codes in the form: [[net0,stat0], [net1,stat1]].
Return type:: list

read_inventory() → Inventory[source]#

Reads an up-to-date version of the corresponding Station Inventory.

Returns:: Station Inventory
Return type:: Inventory

select_inventory_or_load_remote(network: str, station: str) → Inventory[source]#

Checks whether the response for the provided station is in self.inventory if not the response will be fetched from remote.

Parameters:

network (str) – Network code
station (str) – Station Code

Returns:

An Obspy Inventory object holding the response for the requested station.

Return type:

Inventory

compute_spectrogram(network: str, station: str, channel: str, starttime: UTCDateTime, endtime: UTCDateTime, win_len: int, freq_max: float = 25, read_increment: int = 86400, remove_response: bool = True) → Tuple[ndarray, ndarray, ndarray][source]#

Computes a time series of spectrograms for the requested station and channel.

# Enter plotting function here later .. seealso:: Use function plot() to generate

Note

MPI support since version 0.4.2. Then the result is only available on rank 0. All other ranks return None. (to save RAM)

Parameters:

network (str) – network code
station (str) – station code
channel (str) – channel code
starttime (UTCDateTime) – starttime of the first window
endtime (UTCDateTime) – endtime of the last window
win_len (int) – Length of each time window in seconds
freq_max (float, optional) – Maximum frequency to compute, defaults to 25. Also determines the Nyquist frequency.
read_increment (int, optional) – Increment to read data in - does not influence the final result but has to be >= win_len, defaults to 86400

Returns:

Frequency vector, time vector and spectrogram matrix

Return type:

Tuple[np.ndarray, np.ndarray, np.ndarray]

class seismic.trace_data.waveform.FS_Client(fs='SDS_archive')[source]#

Bases: object

Request Client for reading MSEED files from file system

Initialize the client

Refer to read_from_filesystem for documentation

get_waveforms(network: str, station: str, location: str, channel: str, starttime: UTCDateTime, endtime: UTCDateTime, trim: bool = True, debug: bool = False, **kwargs) → Stream[source]#: Read data from the local file system.

seismic.trace_data.waveform.read_from_filesystem(ID: str, starttime: datetime, endtime: datetime, fs: str = 'SDS_archive', trim: bool = True, debug: bool = False) → Stream[source]#

Read data from a filesystem

Function to read miniSEED data from a given file structure by specifying ID and time interval plus a description of the file system structure.

Parameters:

ID (string) – seedID of the channel to read (NET.STA.LOC.CHA)
starttime (datetime.datetime) – start time
endtime (datetime.datetime) – end time
fs (list) – file structure descriptor
trim (bool) – switch for trimming of the stream
debug (bool) – print debugging information

Return type:

obspy.Stream

Returns:

data stream of requested data

If the switch for trimming is False the whole stream in the files is returned.

File structure descriptor fs is a list of strings or other lists indicating the elements in the file structure. Each item of this list is translated in one level of the directory structure. The first element is a string indicating the base_directory. The following elements can be strings to indicate one of the following:

%X as defined by datetime.strftime indicating an element of t
the time. e.g. %H
%NET: network name or %net for lower case network name
%STA: station name or %sta for lower case
%CHA: channel name or %cha for lower case
%LOC: location or %loc for lower case location code
string with out %

The format strings are replaced either with an element of the starttime if they correspond to a datetime specifyer or with the respective part of the seedID. A string withouta % sign is not changed. If more than one format string is required in one directory level the need to be separated within a sublist.

A format string for the ID can be followed by a pair of braces including two strings that will be used to replace the first string with the second. This can be used if the filename contains part of the ID in a different form.

If fs is a single string it is interpreted as the base directory ‘SDSdir’ of a SeisComP Data Structure (SDS) with TYPE fixed to D <SDSdir>/Year/NET/STA/CHAN.TYPE/NET.STA.LOC.CHAN.TYPE.YEAR.DAY This usage should be equivalent to obspy.clients.filesystem.sds client.

Example:: Example for a station ‘GJE’ in network ‘HEJSL’ with channel ‘BHZ’ and location ‘00’ with the start time 2010-12-24_11:36:30 and fs = ['base_dir','%Y','%b','%NET,['%j','_','%STA'','_T_', "%CHA('BH','')", '.mseed']] will be translated in a linux filename base_dir/2010/Nov/HEJSL/317_GJE_T_Z.mseed

Note

If the data contain traces of different channels in the same file with different start and endtimes the routine will not work properly when a period spans multiple files.

seismic.trace_data.waveform.get_day_in_folder(root: str, dirlist: List[str], network: str, station: str, channel: str, type: str) → UTCDateTime[source]#

Assuming that mseed files are by day (i.e., one file per day), this function will return the earliest or the latest day available (depending upon the argument passed as type).

Parameters:

root (str) – The path to the sds root folder
dirlist (List[str]) – alphabetically sorted list of available years
network (str) – The queried network code
station (str) – The queried station code
channel (str) – The queried channel code (wildcards allowed)
type (str) – either start or end

Raises:

NotImplementedError – Unknown argument for type
FileNotFoundError – No files in db

Returns:

The earliest starttime or latest endtime in utc

Return type:

UTCDateTime

seismic.correlate#

seismic.correlate.correlate#

Module to compute correlations

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Monday, 29th March 2021 07:58:18 am Last Modified: Wednesday, 19th June 2024 03:39:46 pm

class seismic.correlate.correlate.Correlator(store_client: Store_Client, options: dict)[source]#

Bases: object

Object to manage the actual Correlation (i.e., Green’s function retrieval) for the database.

Initiates the Correlator object. When executing pxcorr(), it will actually compute the correlations and save them in an hdf5 file that can be handled using CorrelationDataBase. Data has to be preprocessed before calling this (i.e., the data already has to be given in an ASDF format). Consult Preprocessor for information on how to proceed with this.

Parameters:: options (dict or str) – Dictionary containing all options for the correlation. Can also be a path to a yaml file containing all the keys required in options.

find_interstat_dist(dis: float)[source]#

Find stations in database with interstation distance smaller than dis.

If no station inventories are available, they will be downloaded.

Parameters:: dis (float) – Find all Stations with distance less than dis [in m]

Note

only the subset of the in params.yaml defined networks and stations will be queried.

find_existing_times(tag: str, channel: str = '*') → dict[source]#

Returns the already existing starttimes in form of a dictionary (see below)

Parameters:

tag (str) – The tag that the waveforms are saved under
channel (str, optional) – Channel Combination Code (e.g., CH0-CH1), wildcards accepted. Defaults to ‘*’

Returns:

Dictionary that is structured as in the example below

Return type:

dict

Examples

>>> out_dict = my_correlator.find_existing_times('mytag', 'BHZ-BHH')
>>> print(out_dict)
{'NET0.STAT0': {
    'NET1.STAT1': {'BHZ-BHH': [%list of starttimes] ,
    'NET2.STAT2': {'BHZ-BHH':[%list of starttimes]}}}

pxcorr()[source]#: Start the correlation with the parameters that were defined when initiating the object.

seismic.correlate.correlate.st_to_np_array(st: Stream, npts: int) → Tuple[ndarray, Stream][source]#

Converts an obspy stream to a matrix with the shape (npts, st.count()). Also returns the same stream but without the data arrays in tr.data.

Parameters:

st (Stream) – Input Stream
npts (int) – Maximum number of samples per Trace

Returns:

A stream and a matrix

Return type:

np.ndarray

seismic.correlate.correlate.is_in_xcombis(id1: str, id2: str, rcombis: List[str] = None) → bool[source]#: Check if the specific combination is to be calculated according to xcombinations including the channel. xcombination are expected as Net1-Net2.Sta1-Sta2.Cha1-Cha2. (Channel information can be omitted)

seismic.correlate.correlate.calc_cross_combis(st: Stream, ex_corr: dict, method: str = 'betweenStations', rcombis: List[str] = None) → list[source]#

Calculate a list of all cross correlation combination of traces in the stream: i.e. all combination with two different stations involved.

Parameters:

st (Stream) – Stream holding the tracecs to be correlated
ex_corr (dict) – dict holding the correlations that already exist in db
method (stringf) – Determines which traces of the strem are combined.
rcombis – requested combinations, only works if method==betweenStations.

seismic.correlate.correlate.sort_comb_name_alphabetically(network1: str, station1: str, network2: str, station2: str, location1: Optional[str] = '', location2: Optional[str] = '', channel1: Optional[str] = '', channel2: Optional[str] = '') → Tuple[list, list][source]#

Returns the alphabetically sorted network and station codes from the two station.

Parameters:

network1 (str) – network code of first station
station1 (str) – station code of first station
network2 (str) – Network code of second station
station2 (str) – Station Code of second Station

Returns:

A tuple containing the list of the network codes sorted and the list of the station codes sorted.

Return type:

Tuple[ list, list]

Examples

>>> net1 = 'IU'  # Network Code of first station
>>> stat1 = 'HRV'  # Station Code of first station
>>> net2 = 'XN'
>>> stat2 = 'NEP06'
>>> print(sort_comb_name_aphabetically(
        net1, stat1, net2, stat2))
(['IU', 'XN'], ['HRV', 'NEP06'])
>>> print(sort_comb_name_aphabetically(
        net2, stat2, net1, stat1))
(['IU', 'XN'], ['HRV', 'NEP06'])
>>> # Different combination
>>> net1 = 'YP'  # Network Code of first station
>>> stat1 = 'AB3'  # Station Code of first station
>>> net2 = 'XN'
>>> stat2 = 'NEP06'
>>> print(sort_comb_name_aphabetically(
        net1, stat1, net2, stat2))
(['XN', 'YP'], ['NEP06', 'AB3'])
>>> # Different combination
>>> net1 = 'XN'  # Network Code of first station
>>> stat1 = 'NEP07'  # Station Code of first station
>>> net2 = 'XN'
>>> stat2 = 'NEP06'
>>> print(sort_comb_name_aphabetically(
        net1, stat1, net2, stat2))
(['XN', 'XN'], ['NEP06', 'NEP07'])

seismic.correlate.correlate.compute_network_station_combinations(netlist: list, statlist: list, method: str = 'betweenStations', combis: List[str] = None) → Tuple[list, list][source]#

Return the network and station codes of the correlations for the provided lists of networks and stations and the queried combination method.

Parameters:

netlist (list) – List of network codes
statlist (list) – List of Station Codes
method (str, optional) – The combination method to use. Has to be one of the following: betweenStation, betweenComponents, autoComponents, allSimpleCombinations, or allCombinations, defaults to ‘betweenStations’.
combis (List[str]) – List of desired station combinations. Given as [net0-net1.stat0-stat1]. Optional.

Raises:

ValueError – for unkown combination methods.

Returns:

A tuple containing the list of the correlation network code and the list of the correlation station code.

Return type:

Tuple[list, list]

seismic.correlate.correlate.preprocess_stream(st: Stream, store_client: Store_Client, startt: UTCDateTime, endt: UTCDateTime, taper_len: float, remove_response: bool, subdivision: dict, preProcessing: List[dict] = None, **kwargs) → Stream[source]#

Does the preprocessing on a per stream basis. Most of the parameters can be fed in by using the “yaml” dict as kwargs.

Parameters:

st (obspy.core.stream.Stream) – Input Stream to be processed
store_client (Store_Client) – Store Client for the database
inv (Inventory or None) – Station response, can be None if remove_response=False.
startt (obspy.UTCDateTime) – Starttime that the stream should be clipped / padded to.
endt (obspy.UTCDateTime) – Endtime that the stream should be clipped / padded to.
taper_len (float) – If the instrument response is removed, one might want to taper to mitigate the acausal effects of the deconvolution. This is the length of such a taper in seconds.
remove_response (bool) – Should the instrument response be removed?
subdivision (dict) – Dictionary holding information about the correlation lenghts and increments.
preProcessing (List[dict], optional) – List holding information about the different external preprocessing functions to be applied, defaults to None

Raises:

ValueError – For sampling rates higher than the stream’s native sampling rate (upsampling is not permitted).

Returns:

The preprocessed stream.

Return type:

obspy.core.stream.Stream

seismic.correlate.correlate.generate_corr_inc(st: Stream, subdivision: dict, read_len: int, **kwargs) → Iterator[Stream][source]#

Subdivides the preprocessed streams into parts of equal length using the parameters cor_inc and cor_len in subdivision. This function can be acessed by several processes in parallel

Parameters:

st (obspy.core.stream.Stream) – The preprocessed input stream
subdivision (dict) – Dictionary holding the information about the correlation length and increment.
read_len (int) – Length to be read from disk in seconds

Yield:

Equal length windows, padded with nans / masked if data is missing.

Return type:

Generator[Stream]

seismic.correlate.stream#

Module to handle and access correlations as pythonic objects

Manage objects holding correlations.

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Tuesday, 20th April 2021 04:19:35 pm Last Modified: Wednesday, 19th June 2024 04:06:21 pm

class seismic.correlate.stream.CorrBulk(A: ndarray, stats: CorrStats = None, statlist: List[CorrStats] = None)[source]#

Bases: object

An object for faster computations on several correlations. The input correlation contain data from only one Station-Channel pair.

Parameters:

A (np.ndarray) – the Correlation Matrix. Each line corresponds to one correlation start. The column distance is lag time
stats (CorrStats, optional) – CorrStats file holding the header for this object, only relevant if you are reloading this file. Otherwise, use statlist as input to create this header. defaults to None
statlist (List[CorrStats], optional) – Header of each CorrTrace used to create this object, defaults to None

plot(**kwargs)[source]#: Plot the correlation traces contained.

See also

See plot_corr_bulk() for accepted parameters.

normalize(starttime: float = None, endtime: float = None, normtype: str = 'energy')[source]#

Correct amplitude variations with time.

Measure the maximum of the absolute value of the correlation matrix in a specified lapse time window and normalize the correlation traces by this values. A coherent phase in the respective lapse time window will have constant ampitude afterwards..

Parameters:

starttime (float) – Beginning of time window in seconds with respect to the zero position.
endtime (float) – end time window in seconds with respect to the zero position.
normtype (string) – one of the following ‘energy’, ‘max’, ‘absmax’, ‘abssum’ to decide about the way to calculate the normalization.

Return type:

CorrBulk

Returns:

Same object as in self, but normalised.

..note:: This action is performed in-place. If you would like to: keep the original data use copy().

clip(thres: float, axis=1)[source]#

Clip the correlation data’s upper and lower bounds to a multiple of its standard deviation. thres determines the factor and axis the axis the std should be computed over

Parameters:

thres (float) – factor of the standard deviation to clip by
axis (int) – Axis to compute the std over and, subsequently clip over. Can be None, if you wish to compute floating point rather than a vector. Then, the array will be clipped evenly.

..note:: This action is performed in-place. If you would like to: keep the original data use copy().

copy()[source]#

Returns a copy of self

Returns:: A copy of self

correct_decay()[source]#

Correct for the amplitude decay in a correlation matrix.

Due to attenuation and geometrical spreading the amplitude of the correlations decays with increasing lapse time. This decay is corrected by dividing the correlation functions by an exponential function that models the decay.

Returns:: Self but with data corrected for amplitude decay

..note:: This action is performed in-place. If you would like to: keep the original data use copy().

correct_stretch(dv: DV)[source]#

Correct stretching of correlation matrix

In the case of a homogeneous subsurface velocity change the correlation traces are stretched or compressed. This stretching can be measured with self.stretch. The resulting DV object can be passed to this function to remove the stretching from the correlation matrix.

Parameters:: dv (DV) – Velocity Change object

..note:: This action is performed in-place. If you would like to: keep the original data use copy().

correct_shift(dt: DV)[source]#

Correct a shift of the traces.

If time shifts in the (correlation) traces occur due to clock drifts or time offsets in active measurements these can be measured with measure_shift(). If the resulting time shift is passed to this function the shift is corrected for, such that if the measurement is done again no shift will be detected.

create_corr_stream(ind: List[int] = None)[source]#

Creates a CorrStream object from the current CorrBulk object. This can be useful if you want to save your postprocessed data in a hdf5 file again.

Parameters:: ind (Iterable[ind], optional) – Indices to extract. If None, all will be used. Defaults to None.
Returns:: Correlation Stream holding the same data
Return type:: CorrStream

envelope()[source]#

Calculate the envelope of a correlation matrix.

The correlation data of the correlation matrix are replaced by their Hilbert envelopes.

Returns:: self with the envelope in data

..note:: This action is performed in-place. If you would like to: keep the original data use copy().

filter(freqs: Tuple[float, float], order: int = 3)[source]#

Filters the correlation matrix in the frequency band specified in freqs using a zero phase filter of twice the order given in order.

Parameters:

freqs (Tuple) – lower and upper limits of the pass band in Hertz
order (int) – half the order of the Butterworth filter

Returns:

self

..note:: This action is performed in-place. If you would like to: keep the original data use copy().

extract_trace(method: str = 'mean', percentile: float = 50.0) → ndarray[source]#

Extract a representative trace from a correlation matrix.

Extract a correlation trace from the that best represents the correlation matrix. Method decides about method to extract the trace. The following possibilities are available

mean averages all traces in the matrix
median takes the median of all traces in the matrix
norm_mean averages the traces normalized after normalizing for
maxima
similarity_percentile averages the percentile % of traces
that best correlate with the mean of all traces. This will exclude abnormal traces. percentile = 50 will return an average of traces with correlation (with mean trace) above the median.

Parameters:

method (string) – method to extract the trace
percentile (float) – only used for method==’similarity_percentile’

Return type:

np.ndarray

Returns:

extracted trace

..note:: the extracted trace will also be saved in self.ref_trc

extract_multi_trace(win_inc: int, method: str = 'mean', percentile: float = 50.0) → List[ndarray][source]#

Extract several representative traces from a correlation matrix. (one per time window win_inc with the length = 2*win_inc). That is, the reference will be extracted with half increment overlap in each direction.

Extract a correlation trace from the one that best represents the correlation matrix. Method decides about method to extract the trace. The following possibilities are available

mean averages all traces in the matrix
median extract the median of all traces in the matrix
norm_mean averages the traces normalized after normalizing for
maxima
similarity_percentile averages the percentile % of traces
that best correlate with the mean of all traces. This will exclude abnormal traces. percentile = 50 will return an average of traces with correlation (with mean trace) above the median.

Parameters:

win_inc (int or List[int] (number of days)) – Increment between each window that a Trace should be extracted from. Given in days. Can be a list/np.ndarray or an int (constant increment). The windows’ length will be twice the increment. If set to 0, only one trace for the whole time will be used.
method (string) – method to extract the trace
percentile (float) – only used for method==’similarity_percentile’

Return type:

np.ndarray

Returns:

extracted trace

..note:: the extracted traces will also be saved in self.ref_trc

find_clock_shift(ref_trc: ndarray = None, tw: List[ndarray] = None, shift_range: int = 10, shift_steps: int = 101, sides: str = 'both', return_sim_mat: bool = False) → DV[source]#

Compute the shift of correlations as they can occur due to a clock drift.

Parameters:

ref_trc (np.ndarray, optional) – Reference trace to use for the computation, defaults to None. Will extract a single trace if = None.
tw (List[np.ndarray], optional) – Lapse Time window to use for the computation, defaults to None
shift_range (int, optional) – Maximum shift value to test (in n samples not seconds!). Defaults to 10.
shift_steps (int, optional) – Number of shift steps, defaults to 101
sides (str, optional) – Which sides to use. Can be ‘right’, or ‘both’. Defaults to ‘both’
return_sim_mat (bool, optional) – Return the similarity matrix, defaults to False

Returns:

The shift as DV object.

Return type:

measure_shift(ref_trc: Optional[ndarray] = None, tw: Optional[List[Tuple[float, float]]] = None, shift_range: float = 10, shift_steps: int = 101, sides: str = 'both', return_sim_mat: bool = False) → DV[source]#

Time shift estimate through shifting and comparison.

This function estimates shifting of the time axis of traces as it can occur if the clocks of digitizers drift.

Time shifts are estimated comparing each trace (e.g. correlation function stored in the corr_data matrix (one for each row) with shifted versions of reference trace stored in ref_trc. The range of shifting to be tested is given in shift_range in seconds. It is used in a symmetric way from -shift_range to +``shift_range``. Shifting ist tested shift_steps times. shift_steps should be an odd number to test zero shifting. The best match (shifting amount and corresponding correlation value) is calculated in specified time windows. Multiple time windows may be specified in tw.

Parameters:

ref_trc (Optional[np.ndarray], optional) – Refernce trace for the shifting, defaults to None
tw (Optional[List[float]], optional) – Time window(s) to check the shifting in, defaults to None.
shift_range (float, optional) – Maximum shift range in seconds, defaults to 10
shift_steps (int, optional) – Number of shift steps, defaults to 101
sides (str) – Side of the traces to be used for the shifting estimate (‘both’ | ‘single’). single is used for one-sided signals from active sources or if the time window shall not be symmetric. For both the time window will be mirrowd about zero lag time, e.g. [start,end] will result in time windows [-end:-start] and [start:end] being used simultaneousy
return_sim_mat (bool, optional) – Return simmilarity matrix?, defaults to False

Returns:

A DV object holding a shift value.

Return type:

mirror()[source]#

Average the causal and acausal (i.e., right and left) parts of the correlation.

..note:: This action is performed in-place. If you would like to: keep the original data use copy().

resample(starttimes: List[UTCDateTime], endtimes: List[UTCDateTime] = [])[source]#

Function to create correlation matrices with constant sampling

When created from a CorrStream the correlation matrix contains all available correlation traces but homogeneous sampling is not guaranteed as correlation functions may be missing due to data gaps. This function restructures the correlation matrix by inserting or averaging correlation functions to provide temporally homogeneous sampling. Inserted correlation functions consist of ‘nan’ if gaps are present and averaging is done if more than one correlation function falls in a bin between start_times[i] and end_times[i]. If end_time is an empty list (default) end_times[i] is set to start_times[i] + (start_times[1] - start_times[0])

Parameters:

start_times (list of class:~obspy.core.UTCDateTime objects) – list of starting times for the bins of the new sampling
end_times (list of class:~obspy.core.UTCDateTime objects) – list of end times for the bins of the new sampling

..note:: This action is performed in-place. If you want to keep: the original data use copy()

resample_time_axis(freq: float)[source]#

Resample the lapse time axis of a correlation matrix. The correlations are automatically filtered with a highpass filter of 0.4*sampling frequency to avoid aliasing. The function decides automatically whether to decimate or resample depending on the desired frequency

Parameters:: freq (float) – new sampling frequency

..note:: This action is performed in-place. If you want to keep: the original data use copy()

smooth(wsize: int, wtype: str = 'flat', axis: int = 1) → ndarray[source]#

Smoothes the correlation matrix with a given window function of the given width along the given axis. This method is based on the convolution of a scaled window with the signal. Each row/col (i.e. depending on the selected axis) is “prepared” by introducing reflected copies of it (with the window size) in both ends so that transient parts are minimized in the beginning and end part of the resulting array.

Parameters:

wsize (int) – Window size
wtype (string) – Window type. It can be one of: [‘flat’, ‘hanning’, ‘hamming’, ‘bartlett’, ‘blackman’] defaults to ‘flat’
axis (int) – Axis along with apply the filter. O: smooth along correlation lag time axis 1: smooth along time axis

Return type:

ndarray

Returns:

Filtered matrix

..note:: This action is performed in-place. If you want to keep: the original data use copy()

stretch(ref_trc: ndarray = None, tw: List[ndarray] = None, stretch_range: float = 0.1, stretch_steps: int = 101, sides: str = 'both', return_sim_mat: bool = False, ref_tr_trim: Optional[Tuple[float, float]] = None, ref_tr_stats=None, processing: Optional[dict] = None) → DV[source]#

Compute the velocity change with the stretching method (see Sens-Schönfelder and Wegler, 2006).

Parameters:

ref_trc (np.ndarray, optional) – Reference trace(s) to use for the computatio, defaults to None. Will extract a single trace if = None.
tw (List[np.ndarray], optional) – Lapse Time window to use for the computation, defaults to None
stretch_range (float, optional) – Maximum stretch value to test, defaults to 0.1
stretch_steps (int, optional) – Number of stretch steps, defaults to 101
sides (str, optional) – Which sides to use. Can be ‘left’, ‘right’ (or ‘single’), or ‘both’. Defaults to ‘both’
return_sim_mat (bool, optional) – Return the similarity matrix, defaults to False
processing (dict) – dictionary holding processing information.

Returns:

The velocity change as DV object.

Return type:

save(path: str)[source]#

Save the object to a numpy binary format (.npz)

Parameters:: path (str) – Output path

select_corr_time(starttime: UTCDateTime, endtime: UTCDateTime, include_partially_selected: bool = True)[source]#

Selects correlations that are inside of the requested time window.

refer to CorrStream.select_corr_time for details

slice(starttime: UTCDateTime, endtime: UTCDateTime, include_partial: bool = True)[source]#

return a subset of the current CorrelationBulk with data inside the requested time window.

Parameters:

starttime (UTCDateTime) – Earliest Starttime
endtime (UTCDateTime) – Latest endtime
include_partial (bool, optional) – Include time window if it’s only partially covered, defaults to True

Returns:

Another CorrBulk object

Return type:

CorrBulk

taper(width: float)[source]#

Taper the data.

Parameters:: width (float) – width to be tapered in seconds (per side)

..note:: This action is performed in-place. If you want to keep: the original data use copy()

taper_center(width: float, slope_frac: float = 0.05)[source]#

Taper the central part of a correlation matrix.

Due to electromagnetic cross-talk, signal processing or other effects the correlaton matrices are often contaminated around the zero lag time. This function tapers (multiples by zero) the central part of width width. To avoid problems with interpolation and filtering later on this is done with cosine taper.

Parameters:

width (float) – width of the central window to be tapered in seconds (total length i.e., not per side).
slope_frac (float, optional) – fraction of width used for soothing of edges, defaults to 0.05

Returns:

self

Return type:

CorrBulk

..note:: This action is performed in-place. If you want to keep: the original data use copy()

trim(starttime: float, endtime: float)[source]#

Trim the correlation matrix to the period from starttime to endtime given in seconds from the zero position, so both can be positive and negative.

Parameters:

starttime (float) – start time in seconds with respect to the zero position
order – end time in seconds with respect to the zero position

Returns:

self

..note:: This action is performed in-place. If you want to keep: the original data use copy()

wfc(ref_trc: ndarray, time_window: ndarray, sides: str, tw_start: float, tw_len: float, freq_min: float, freq_max: float, remove_nans: bool = True) → WFC[source]#

Computes the waveform coherency (WFC) between the given reference correlation(s) and correlation matrix. If several references are given, it will loop over these

See Steinmann, et. al. (2021) for details.

Parameters:

refcorr (np.ndarray) – 1 or 2D reference Correlation Trace extracted from the correlation data. If 2D, it will be interpreted as one refcorr trace per row.
tw (np.ndarray) – Lag time window to use for the computation
sides (str) – Which sides to use. Can be both, right, left, or single.
tw_start (float) – Time window start in seconds lag time.
tw_len (float) – Length of the Lapse time window in seconds.
freq_min (float) – Highpass frequency used for the computed values. (Hz)
freq_max (float) – Lowpass frequency used for the computed values. (Hz)
remove_nans – Remove nans from CorrMatrix, defaults to True

Returns:

A dictionary with one correlation array for each reference trace. The keys are using the syntax reftr_%n.

Return type:

dict

seismic.correlate.stream.read_corr_bulk(path: str) → CorrBulk[source]#

Reads a CorrBulk object from an .npz file.

Parameters:: path (str) – Path to file
Returns:: the corresponding and converted CorrBulk object
Return type:: CorrBulk

class seismic.correlate.stream.CorrStream(traces: list = None)[source]#

Bases: Stream

Baseclass to hold correlation traces. Basically just a list of the correlation traces.

create_corr_bulk(network: str = None, station: str = None, channel: str = None, location: str = None, times: Tuple[UTCDateTime, UTCDateTime] = None, inplace=True) → CorrBulk[source]#

Creates a CorrelationBulk object, which offers additional options for faster postprocessing.

Parameters:

network (str, optional) – Select only this network, defaults to None
station (str, optional) – Take data from this station, defaults to None
channel (str, optional) – Take data from this channel. If None all channels will be retained, defaults to None
location (str, optional) – Take data from only this location. Else various locations can be processed together, defaults to None
times (Tuple[UTCDateTime, UTCDateTime], optional) – Only take data from this time window, defaults to None
inplace (bool, optional) – The original data will be deleted to save memory, defaults to True.

Returns:

The CorrelationBulk object

Return type:

CorrBulk

Note

This function will check whether the metadata of the input stream is identical, so that correlations from different stations, components, or differently processed data cannot be mixed.

plot(sort_by: str = 'corr_start', timelimits: Tuple[float, float] = None, ylimits: Tuple[float, float] = None, scalingfactor: float = None, ax: Axes = None, linewidth: float = 0.25, outputfile: str = None, title: str = None, type: str = 'heatmap', cmap: str = 'inferno', vmin: float = None, vmax: float = None, **kwargs)[source]#

Creates a section plot of all correlations in this stream. kwargs will be passed to plot_cst()

Parameters:

sort_by (str, optional) – Which parameter to plot against. Can be either corr_start or distance, defaults to ‘corr_start’.
timelimits (Tuple[float, float], optional) – xlimits (lag time) in seconds, defaults to None
ylimits (Tuple[float, float], optional) – limits for Y-axis (either a datetime.datetime or float in km (if plotted against distance)), defaults to None.
scalingfactor (float, optional) – Which factor to scale the Correlations with. Play around with this if you want to make amplitudes bigger or smaller, defaults to None (automatically chosen).
ax (plt.Axes, optional) – Plot in existing axes? Defaults to None
linewidth (float, optional) – Width of the lines to plot, defaults to 0.25
outputfile (str, optional) – Save the plot? defaults to None
title (str, optional) – Title of the plot, defaults to None
type – Type of plot. Either ‘heatmap’ for a heat plot or ‘section’ for a wiggle type plot. Defaults to heatmap.
cmap – Decides about colormap if type == ‘heatmap’.

Defaults to ‘inferno’. :type cmap: str, optional

Note

If you would like to plot a subset of this stream, use select().

pop_at_utcs(utcs: ndarray[Any, dtype[UTCDateTime]])[source]#

Remove Correlations that contain any time given in utcs.

Parameters:: utcs (npt.ArrayLike[UTCDateTime]) – Array Containing UTC times that should be filtered. If UTC is between any of corr_start and corr_end, the corresponding CorrTrace will be removed from the stream.
Returns:: filtered CorrStream
Return type:: CorrStream

remove_duplicates()[source]#: Removes identical CorrTraces.

Note

This action is performed in-place. If you want to keep the original data use copy()

select_corr_time(starttime: UTCDateTime, endtime: UTCDateTime, include_partially_selected: bool = True)[source]#

Selects correlations that are inside of the requested time window.

Parameters:

starttime (UTCDateTime) – Requested start
endtime (UTCDateTime) – Requested end
include_partially_selected (bool, optional) –
If set to True, also the half selected time window before the requested time will be attached Given the following stream containing 6 correlations, “|” are the correlation starts and ends, “A” is the requested starttime and “B” the corresponding endtime:
```
|         |A        |         |       B |         |
1         2         3         4         5         6
```
include_partially_selected=True will select samples 2-4, include_partially_selected=False will select samples 3-4 only. Defaults to True

Returns:

Correlation Stream holding all selected traces

Return type:

CorrStream

select_time(start: Tuple[int, int, float], end: Tuple[int, int, float], exclude: bool = False)[source]#

Selects correlations that are inside of the requested time window.

Parameters:

start (Tuple[int, int, float]) – Start time given as Tuple in the form (h, m, s)
end (Tuple[int, int, float]) – End time (refers to corr_end) given as Tuple in the form (h, m, s).
exclude (bool, optional) – Exclude the times inside of the (start, end) interval, defaults to False

Returns:

A CorrStream holding the selected CorrTraces

Return type:

CorrStream

slide(window_length: float, step: float, include_partially_selected: bool = True, starttime: UTCDateTime = None, endtime: UTCDateTime = None)[source]#

Generator yielding correlations that are inside of each requested time window and inside of this stream.

Please keep in mind that it only returns a new view of the original data. Any modifications are applied to the original data as well. If you don’t want this you have to create a copy of the yielded windows. Also be aware that if you modify the original data and you have overlapping windows, all following windows are affected as well.

Not all yielded windows must have the same number of traces. The algorithm will determine the maximal temporal extents by analysing all Traces and then creates windows based on these times.

Parameters:

window_length (float) – The length of the requested time window in seconds. Note that the window length has to correspond at least to the length of the longest correlation window (i.e., the length of the correlated waveforms). This is because the correlations cannot be sliced.
step (float) – The step between the start times of two successive windows in seconds. Has to be greater than 0
include_partially_selected (bool, optional) –
If set to True, also the half selected time window before the requested time will be attached Given the following stream containing 6 correlations, “|” are the correlation starts and ends, “A” is the requested starttime and “B” the corresponding endtime:
```
|         |A        |         |       B |         |
1         2         3         4         5         6
```
include_partially_selected=True will select samples 2-4, include_partially_selected=False will select samples 3-4 only. Defaults to True.
starttime (UTCDateTime) – Start the sequence at this time instead of the earliest available starttime.
endtime (UTCDateTime) – Start the sequence at this time instead of the latest available endtime.

stack(weight: str = 'by_length', starttime: UTCDateTime = None, endtime: UTCDateTime = None, stack_len: int | str = 0, regard_location=True, norm_traces: bool = False)[source]#

Average the data of all traces in the given time windows. Will only stack data from the same network/channel/station combination. Location codes will only optionally be regarded.

Parameters:

starttime (UTCDateTime, optional) – starttime of the stacking time windows. If None, the earliest available is chosen, defaults to None.
endtime (UTCDateTime, optional) – endtime of the stacking time windows. If None, the latest available is chosen, defaults to None
stack_len (intorstr, optional) – Length of one stack. Is either a value in seconds, the special option “daily” (creates 24h stacks that always start at midnight), or 0 for a single stack over the whole time period, defaults to 0.
regard_location (bool, optional) – Don’t stack correlations with varying location code combinations, defaults to True.
norm_traces (bool, optional) – norm each trace by its maximum before stacking. Defaults to false

Returns:

A stream holding the stacks.

Return type:

:class`~seismic.correlate.stream.CorrStream`

class seismic.correlate.stream.CorrTrace(data: ndarray, header1: Stats = None, header2: Stats = None, inv: Inventory = None, start_lag: float = None, end_lag: float = None, _header: dict = None)[source]#

Bases: Trace

Baseclass to hold correlation data. Derived from the class Trace.

Initialise the correlation trace. Is done by combining the stats of the two Trace objects’ headers. If said headers do not contain Station information (i.e., coordinates), an Inventory with information about both stations should be provided as well.

Parameters:

data (np.ndarray) – The correlation data
header1 (Stats, optional) – header of the first trace, defaults to None
header2 (Stats, optional) – header of the second trace, defaults to None
inv (Inventory, optional) – Inventory object for the stations, defaults to None
start_lag (float) – The lag of the first sample of the correlation given in seconds.
end_lag (float) – The lag of the last sample of the correlation in seconds.
_header (dict, optional) – Already combined header, used when reading correlations from a file, defaults to None

plot(tlim: Tuple[float, float] = None, ax: Axes = None, outputdir: str = None, clean: bool = False) → Axes[source]#

Plots thios CorrelationTrace.

Parameters:

tlim (Tuple[float, float], optional) – Limits for the lapse axis in seconds, defaults to None
ax (plt.Axes, optional) – Plot in existing axes, defaults to None
outputdir (str, optional) – Save this plot? Defaults to None
clean (bool, optional) – Make a clean plot without labels & axes, defaults to False.

times() → ndarray[source]#

Convenience Function that returns an array holding the lag times of the correlation.

Returns:: Array with lag times
Return type:: np.ndarray

seismic.correlate.stream.alphabetical_correlation(header1: Stats, header2: Stats, start_lag: float, end_lag: float, data: ndarray, inv: Inventory) → Tuple[CorrStats, ndarray][source]#

Make sure that Correlations are always created in alphabetical order, so that we won’t have both a correlation for AB-CD and CD-AB. If the correlation was computed in the wrong order, the corr-data will be flipped along the t-axis.

Parameters:

header1 (Stats) – Header of the first trace.
header2 (Stats) – Header of the second trace
start_lag (float) – start lag in s
end_lag (float) – end lag in s
data (np.ndarray) – The computed cross-correlation for header1-header2
inv (Inventory) – The inventory holding the station coordinates. Only needed if coords aren’t provided in stats.

Returns:

the header for the CorrTrace and the data (will also be modified in place)

Return type:

Tuple[Stats, np.ndarray]

seismic.correlate.stream.combine_stats(stats1: Stats, stats2: Stats, start_lag: float, inv: Inventory = None) → CorrStats[source]#

Combine the meta-information of two ObsPy Trace.Stats objects

This function returns a ObsPy Stats object obtained combining the two associated with the input Traces. Namely stats1 and stats2.

The fields [‘network’,’station’,’location’,’channel’] are combined in a - separated fashion to create a “pseudo” SEED like id.

For all the others fields, only “common” information are retained: This means that only keywords that exist in both dictionaries will be included in the resulting one.

Parameters:

stats1 (Stats) – First Trace’s stats
stats2 (Stats) – Second Trace’s stats
start_lag (float) – The lag of the first sample of the correlation given in seconds (usually negative).
inv (Inventory, optional) – Inventory containing the station coordinates. Only needed if station coordinates are not in Trace.Stats. Defaults to None.

Return type:

Stats

Returns:

stats: combined Stats object

seismic.correlate.stream.compare_tr_id(tr0: Trace, tr1: Trace, regard_loc: bool = True) → bool[source]#

Check whether two traces are from the same channel, station, network, and, optionally, location. Useful for stacking

Parameters:

tr0 (Trace) – first trace
tr1 (Trace) – second trace
regard_loc (bool) – Regard the location code or not

Returns:

Bool whether the two are from the same (True) or not (False)

Return type:

bool

seismic.correlate.stream.stack_st_by_group(st: Stream, regard_loc: bool, weight: str, norm_traces: bool = False) → CorrStream[source]#

Stack all traces that belong to the same network, station, channel, and (optionally) location combination in the input stream.

Parameters:

st (Stream) – input Stream
regard_loc (bool) – Seperate data with different location code
norm_traces –

norm each trace by its maximum before stacking.
Defaults to false

type norm_traces:

bool, optional

Returns:

CorrStream

Return type:

CorrStream

seismic.correlate.stream.stack_st(st: CorrStream, weight: str, norm: bool = True) → CorrTrace[source]#

Returns an average of the data of all traces in the stream. Also adjusts the corr_start and corr_end parameters in the header.

Parameters:

st (CorrStream) – input Stream
weight – type of weigthing to use. Either mean or by_length
norm (bool) – Should the traces be normalised by their absolute maximum prior to stacking?

Returns:

Single trace with stacked data

Return type:

CorrTrace

seismic.correlate.stream.convert_statlist_to_bulk_stats(statlist: List[CorrStats], varying_loc: bool = True, varying_channel: bool = False) → CorrStats[source]#

Converts a list of CorrTrace stats objects to a single stats object that can be used for the creation of a CorrBulk object

Parameters:

statlist (List[Stats]) – list of Stats
varying_loc (bool) – Set true if the location codes vary. Defaults to True.
varying_loc – Set true if the channel codes vary (e.g., EHZ to BHZ), defaults to False.

Raises:

ValueError – raised if data does not fit together

Returns:

single Stats object

Return type:

Stats

seismic.correlate.stats#

Managing headers

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Monday, 5th July 2021 02:44:13 pm Last Modified: Monday, 30th January 2023 02:26:32 pm

class seismic.correlate.stats.CorrStats(header={})[source]#

Bases: AttribDict

From Obspy, but with the difference that some items can be lists and that corr_start and corr_end are introduced, some other differences only for correlations.

A container for additional header information of a ObsPy Trace object.

A Stats object may contain all header information (also known as meta data) of a CorrTrace object. Those headers may be accessed or modified either in the dictionary style or directly via a corresponding attribute. There are various default attributes which are required by every waveform import and export modules within ObsPy such as obspy.io.mseed.

Parameters:: header (dict or Stats, optional) – Dictionary containing meta information of a single Trace object. Possible keywords are summarized in the following Default Attributes section.

Basic Usage

>>> stats = CorrStats()
>>> stats.network = 'BW'
>>> print(stats['network'])
BW
>>> stats['station'] = 'MANZ'
>>> print(stats.station)
MANZ

Default Attributes

sampling_ratefloat, optional: Sampling rate in hertz (default value is 1.0).
deltafloat, optional: Sample distance in seconds (default value is 1.0).
calibfloat, optional: Calibration factor (default value is 1.0).
nptsint, optional: Number of sample points (default value is 0, which implies that no data is present).
networkstring, optional: Network code (default is an empty string).
locationstring, optional: Location code (default is an empty string).
stationstring, optional: Station code (default is an empty string).
channelstring, optional: Channel code (default is an empty string).
starttimeUTCDateTime, optional: Date and time of the first data sample used for correlation in UTC (default value is “1970-01-01T00:00:00.0Z”).
endtimeUTCDateTime, optional: Date and time of the last data sample used for correlation in UTC (default value is “1970-01-01T00:00:00.0Z”).
corr_start: UTCDateTime, optional: Date and time of the first data sample used for correlation in UTC (default value is “1970-01-01T00:00:00.0Z”).
corr_end: UTCDateTime, optional: Date and time of the last data sample used for correlation in UTC (default value is “1970-01-01T00:00:00.0Z”).
start_lag: float, optional: Lag of the first sample in seconds (usually negative)
end_lag: float, optional: Lag of the last sample in seconds.

Notes

The attributes sampling_rate and delta are linked to each other. If one of the attributes is modified the other will be recalculated.
```
>>> stats = Stats()
>>> stats.sampling_rate
1.0
>>> stats.delta = 0.005
>>> stats.sampling_rate
200.0
```

The attributes start_lag, npts, sampling_rate and delta are monitored and used to automatically calculate the end_lag.

>>> stats = Stats()
>>> stats.npts = 61
>>> stats.delta = 1.0
>>> stats.start_lag = -30
>>> stats.end_lag
30
>>> stats.delta = 0.5
>>> stats.end_lag
0

The attribute endtime, end_lag, and starttime are read only and cannot be modified. starttime and endtime are just simply aliases of corr_start and corr_end.

>>> stats = Stats()
>>> stats.endtime = UTCDateTime(2009, 1, 1, 12, 0, 0)
Traceback (most recent call last):
...
AttributeError: Attribute "endtime" in Stats object is read only!
>>> stats['endtime'] = UTCDateTime(2009, 1, 1, 12, 0, 0)
Traceback (most recent call last):
...
AttributeError: Attribute "endtime" in Stats object is read only!

The attribute npts will be automatically updated from the CorrTrace object.

>>> trace = CorrTrace()
>>> trace.stats.npts
0
>>> trace.data = np.array([1, 2, 3, 4])
>>> trace.stats.npts
4

The attribute component can be used to get or set the component, i.e. the last character of the channel attribute.

>>> stats = Stats()
>>> stats.channel = 'HHZ'
>>> stats.component  
'Z'
>>> stats.component = 'L'
>>> stats.channel  
'HHL'

readonly = ['endtime', 'end_lag', 'starttime']#

defaults = {'calib': 1.0, 'channel': '', 'corr_end': UTCDateTime(1970, 1, 1, 0, 0), 'corr_start': UTCDateTime(1970, 1, 1, 0, 0), 'delta': 1.0, 'end_lag': 0, 'endtime': UTCDateTime(1970, 1, 1, 0, 0), 'location': '', 'network': '', 'npts': 0, 'sampling_rate': 1.0, 'start_lag': 0, 'starttime': UTCDateTime(1970, 1, 1, 0, 0), 'station': ''}#

seismic.correlate.preprocessing_stream#

Preprocessing functions that are executed on Obspy Streams

Module that contains functions for preprocessing on obspy streams

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Tuesday, 20th July 2021 03:47:00 pm Last Modified: Tuesday, 26th September 2023 05:42:28 pm

seismic.correlate.preprocessing_stream.cos_taper_st(st: Stream, taper_len: float, lossless: bool = False, taper_at_masked: bool = False) → Stream[source]#

Applies a cosine taper to the input Stream.

Parameters:

tr (Stream) – Input Stream
taper_len (float) – Length of the taper per side
taper_at_masked (bool) – applies a split to each trace and merges again afterwards
lossless – Lossless tapering pads the trace’s ends with a copy of the trace’s data before tapering. Note that you will want to trim the trace later to remove this artificial ends.

type lossless: bool :return: Tapered Stream :rtype: Stream

Note

This action is performed in place. If you want to keep the original data use copy().

seismic.correlate.preprocessing_stream.cos_taper(tr: Trace, taper_len: float, taper_at_masked: bool, lossless: bool) → Trace[source]#

Applies a cosine taper to the input trace.

Parameters:

tr (Trace) – Input Trace
taper_len (float) – Length of the taper per side in seconds
taper_at_masked (bool) – applies a split to each trace and merges again afterwards. Lossless tapering is not supporting if the trace contains masked values in the middle of the trace
lossless – Lossless tapering pads the trace’s ends with a copy of the trace’s data before tapering. Note that you will want to trim the trace later to remove this artificial ends.

type lossless: bool :return: Tapered Trace :rtype: Trace

Note

This action is performed in place. If you want to keep the original data use copy().

seismic.correlate.preprocessing_stream.detrend_st(st: Stream, *args, **kwargs) → Stream[source]#

Detrends a stream while dealing with data gaps

Parameters:: st (Stream) – input Stream
Returns:: the obspy Stream detrended
Return type:: Stream

Note

This action is performed in place. If you want to keep the original data use copy().

seismic.correlate.preprocessing_stream.stream_filter(st: Stream, ftype: str, filter_option: dict) → Stream[source]#

Filter each trace of a Stream according to the given parameters

This faction apply the specified filter function to all the traces in the present in the input Stream.

Parameters:

ftype (str) – String that specifies which filter is applied (e.g. "bandpass"). See the Supported Filter section below for further details.
filter_option (dict) – Necessary arguments for the respective filter that will be passed on. (e.g. freqmin=1.0, freqmax=20.0 for "bandpass")

Pram parallel:

If the filtering will be run in parallel or not

Pram processes:

Number of processes to start (if None it will be equal to the number of cores available in the hosting machine)

Note

This operation is performed in place on the actual data arrays. The raw data is not accessible anymore afterwards. To keep your original data, use stream_copy() to create a copy of your stream object. This function can also work in parallel an all or a specified number of cores available in the hosting machine.

Supported Filter

'bandpass': Butterworth-Bandpass (uses obspy.signal.filter.bandpass()).
'bandstop': Butterworth-Bandstop (uses obspy.signal.filter.bandstop()).
'lowpass': Butterworth-Lowpass (uses obspy.signal.filter.lowpass()).
'highpass': Butterworth-Highpass (uses obspy.signal.filter.highpass()).
'lowpassCheby2': Cheby2-Lowpass (uses obspy.signal.filter.lowpassCheby2()).
'lowpassFIR' (experimental): FIR-Lowpass (uses obspy.signal.filter.lowpassFIR()).
'remezFIR' (experimental): Minimax optimal bandpass using Remez algorithm (uses obspy.signal.filter.remezFIR()).

seismic.correlate.preprocessing_stream.stream_mask_at_utc(st: Stream, starts: List[UTCDateTime], ends: List[UTCDateTime] = None, masklen: float = None, reverse: bool = False) → Stream[source]#

Mask the Data in the Stream between the times given by starts and ends or between starts and starts``+``masklen.

Parameters:

st (Stream) – Input Strem to be tapered
starts (List[UTCDateTime]) – Start-time (in UTC) that the masked values should start from
ends (List[UTCDateTime], optional) – End times (in UTC) of the masked values. Has to have the same length as starts. If None, masklen has to be defined, defaults to None.
masklen (float, optional) – Alternatively to providing ends, one can provide a constant length (in s) per mask, defaults to None.
reverse (bool, optional) – Only keep the data in the mask unmasked and mask everything else. Defaults to False.

Raises:

ValueError – If ends, starts, and masklen are incompatible with each other.

Warning

This function will not taper before and after the mask. If you should desire to do so use. cos_taper_st() and set taper_at_mask to True.

seismic.correlate.preprocessing_stream.trace_mask_at_utc(tr: Trace, starts: ndarray, ends: ndarray, reverse: bool)[source]#: See also

stream_mask_at_utc()

seismic.correlate.preprocessing_td#

Preprocessing functions that are executed in Time Domain

Module that contains functions for preprocessing in the time domain

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Tuesday, 20th July 2021 03:24:01 pm Last Modified: Thursday, 29th June 2023 11:15:56 am

seismic.correlate.preprocessing_td.clip(A: ndarray, args: dict, params: dict) → ndarray[source]#

Clip time series data at a multiple of the standard deviation

Set amplitudes exeeding a certain threshold to this threshold. The threshold for clipping is estimated as the standard deviation of each trace times a factor specified in args.

Note:

Traces should be demeaned before clipping.

Parameters:

A (numpy.ndarray) – time series data with time oriented along the first dimension (columns)
args (dictionary) – the only keyword allowed is std_factor describing the scaling of the standard deviation for the clipping threshold
params (dictionary) – not used here

Return type:

numpy.ndarray

Returns:

clipped time series data

seismic.correlate.preprocessing_td.detrend(A: ndarray, args: dict, params: dict) → ndarray[source]#

Detrend wrapper. Can call scipy.signal.detrend or a faster custom function.

Parameters:

A (np.ndarray) – 1 or 2D array
args (dict) – arguments for detrend function
params (dict) – params file

Raises:

ValueError – For unknown detrend method
ValueError – For unknown detrend type

Returns:

detrended data

Return type:

np.ndarray

seismic.correlate.preprocessing_td.detrend_scipy(A: ndarray, args: dict, params: dict) → ndarray[source]#: Remove trend from data

See also

scipy.signal.detrend() for input arguments

seismic.correlate.preprocessing_td.detrendqr(data: ndarray) → ndarray[source]#

Remove trend from data using QR decomposition. Faster than scipy.signal.detrend. Shamelessly adapted from NoisePy

See also

Adapted from NoisePy

Parameters:: data (np.ndarray) – 1 or 2D array
Raises:: ValueError – For data with more than 2 dimensions
Returns:: detrended data
Return type:: np.ndarray

seismic.correlate.preprocessing_td.demean(data: ndarray) → ndarray[source]#

Demean data

Parameters:: data (np.ndarray) – 1 or 2D array
Returns:: demeaned data
Return type:: np.ndarray

seismic.correlate.preprocessing_td.mute(A: ndarray, args: dict, params: dict) → ndarray[source]#

Mute parts of data that exceed a threshold

To completely surpress the effect of data with high amplitudes e.g. after aftershocks these parts of the data are muted (set to zero). The respective parts of the signal are identified as those where the envelope in a given frequency exceeds a threshold given directly as absolute numer or as a multiple of the data’s standard deviation. A taper of length taper_len is applied to smooth the edges of muted segments. Setting extend_gaps to True will ensure that the taper is applied outside the segments and data inside these segments will all zero. Edges of the data will be tapered too in this case.

Parameters:

A (numpy.ndarray) – time series data with time oriented along the first dimension (columns)
args (dictionary) –
the following keywords are allowed:
- filter:
  (dictionary) description of filter to be applied before calculation of the signal envelope. If not given the envelope is calculated from raw data. The value of the keyword filter is the same as the args for the function TDfilter.
- threshold:
  (float) absolute amplitude of threshold for muting
- std_factor:
  (float) alternativly to an absolute number the threhold can be estimated as a multiple of the standard deviation if the scaling is given in as value of the keyword std_factor. If neither threshold nor std_factor are given std_factor=1 is assumed.
- extend_gaps:
  (boolean) if True date above the threshold is guaranteed to be muted, otherwise tapering will leak into these parts. This step involves an additional convolution.
- taper_len:
  (float) length of taper for muted segments in seconds
params (dictionary) – filled automatically by pxcorr

Return type:

numpy.ndarray

Returns:

clipped time series data

Example:

args={'filter':{'type':'bandpass', 'freqmin':1., 'freqmax':6.}, 'taper_len':1., 'threshold':1000, 'std_factor':1, 'extend_gaps':True}

seismic.correlate.preprocessing_td.normalizeStandardDeviation(A: ndarray, args: dict, params: dict) → ndarray[source]#

Divide the time series by their standard deviation

Divide the amplitudes of each trace by its standard deviation.

Parameters:

A (numpy.ndarray) – time series data with time oriented along the first dimension (columns)
args (dictionary) – not used here
params (dictionary) – not used here

Return type:

numpy.ndarray

Returns:

normalized time series data

seismic.correlate.preprocessing_td.signBitNormalization(A: ndarray, args: dict, params: dict) → ndarray[source]#

One bit normalization of time series data

Return the sign of the samples (-1, 0, 1).

Parameters:

A (numpy.ndarray) – time series data with time oriented along the first dimension (columns)
args (dictionary) – not used here
params (dictionary) – not used here

Return type:

numpy.ndarray

Returns:

1-bit normalized time series data

seismic.correlate.preprocessing_td.taper(A: ndarray, args: dict, params: dict) → ndarray[source]#

Taper to the time series data

Apply a simple taper to the time series data.

args has the following structure:

args = {‘type’:type of taper,taper_args}``

type may be cosine_taper with the corresponding taper_args p the percentage of the traces to taper. Possibilities of type are given by obspy.signal.

Example:

args = {'type':'cosine_taper','p':0.1}

Parameters:

A (numpy.ndarray) – time series data with time oriented along the first dimension (columns)
args (dictionary) – arguments dictionary as described above
params (dictionary) – not used here

Return type:

numpy.ndarray

Returns:

tapered time series data

seismic.correlate.preprocessing_td.TDnormalization(A: ndarray, args: dict, params: dict) → ndarray[source]#

Amplitude dependent time domain normalization

Calculate the envelope of the filtered trace, smooth it in a window of length windowLength and normalize the waveform by this trace. The two used keywords in args are filter and `windowLength that describe the filter and the length of the envelope smoothing window, respectively.

args has the following structure:

args = {‘windowLength’:length of the envelope smoothing window in [s],’filter’:{‘type’:filterType, fargs}}``

type may be bandpass with the corresponding fargs freqmin and freqmax or highpass/lowpass with the fargs freqmin/freqmax

Example:

args = {'windowLength':5,'filter':{'type':'bandpass', 'freqmin':0.5, 'freqmax':2.}}

Parameters:

A (numpy.ndarray) – time series data with time oriented along the first dimension (columns)
args (dictionary) – arguments dictionary as described above
params (dictionary) – not used here

Return type:

numpy.ndarray

Returns:

normalized time series data

seismic.correlate.preprocessing_td.TDfilter(A: ndarray, args: dict, params: dict) → ndarray[source]#

Filter time series data

Filter in time domain. Types of filters are defined by obspy.signal.

args has the following structure:

args = {‘type’:filterType, fargs}

type may be bandpass with the corresponding fargs freqmin and freqmax or highpass/lowpass with the fargs freqmin/freqmax

Example:

args = {'type':'bandpass','freqmin':0.5,'freqmax':2.}

Parameters:

A (numpy.ndarray) – time series data with time oriented along the first dimension (columns)
args (dictionary) – arguments dictionary as described above
params (dictionary) – not used here

Return type:

numpy.ndarray

Returns:

filtered time series data

seismic.correlate.preprocessing_td.zeroPadding(A: ndarray, args: dict, params: dict, axis=1) → ndarray[source]#

Append zeros to the traces

Pad traces with zeros to increase the speed of the Fourier transforms and to avoid wrap around effects. Three possibilities for the length of the padding can be set in args['type']:

nextFastLen:
Traces are padded to a length that is the next fast fft length
avoidWrapAround:
depending on length of the trace that is to be used the padded part is just long enough to avoid wrap around
avoidWrapFastLen:
Use the next fast length that avoids wrap around

Example:

args = {'type':'avoidWrapPowerTwo'}

Parameters:

A (numpy.ndarray) – time series data with time oriented along the first dimension (columns)
args (dictionary) – arguments dictionary as described above
params (dictionary) – not used here
axis (tuple, optional) – axis to pad on

Return type:

numpy.ndarray

Returns:

zero padded time series data

seismic.correlate.preprocessing_fd#

Preprocessing functions that are executed in Frequency Domain

Module containing functions for preprocessing in the frequency domain

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Tuesday, 20th July 2021 03:40:11 pm Last Modified: Wednesday, 28th June 2023 02:00:00 pm

seismic.correlate.preprocessing_fd.FDfilter(B: ndarray, args: dict, params: dict) → ndarray[source]#

Filter Fourier-transformed data

Filter Fourier tranformed data by tapering in frequency domain. The args dictionary is supposed to contain the key flimit with a value that is a four element list or tuple defines the corner frequencies (f1, f2, f3, f4) in Hz of the cosine taper which is one between f2 and f3 and tapers to zero for f1 < f < f2 and f3 < f < f4.

Parameters:

B (numpy.ndarray) – Fourier transformed time series data with frequency orientedalong the first dimension (columns)
args (dictionary) – arguments dictionary as described above
params (dictionary) – params[‘freqs’] contains an array with the freqency values of the samples in B

Return type:

numpy.ndarray

Returns:

filtered spectal data

seismic.correlate.preprocessing_fd.FDsignBitNormalization(B: ndarray, args: dict, params: dict) → ndarray[source]#

Sign bit normalization of frequency transformed data

Divides each sample by its amplitude resulting in trace with amplidues of (-1, 0, 1). As this operation is done in frequency domain it requires two Fourier transforms and is thus quite costly but alows to be performed after other steps of frequency domain procesing e.g. whitening.

Parameters:

B (numpy.ndarray) – Fourier transformed time series data with frequency orientedalong the first dimension (columns)
args (dictionary) – not used in this function
params (dictionary) – not used in this function

Return type:

numpy.ndarray

Returns:

frequency transform of the 1-bit normalized data

seismic.correlate.preprocessing_fd.spectralWhitening(B: ndarray, args: dict, params) → ndarray[source]#

Spectal whitening of Fourier-transformed date

Normalize the amplitude spectrum of the complex spectra in B. The args dictionary may contain the keyword joint_norm. If its value is True the normalization of sets of three traces are normalized jointly by the mean of their amplitude spectra. This is useful for later rotation of correlated traces in the ZNE system into the ZRT system.

Parameters:

B (numpy.ndarray) – Fourier transformed time series data with frequency orientedalong the first dimension (columns)
args (dictionary) – arguments dictionary as described above
params (dictionary) – not used here

Return type:

numpy.ndarray

Returns:

whitened spectal data

seismic.db#

seismic.db.corr_hdf5#

Save your Correlations in h5 format.

Manages the file format and class for correlations.

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Friday, 16th April 2021 03:21:30 pm Last Modified: Wednesday, 2nd October 2024 11:28:48 am

class seismic.db.corr_hdf5.DBHandler(path, mode, compression, co, force)[source]#

Bases: File

The actual file handler of the hdf5 correlation files.

Warning

Should not be accessed directly. Access :class:`~seismic.db.corr_hdf5.CorrelationDataBase` instead.

Child object of h5py.File and inherets all its attributes and functions in addition to functions that are particularly useful for noise correlations.

Create a new file object.

See the h5py user guide for a detailed explanation of the options.

name

Name of the file on disk, or file-like object. Note: for files created with the ‘core’ driver, HDF5 still requires this be non-empty.

mode

r Readonly, file must exist (default) r+ Read/write, file must exist w Create file, truncate if exists w- or x Create file, fail if exists a Read/write if exists, create otherwise

driver

Name of the driver to use. Legal values are None (default, recommended), ‘core’, ‘sec2’, ‘direct’, ‘stdio’, ‘mpio’, ‘ros3’.

libver

Library version bounds. Supported values: ‘earliest’, ‘v108’, ‘v110’, ‘v112’ and ‘latest’. The ‘v108’, ‘v110’ and ‘v112’ options can only be specified with the HDF5 1.10.2 library or later.

userblock_size

Desired size of user block. Only allowed when creating a new file (mode w, w- or x).

swmr

Open the file in SWMR read mode. Only used when mode = ‘r’.

rdcc_nbytes

Total size of the dataset chunk cache in bytes. The default size is 1024**2 (1 MiB) per dataset. Applies to all datasets unless individually changed.

rdcc_w0

The chunk preemption policy for all datasets. This must be between 0 and 1 inclusive and indicates the weighting according to which chunks which have been fully read or written are penalized when determining which chunks to flush from cache. A value of 0 means fully read or written chunks are treated no differently than other chunks (the preemption is strictly LRU) while a value of 1 means fully read or written chunks are always preempted before other chunks. If your application only reads or writes data once, this can be safely set to 1. Otherwise, this should be set lower depending on how often you re-read or re-write the same data. The default value is 0.75. Applies to all datasets unless individually changed.

rdcc_nslots

The number of chunk slots in the raw data chunk cache for this file. Increasing this value reduces the number of cache collisions, but slightly increases the memory used. Due to the hashing strategy, this value should ideally be a prime number. As a rule of thumb, this value should be at least 10 times the number of chunks that can fit in rdcc_nbytes bytes. For maximum performance, this value should be set approximately 100 times that number of chunks. The default value is 521. Applies to all datasets unless individually changed.

track_order

Track dataset/group/attribute creation order under root group if True. If None use global default h5.get_config().track_order.

fs_strategy

The file space handling strategy to be used. Only allowed when creating a new file (mode w, w- or x). Defined as: “fsm” FSM, Aggregators, VFD “page” Paged FSM, VFD “aggregate” Aggregators, VFD “none” VFD If None use HDF5 defaults.

fs_page_size

File space page size in bytes. Only used when fs_strategy=”page”. If None use the HDF5 default (4096 bytes).

fs_persist

A boolean value to indicate whether free space should be persistent or not. Only allowed when creating a new file. The default value is False.

fs_threshold

The smallest free-space section size that the free space manager will track. Only allowed when creating a new file. The default value is 1.

page_buf_size

Page buffer size in bytes. Only allowed for HDF5 files created with fs_strategy=”page”. Must be a power of two value and greater or equal than the file space page size when creating the file. It is not used by default.

min_meta_keep

Minimum percentage of metadata to keep in the page buffer before allowing pages containing metadata to be evicted. Applicable only if page_buf_size is set. Default value is zero.

min_raw_keep

Minimum percentage of raw data to keep in the page buffer before allowing pages containing raw data to be evicted. Applicable only if page_buf_size is set. Default value is zero.

locking

The file locking behavior. Defined as:

False (or “false”) – Disable file locking
True (or “true”) – Enable file locking
“best-effort” – Enable file locking but ignore some errors
None – Use HDF5 defaults

Warning

The HDF5_USE_FILE_LOCKING environment variable can override this parameter.

Only available with HDF5 >= 1.12.1 or 1.10.x >= 1.10.7.

alignment_threshold

Together with alignment_interval, this property ensures that any file object greater than or equal in size to the alignement threshold (in bytes) will be aligned on an address which is a multiple of alignment interval.

alignment_interval

This property should be used in conjunction with alignment_threshold. See the description above. For more details, see https://portal.hdfgroup.org/display/HDF5/H5P_SET_ALIGNMENT

meta_block_size

Set the current minimum size, in bytes, of new metadata block allocations. See https://portal.hdfgroup.org/display/HDF5/H5P_SET_META_BLOCK_SIZE

Additional keywords

Passed on to the selected file driver.

add_corr_options(co: dict)[source]#

add_correlation(data: CorrTrace, tag='subdivision')[source]#

Add correlation data to the hdf5 file. Can be accessed using the get_data() method.

Parameters:

data (CorrTrace or CorrStream) – Data to save. Either a CorrTrace object or a CorrStream holding one or several traces.
tag – The tag that the data should be saved under. By convention, unstacked correlations are saved with the tag ‘subdivision’, whereas stacks are saved with the tag stack_$stacklen$, where $stacklen$ is to be replaced by the length of the stack in seconds.

Raises:

TypeError – for wrong data type.

remove_data(network: str, station: str, location: str, channel: str, tag: str, corr_start: obspy.core.utcdatetime.UTCDateTime | str)[source]#

Deletes the correlation from the file.

Parameters:

network (str) – network comb string
station (str) – station comb string
location (str) – location combination code
channel (str) – channel comb string
tag (str) – tag
corr_start (UTCDateTime | str) – Correlation Start, may either be a UTCDateTime object or a string. String is allowed to contain wildcards.

Raises:

TypeError – if corr_start is not string or UTCDateTime

get_corr_options() → dict[source]#

get_data(network: str, station: str, location: str, channel: str, tag: str, corr_start: UTCDateTime = None, corr_end: UTCDateTime = None) → CorrStream[source]#

Returns a CorrStream holding all the requested data.

Note

Wildcards are allowed for all parameters.

Parameters:

network (str) – network (combination), e.g., IU-YP
station (str) – station (combination), e.g., HRV-BRK
channel (str) – channel (combination), e.g., BZ-BR
corr_start (UTCDateTime, optional) – starttime of the time windows used to computed this correlation, defaults to None
corr_end (UTCDateTime, optional) – endtime of the time windows used to computed this correlation, defaults to None

Returns:

a CorrStream holding all the requested data.

Return type:

CorrStream

get_available_starttimes(network: str, station: str, tag: str, location: str, channel: str = '*') → dict[source]#

Returns a dictionary with channel codes as keys and available correlation starttimes as values.

..note::: Wildcards are only allowed for channel.

Parameters:

network (str) – Network code (combined code, e.g., IU-YP)
station (str) – Station code (combined)
tag (str) – Tag
location (str) – Combined Location code.
channel (str, optional) – Channel code (combined), wildcards allowed, defaults to ‘*’

Returns:

A dictionary holding the availabe starttimes for each channel combination

Return type:

dict

get_available_channels(tag: str, network: str, station: str, location: str) → List[str][source]#

Returns a list of all available channels (i.e., their combination codes) in this file.

Parameters:

tag (str) – The tag that this data was save as.
network (str) – Network combination code in the form net0-net1
station (str) – Network combination code in the form stat0-stat1
location (str) – Location combination code

Returns:

A list of all channel combinations.

Return type:

List[str]

class seismic.db.corr_hdf5.CorrelationDataBase(path: str, corr_options: dict = None, mode: str = 'a', compression: str = 'gzip3', _force: bool = False)[source]#

Bases: object

Base class to handle the hdf5 files that contain noise correlations.

Access an hdf5 file holding correlations. The resulting file can be accessed using all functionalities of h5py (for example as a dict).

Parameters:

path (str) – Full path to the file
corr_options (dict or None) – The dictionary holding the parameters for the correlation. If set to None. The mode will be set to read-only = ‘r’. Defaults to None.
mode (str, optional) – Mode to access the file. Options are: ‘a’ for all, ‘w’ for write, ‘r+’ for writing in an already existing file, or ‘r’ for read-only , defaults to ‘a’.
compression (str, optional) – The compression algorithm and compression level that the arrays should be saved with. ‘gzip3’ tends to perform well, else you could choose ‘gzipx’ where x is a digit between 1 and 9 (i.e., 9 is the highest compression) or None for fastest perfomance, defaults to ‘gzip3’.
_force (bool, optional) – allow differnt correlation options, defaults to False

Warning

Access only through a context manager (see below):

>>> with CorrelationDataBase(myfile.h5) as cdb:
>>>     type(cdb)  # This is a DBHandler
<class 'seismic.db.corr_hdf5.DBHandler'>

Example:

>>> with CorrelationDataBase(
            '/path/to/db/XN-XN.NEP06-NEP06.h5') as cdb:
>>>     # find the available tags for existing db
>>>     print(list(cdb.keys()))
['co', 'recombined', 'stack_86398', 'subdivision']
>>>     # find available channels with tag subdivision
>>>     print(cdb.get_available_channels(
>>>         'subdivision', 'XN-XN', 'NEP06-NEP06'))
['HHE-HHN', 'HHE-HHZ', 'HHN-HHZ']
>>>     # Get Data from all times, specific channel and tag
>>>     st = cdb.get_data(
>>>         'XN-XN', 'NEP06-NEP06', 'HHE-HHN', 'subdivision')
>>> print(st.count())
250

seismic.db.corr_hdf5.all_traces_recursive(group: Group, stream: CorrStream, pattern: str) → CorrStream[source]#

Recursively, appends all traces in a h5py group to the input stream. In addition this will check whether the data matches a certain pattern.

Parameters:

group (class:h5py._hl.group.Group) – group to search through
stream (CorrStream) – Stream to append the traces to
pattern (str) – pattern for the path in the hdf5 file, see fnmatch for details.

Returns:

Stream with appended traces

Return type:

CorrStream

seismic.db.corr_hdf5.convert_header_to_hdf5(dataset: Dataset, header: Stats)[source]#

Convert an Stats object and adds it to the provided hdf5 dataset.

Parameters:

dataset (h5py.Dataset) – the dataset that the header should be added to
header (Stats) – The trace’s header

seismic.db.corr_hdf5.read_hdf5_header(dataset: Dataset) → Stats[source]#

Takes an hdft5 dataset as input and returns the header of the CorrTrace.

Parameters:: dataset (h5py.Dataset) – The dataset to be read from
Returns:: The trace’s header
Return type:: Stats

seismic.db.corr_hdf5.co_to_hdf5(co: dict) → dict[source]#

seismic.monitor#

seismic.monitor.monitor#

Compute seismic velocity changes

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Thursday, 3rd June 2021 04:15:57 pm Last Modified: Tuesday, 9th July 2024 04:48:54 pm

class seismic.monitor.monitor.Monitor(options: dict | str)[source]#

Bases: object

Object that handles the computation of seismic velocity changes. This will access correlations that have been computed previously with the Correlator class.

Takes the parameter file as input (either as yaml or as dict).

Parameters:: options (dict or str) – Parameters for the monitor. Usually, they come in form of a yaml file. If this parameter is a str, it will be interpreted to be the path to said yaml file.

compute_velocity_change(corr_file: str, tag: str, network: str, station: str, location: str, channel: str, ref_trcs: ndarray = None)[source]#

Computes the velocity change for a cross (or auto) correlation time-series. This process is executed “per station combination”. Meaning, several processes of this method can be started via MPI and the method compute_velocity_change_bulk(), which is the function that should rather be acessed than this one. Data will be written into the defined dv folder and plotted if this option was previously set to True.

Parameters:

corr_file (str) – File to read the correlation data from.
tag (str) – Tag of the data in the file (almost always ‘subdivision’)
network (str) – Network combination code
station (str) – Station Combination Code
location (str) – Location combination code
channel (str) – Channel combination code
ref_trcs (np.ndarray, optional) – Feed in on or several custom reference traces. If None the program will determine a reference trace from the chosen method in the config. Defaults to None

compute_velocity_change_bulk()[source]#

Compute the velocity change for all correlations using MPI. The parameters for the correlation defined in the yaml file will be used.

This function will just call compute_velocity_change() several times.

compute_components_average(method: str = 'AutoComponents')[source]#

Averages the Similarity matrix of different velocity changes in the whole dv folder. Based upon those values, new dvs and correlations are computed.

Parameters:: method (str, optional) – Which components should be averaged? Can be ‘StationWide’ if all component-combinations for one station should be averaged, ‘AutoComponents’ if only the dv results of autocorrelations are to be averaged, ‘CrossComponents’ if you wish to average dvs from the same station but only intercomponent corrrelations, or ‘CrossStations’ if cross-station correlations should be averaged (same station combination and all component combinations). Defaults to ‘AutoComponents’
Raises:: ValueError – For Unknown combination methods.

compute_waveform_coherence_bulk()[source]#

Compute the WFC for all specified (params file) correlations using MPI. The parameters for the correlation defined in the yaml file will be used.

This function will just call compute_waveform_coherence() several times. Subsequently, the average of the different component combinations will be computed.

compute_waveform_coherence(corr_file: str, tag: str, network: str, station: str, location, channel: str) → WFC[source]#

Computes the waveform coherence corresponding to one correlation (i.e., one single file).

The waveform coherence can be used as a measure of stability of a certain correlation. See Steinmann, et. al. (2021) for details.

Parameters:

corr_file (str) – File to compute the wfc from
tag (str) – Tag inside of hdf5 file to retrieve corrs from.
network (str) – Network combination code.
station (str) – Station combination code
location (str) – Location combination code
channel (str) – Channel combination code.

Raises:

ValueError – For short correlations

Returns:

An object holding the waveform coherence

Return type:

WFC

See also

To compute wfc for several correlations use: waveform_coherence _bulk().

seismic.monitor.monitor.make_time_list(start_date: str, end_date: str, date_inc: int, win_len: int) → Tuple[ndarray, ndarray][source]#

Returns numpy array of starttimes and endtimes, for which each data point of velocity change should be computed.

Parameters:

start_date (str) – Date (including time) to start monitoring
end_date (str) – Last date (including time) of monitoring
date_inc (int) – Increment in seconds between each datapoint
win_len (int) – Window Length, for which each data point is computed.

Returns:

A numpy array holding the starttimes and endtimes of the computing time windows in seconds (as UTC timestamps)

Return type:

Tuple[ np.ndarray, np.ndarray]

..note:

see `obspy's documentation
<https://docs.obspy.org/packages/autogen/obspy.core.utcdatetime.UTCDateTime.html>`
for compatible input strings.

seismic.monitor.monitor.corr_find_filter(indir: str, net: dict, **kwargs) → Tuple[List[str], List[str], List[str]][source]#

Finds the files of available correlations.
Out of those files, it will extract the ones that are requested.

Parameters:

indir (str) – Directory that the correlations are saved in.
net (dict) – Network dictionary from yaml file.

Returns:

list of network combination codes, list of station combination codes, and list of the file paths to their corresponding correlation files.

Return type:

Tuple[List[str], List[str], List[str]]

seismic.monitor.monitor.correct_dv_shift(dv0: DV, dv1: DV, method: str = 'mean', n_overlap: int = 0) → Tuple[DV, DV][source]#

Shift dv0 with respect to dv1, so that the same times will have the same value.

Note

in-place operation.

Parameters:

dv0 (DV) – DV object to be shifted
dv1 (DV) – DV object to compare to
method (str, optional) – either mean or median, defaults to ‘mean’. If mean is chosen the points will be weighted by the correlation coefficient.
n_overlap (int, optional) – Number of points to compare to beyond the already overlapping point. Makes sense if, e.g., one station replaces the other, defaults to 0.

Raises:

ValueError – If the two dvs only have non-nan value more than n_overlap apart.

Returns:

The two dv object. Only the first one is modified. Modification also happens in-place.

Return type:

Tuple[DV, DV]

seismic.monitor.monitor.correct_time_shift_several(dvs: List[DV], method: str, n_overlap: int)[source]#

Suppose you have a number of dv/v time series from the same location, but station codes changed several times. To stitch these together, this code will iteratively shift each progressive time-series to the shift of the end of its precessor.

Note

The shifting correction happens in-place.

Parameters:

dvs (List[DV]) – List of dvs to be shifted
method (str) – method to be used to measure the shift. Can be ‘mean’ or ‘median’. If ‘mean’ is chosen, the points will be weighted by their respective correlation coefficient.
n_overlap (int) – amount of points beyond the overlap to use to compute the shift.

seismic.monitor.monitor.average_components_mem_save(dvs: Generator[DV, None, None], save_scatter: bool = True) → DV[source]#

Averages the Similariy matrix of the DV objects. Based on those, it computes a new dv value and a new correlation value. Less memory intense but slower than average_components(). Should be used if you run out of memory using the former. Uses the sum of squares rule for the standard deviation.

Parameters:

dvs (Iterator[class:~seismic.monitor.dv.DV]) – Iterator over dvs from the different components to compute an average from. Note that it is possible to use almost anything as input as long as the similarity matrices of the dvs have the same shape
save_scatter (bool, optional) – Saves the scattering of old values and correlation coefficients to later provide a statistical measure. Defaults to True.

Raises:

TypeError – for DVs that were computed with different methods

Returns:

A single dv with an averaged similarity matrix.

Return type:

seismic.monitor.monitor.average_components(dvs: List[DV], save_scatter: bool = True, correct_shift: bool = False, correct_shift_method: str = 'mean', correct_shift_overlap: int = 0) → DV[source]#

Averages the Similariy matrix of the DV objects. Based on those, it computes a new dv value and a new correlation value.

Parameters:

dvs (List[class:~seismic.monitor.dv.DV]) – List of dvs from the different components to compute an average from. Note that it is possible to use almost anything as input as long as the similarity matrices of the dvs have the same shape
save_scatter (bool, optional) – Saves the scattering of old values and correlation coefficients to later provide a statistical measure. Defaults to True.
correct_shift – Shift the dvs with respect to each other. This is relevant if not all stations are active all the time and the references are not describing the same state. This only works if stations were active almost simultaneously at least for a bit (i.e., corr-start[n+n_overlap]-corr-start[n])

Raises:

TypeError – for DVs that were computed with different methods

Returns:

A single dv with an averaged similarity matrix.

Return type:

See also

If you should get an OutOfMemoryError consider using average_components_mem_save().

seismic.monitor.monitor.average_dvs_by_coords(dvs: List[DV], lat: Tuple[float, float], lon: Tuple[float, float], el: Tuple[float, float] = (-1000000.0, 1000000.0), return_std: bool = False) → DV[source]#

Averages the Similariy matrix of the DV objects if they are inside of the queried coordionates. Based on those, it computes a new dv value and a new correlation value.

Parameters:

dvs (List[class:~seismic.monitor.dv.DV]) – List of dvs from the different components to compute an average from. Note that it is possible to use almost anything as input as long as the similarity matrices of the dvs have the same shape
lat (Tuple[float, float]) – Minimum and maximum latitude of the stations(s). Note that for cross-correlations both stations must be inside of the window.
lon (Tuple[float, float]) – Minimum and maximum longitude of the stations(s). Note that for cross-correlations both stations must be inside of the window.
el (Tuple[float, float], optional) – Minimum and maximum elevation of the stations(s) in m. Note that for cross-correlations both stations must be inside of the window. Defaults to (-1e6, 1e6) - i.e., no filter.
return_std (bool, optional) – Save the standard deviation of the similarity matrices into the dv object. Defaults to False.

Raises:

TypeError – for DVs that were computed with different methods

Returns:

A single dv with an averaged similarity matrix.

Return type:

seismic.monitor.monitor.average_components_wfc(wfcs: List[WFC]) → WFC[source]#

Averages the Similariy matrix of the three DV objects. Based on those, it computes a new dv value and a new correlation value.

Parameters:: dvs (List[class:~seismic.monitor.wfc.WFC]) – List of dvs from the different components to compute an average from. Note that it is possible to use almost anything as input (also from different stations). However, at the time, the function requires the input to be of the same shape
Raises:: TypeError – for DVs that were computed with different methods
Returns:: A single dv with an averaged similarity matrix
Return type:: DV

seismic.monitor.dv#

Module to handle velocity change objects.

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Tuesday, 15th June 2021 04:12:18 pm Last Modified: Wednesday, 19th June 2024 04:07:10 pm

class seismic.monitor.dv.DV(corr: ndarray, value: ndarray, value_type: str, sim_mat: ndarray, second_axis: ndarray, method: str, stats: CorrStats, stretches: ndarray = None, corrs: ndarray = None, n_stat: ndarray = None, dv_processing: dict = {})[source]#

Bases: object

A simple object to save velocity changes.

Creates an object designed to hold and process velocity changes.

Parameters:

corr (np.ndarray) – Array with maximum Correlation coefficients.
value (np.ndarray) – Array with decimal velocity changes / stretches.
value_type (str) – Type of value (e.g., stretch)
sim_mat (np.ndarray) – Similarity Matrix with the dimensions time and stretch.
second_axis (np.ndarray) – second axis (e.g., the stretch axis)
method (str) – Stretching method that was used
stats (CorrStats) – Stats of the correlation object that was used
stretches (np.ndarray, optional) – Scatter of the value (e.g., stretch). In case this dv object holds an average of several dvs, the user can choose to keep all dvs’ stretches. defaults to None
corrs (np.ndarray, optional) – Scatter of the correlation coefficient (e.g., stretch). In case this dv object holds an average of several dvs, the user can choose to keep all dvs’ corrs. defaults to None
n_stat (np.ndarray, optional) – Number of stations used for the stack at corr_start t. Has the same shape as value and corr. Defaults to None

save(path: str)[source]#

Saves the dv object to a compressed numpy binary format. The DV can later be read using read_dv().

Parameters:: path (str) – output path

plot(save_dir: str = '.', figure_file_name: str = None, mark_time: datetime = None, normalize_simmat: bool = False, sim_mat_Clim: List[float] = [], xlim: Tuple[datetime, datetime] = None, ylim: Tuple[int, int] = None, plot_scatter: bool = False, figsize: Tuple[float, float] = (9, 11), dpi: int = 144, title: str = None, return_ax=False, style: str = 'technical', dateformat: str = '%d %b %y', ax: Axes = None) → Tuple[figure, List[axis]][source]#

Plots the dv object into a multi-panel-view of similarity matrix dv-value, and correlation coefficient.

Parameters:

save_dir (str, optional) – Directory to save the figure to, defaults to ‘.’
figure_file_name (str, optional) – Filename of the figure. If None, the figure will not be save, defaults to None.
mark_time (datetime, optional) – times to put ticks to. If not set, they will be determined automatically, defaults to None.
normalize_simmat (bool, optional) – Normalise the similarity matrix? Defaults to False.
sim_mat_Clim (List[float], optional) – Color limits for the similarity matrix , defaults to []. Format is [low_limit, upper_limit]
xlim (Tuple[datetime, datetime], optional) – Time Limits to plot, defaults to None.
ylim (Tuple[int, int], optional) – Y-limits (e.g., stretch) of the plot, defaults to None
plot_scatter (bool, optional) – If set to True, the scatter of all used combinations is plotted in form of a histogram. Defaults to False.
figsize (Tuple[float, float], optional) – Size of the figure/canvas, defaults to (9, 11)
dpi (int, optional) – Pixels per inch, defaults to 144
title (str, optional) – Define a custom title. Otherwise, an automatic title will be created, defaults to None
return_ax (bool, optional) – Return plt.figure and list of axes. Defaults to False. This overwrites any choice to save the figure.
style (str, optional) – Style of the plot. Defaults to ‘technical’. The other option would be ‘publication’ (looks better but is less informative).
dateformat (str, optional) – Format of the date on the x-axis. Defaults to ‘%d %b %y’.

Returns:

If return_ax is set to True it returns fig and axes.

smooth_sim_mat(win_len: int, exclude_corr_below: Optional[float] = None)[source]#

Smoothes the similarity matrix along the correlation time axis.

Parameters:

win_len (int) – Length of the window in number of samples.
exclude_corr_below (float) – Exclude data points with a correlation coefficient below the chosen threshold. Defaults to None (i.e., include all)

Note::

This action is perfomed in-place and irreversible as the original data are not accesible any more.

seismic.monitor.dv.read_dv(path: str) → DV[source]#

Reads a saved velocity change object from an .npz file.

Parameters:: path (str) – Path to file, wildcards allowed
Returns:: the corresponding and converted DV object
Return type:: DV, List[DV] if path contains wildcards

seismic.monitor.spatial#

Module to invert velocity change time series from a set of correlations to a spatial velocity field. Base on the algorithm proposed in [Obermann et al. (2013)](https://doi.org/10.1002/2013JB010399).

Spatial imaging by means of a Boltzmann-equation-based sensitivity kernel (see Obermann et al., 2013 and Paasschens, 1997)

Implementation here is just for the 2D case

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Monday, 16th January 2023 10:53:31 am Last Modified: Tuesday, 5th March 2024 03:03:01 pm

seismic.monitor.spatial.probability(dist: numpy.ndarray | float, t: numpy.ndarray | float, vel: float, mf_path: float, atol: float) → numpy.ndarray | float[source]#

Compute the probability for a wave with the given properties to pass through a point at distance dist. Following the diffusion law.

Implementation as in Obermann et al. (2013).

Parameters:

dist (np.ndarray | float) – distance to the source. Can be a matrix.
t (float or 1D numpy array) – time
vel (float) – velocity for homogeneous velocities
mf_path (float) – Mean free path of the wave.

Returns:

The probability. In the same shape as dist.

Return type:

np.ndarray | float

seismic.monitor.spatial.compute_grid_dist(x: ndarray, y: ndarray, x0: float, y0: float) → ndarray[source]#

Compute the distance of each point on the grid x, y to x0, y0. (2D)

Parameters:

x (np.ndarray) – X coordinate axis
y (np.ndarray) – y coordiante axis
x0 (float) – x-position of point
y0 (float) – y-position of point

Returns:

distance matrix in form of a 2D matrix.

Return type:

np.ndarray

seismic.monitor.spatial.sensitivity_kernel(s1: ndarray, s2: ndarray, x: ndarray, y: ndarray, t: float, dt: float, vel: float, mf_path: float) → ndarray[source]#

Computes a 2D surface-wave sensitivity kernel for an ambient noise cross-correlation between two station s1 and s2. The computation is based on a time-dependendent solution of the Boltzmann equation for a homogeneous medium (see Obermann et al., 2013 and Paasschens, 1997).

Implementation as in Obermann et al. (2013). Limits for nan and inf will return 0, or 1 respectively.

Note

A large dt can lead to artefacts in the kernel. dt should always be smaller than dx/(2*vel).

Parameters:

s1 (np.ndarray) – Position of Station 1, format is [x, y].
s2 (np.ndarray) – Position of Station 2, format is [x, y].
x (np.ndarray) – Monotonously increasing array holding the x-coordinates of the grid the kernel should be computed on.
y (np.ndarray) – Monotonously increasing array holding the y-coordinates of the grid the kernel should be computed on.
t (float) – Time of coda that is probed. Usually set as the middle of the used time window
dt (float) – Sampling interval that the probability will be computed in. Note that dt should not be too large. Otherwise, numerical artefacts can occur.
vel (float) – Surface wave velocity. Note that units between coordinate grid and velocity have to be consistent (i.e., km/s and km or m and m/s)
mf_path (float) – Mean free path of the wave (see Obermann et al. 2013 for details)

Returns:

The sensitivity_kernel of the grid. Can be plotted against. np.meshgrid(x, y)

Return type:

np.ndarray

seismic.monitor.spatial.data_variance(corr: numpy.ndarray | float, bandwidth: float, tw: Tuple[float, float], freq_c: float) → numpy.ndarray | float[source]#

Compute the variance of the velocity change data based on the coherence value between stretched reference trace and correlation function.

(see Obermann et. al. 2013)

Parameters:

corr (np.ndarray | float) – coherence
bandwidth (float) – Bandwidth between low and highpass filter
tw (Tuple[float, float]) – Time Window used for the dv/v computation. In the form Tuple[tw_start, tw_end]
freq_c (float) – Centre frequency.

Returns:

Data Variance in the same shape as corr

Return type:

np.ndarray | float

seismic.monitor.spatial.compute_cm(scaling_factor: float, corr_len: float, std_model: float, dist: numpy.ndarray | float) → numpy.ndarray | float[source]#

Computes the model variance for the dv/v grid.

(see Obermann et. al. 2013)

Parameters:

scaling_factor (float) – Scaling factor, dependent on the cell length. (e.g., 1*cell_length)
corr_len (float) – Length over which the parameters are related (high value means stronger smoothing).
std_model (float) – A-priori standard deviation of the model. Corresponds to the best agreement between 1) the stability of the model for the velocity variations and 2) the minimised difference between the model predictions (output of the forward problem) and the data. Can be estimated via the L-curve criterion (see Hansen 1992)
dist (np.ndarray) – Matrix containing distance between the two cells i and j. Will always be 0 on the diagonal where i==j.

Returns:

Returns the model variance in the same shape as dist. I.e., NxN where N is the number of station combinations.

Return type:

np.ndarray

seismic.monitor.spatial.geo2cart(lat: numpy.ndarray | float, lon: numpy.ndarray | float, lat0: float, lon0: float) → Tuple[numpy.ndarray | float, numpy.ndarray | float][source]#

Convert geographic coordinates to Northing and Easting for a local grid.

Parameters:

lat (np.ndarray | float) – Latitude(s)
lon (np.ndarray | float) – Longitude(s)
lat0 (float) – Latitude of the lower left corner of the grid
lon0 (float) – Latitude of the lower left corner of the grid

Returns:

A tuple holding (Easting, Northing) in km and in the same shape as lat and lon.

Return type:

Tuple[np.ndarray, np.ndarray]

class seismic.monitor.spatial.DVGrid(lat0: float, lon0: float, res: float, x: float, y: float, dt: float, vel: float, mf_path: float)[source]#

Bases: object

Object to handle spatially gridded DV objects.

Initialises the grid to compute dv/v on. The grid will be in Cartesian coordinates with the lower left boarder (lon0, lat0) as reference point.

Parameters:

lat0 (float) – Minimum latitude on the grid.
lon0 (float) – Minimum longitude on the grid
res (float) – Resolution in km
x (float) – Eastwards extent of the grid in km
y (float) – Northwards extent of the grid in km
dt (float) – Sampling interval for the numerical integration of the sensitivity kernels. Coarse spacing might cause artefacts. Fine spacing will result in higher computation times.
vel (float) – Wave velocity
mf_path (float) – Mean free path of a wave

forward_model(dv_grid: ndarray, dvs: Optional[Iterable[DV]] = None, utc: Optional[UTCDateTime] = None, tw: Optional[Tuple[float, float]] = None, stat0: Optional[Iterable[Tuple[float, float]]] = None, stat1: Optional[Iterable[Tuple[float, float]]] = None, verbose: bool = False) → ndarray[source]#

Solves the forward problem associated to the grid computation. I.e., computes the velocity changes theoretically measured at each station for a given dv/v grid.

Parameters:

dv_grid (np.ndarray) – dv/v grid. Must have the same shape as self.xgrid
dvs (Optional[Iterable[DV]], optional) – Iterator over DV to extract the parameters from. If this is defined utc has to be defined as well. If it’s None all other parameters have to be defined, defaults to None
utc (Optional[UTCDateTime], optional) – The time to extract the properties from the dv files from, defaults to None
tw (Optional[Tuple[float, float]], optional) – Lapse time Window used for the dv/v estimation, defaults to None
stat0 (Optional[Iterable[Tuple[float, float]]], optional) – Coordinates of the first station of the cross-correlation. Given as an array with the form: [(latitude0, longitude0), (latitude1, longitude1), …], defaults to None.
stat1 (Optional[Iterable[Tuple[float, float]]], optional) – Coordinates of the second station of the cross-correlation. Given as an array with the form: [(latitude0, longitude0), (latitude1, longitude1), …], defaults to None., defaults to None

Raises:

TypeError – If too few variables are defined

Returns:

A vector with one dv/v value for each station combination in the stat0-stat1 vectors.

Return type:

np.ndarray

compute_dv_grid(dvs: Iterator[DV], utc: UTCDateTime, scaling_factor: float, corr_len: float, std_model: float, tw: Optional[Tuple[float, float]] = None, freq0: Optional[float] = None, freq1: Optional[float] = None, compute_resolution: bool = False, verbose: bool = False) → ndarray[source]#

Perform a linear least-squares inversion of station specific velocity changes (dv/v) at time t.

The algorithm will perform the following steps: 1. Computes the sensitivity kernels for each station pair. 2. Computes model and data variance. 3. Inverts the velocity changes onto a flattened grid.

Note

All computations are performed on a flattened grid, so that the grid is actually just one-dimensional. The return value is reshaped onto the shape of self.xgrid.

Note

The returned value and the value in self.vel_change is strictly speaking not a velocity change, but an epsilon value. That corresponds to -dv/v

Parameters:

dvs (Iterator[DV]) – The :class:~seismic.monitor.dv.DV objects to use. Note that the input can be an iterator, so that not actually all dvs have to be loaded into RAM at the same time.
utc (UTCDateTime) – Time at which the velocity changes should be inverted. Has to be covered by the DV objects.
scaling_factor (float) – Scaling factor, dependent on the cell length. (e.g., 1*cell_length)
corr_len (float) – Length over which the parameters are related (high value means stronger smoothing).
std_model (float) – A-priori standard deviation of the model. Corresponds to the best agreement between 1) the stability of the model for the velocity variations and 2) the minimised difference between the model predictions (output of the forward problem) and the data. Can be estimated via the L-curve criterion (see Hansen 1992)
tw (Tuple[float, float], optional) – Time window used to evaluate the velocity changes. Given in the form Tuple[tw_start, tw_end]. Not required if the used dvs have a dv_processing attribute.
freq0 (float, optional) – high-pass frequency of the used bandpass filter. Not required if the used dvs have a dv_processing attribute.
freq1 (float, optional) – low-pass frequency of the used bandpass filter. Not required the used dvs have a dv_processing attribute.
compute_resolution (bool, defaults to False) – Compute the resolution matrix? Will be saved in self.resolution

Raises:

ValueError – For times that are not covered by the dvs time-series.

Returns:

The dv-grid (with corresponding coordinates in self.xgrid and self.ygrid or self.xaxis and self.yaxis, respectively). Note this is actually an epsilon grid (meaning it represents the stretch value i.e., -dv/v)

Return type:

np.ndarray

compute_dv_tsvd(dvs: Iterator[DV], utc: UTCDateTime, tw: Optional[Tuple[float, float]] = None, freq0: Optional[float] = None, freq1: Optional[float] = None, cutoff: float = 0.2, verbose: bool = False) → ndarray[source]#

Performs a truncated singular value decomposition (TSVD) to invert for the dv/v grid. The advantage of this approach is that no a-priori information and no subjective damping parameter has to be defined. However, the disadvantage is that the inversion is not aware of the data variance (defined by the correlation coefficient). This can be a useful approach for dv estimates with a very high stability.

The algorithm will perform the following steps: 1. Computes the sensitivity kernels for each station pair. 2. Computes model and data variance. 3. Inverts the velocity changes onto a flattened grid.

Note

All computations are performed on a flattened grid, so that the grid is actually just one-dimensional. The return value is reshaped onto the shape of self.xgrid.

Note

The returned value and the value in self.vel_change is strictly speaking not a velocity change, but an epsilon value. That corresponds to -dv/v

Warning

This has not been thoroughly tested and has the disadvantage that no “physical” prior weighting is done (as opposed to the linear-least-squares method), where data weighting is performed based on the coherence values.

Parameters:

dvs (Iterator[DV]) – The :class:~seismic.monitor.dv.DV objects to use. Note that the input can be an iterator, so that not actually all dvs have to be loaded into RAM at the same time.
utc (UTCDateTime) – Time at which the velocity changes should be inverted. Has to be covered by the DV objects.
tw (Tuple[float, float], optional) – Time window used to evaluate the velocity changes. Given in the form Tuple[tw_start, tw_end]. Not required if the used dvs have a dv_processing attribute.
freq0 (float, optional) – high-pass frequency of the used bandpass filter. Not required if the used dvs have a dv_processing attribute.
freq1 (float, optional) – low-pass frequency of the used bandpass filter. Not required the used dvs have a dv_processing attribute.
compute_resolution (bool, defaults to False) – Compute the resolution matrix? Will be saved in self.resolution
cutoff – Cutoff value for the singular values. All singular values below this fraction of the maximum singular value will be set to 0.

Raises:

ValueError – For times that are not covered by the dvs time-series.

Returns:

Return type:

np.ndarray

compute_resolution(dvs: Iterable[DV], utc: UTCDateTime, scaling_factor: float, corr_len: float, std_model: float, tw: Optional[Tuple[float, float]] = None, freq0: Optional[float] = None, freq1: Optional[float] = None, verbose: bool = False) → ndarray[source]#

Compute the model resolution of the dv/v model.

See also

Described in the supplement to Obermann (2013)

Note

Note that it’s significantly faster to perform a joint inversion or a gridded dv/v and resolution rather than performing them separately.

Parameters:

dvs (Iterator[DV]) – The :class:~seismic.monitor.dv.DV objects to use. Note that the input can be an iterator, so that not actually all dvs have to be loaded into RAM at the same time.
utc (UTCDateTime) – Time at which the velocity changes should be inverted. Has to be covered by the DV objects.
scaling_factor (float) – Scaling factor, dependent on the cell length. (e.g., 1*cell_length)
corr_len (float) – Length over which the parameters are related (high value means stronger smoothing).
std_model (float) – A-priori standard deviation of the model. Corresponds to the best agreement between 1) the stability of the model for the velocity variations and 2) the minimised difference between the model predictions (output of the forward problem) and the data. Can be estimated via the L-curve criterion (see Hansen 1992)
tw (Tuple[float, float], optional) – Time window used to evaluate the velocity changes. Given in the form Tuple[tw_start, tw_end]. Not required if dv has attribute dv_processing.
freq0 (float, optional) – high-pass frequency of the used bandpass filter. Not required if dv has attribute dv_processing.
freq1 (float, optional) – low-pass frequency of the used bandpass filter. Not required if dv has attribute dv_processing.

Raises:

ValueError – For times that are not covered by the dvs time-series.

Returns:

A gridded resolution matrix.

Return type:

np.ndarray

compute_posterior_covariance(dvs: Iterable[DV], utc: UTCDateTime, scaling_factor: float, corr_len: float, std_model: float, tw: Optional[Tuple[float, float]] = None, freq0: Optional[float] = None, freq1: Optional[float] = None, verbose: bool = False) → ndarray[source]#

Computes the posterior covariance matrix of the dv/v model.

Parameters:

dvs (Iterator[DV]) – The :class:~seismic.monitor.dv.DV objects to use. Note that the input can be an iterator, so that not actually all dvs have to be loaded into RAM at the same time.
utc (UTCDateTime) – Time at which the velocity changes should be inverted. Has to be covered by the DV objects.
scaling_factor (float) – Scaling factor, dependent on the cell length. (e.g., 1*cell_length)
corr_len (float) – Length over which the parameters are related (high value means stronger smoothing).
std_model (float) – A-priori standard deviation of the model. Corresponds to the best agreement between 1) the stability of the model for the velocity variations and 2) the minimised difference between the model predictions (output of the forward problem) and the data. Can be estimated via the L-curve criterion (see Hansen 1992)
tw (Tuple[float, float], optional) – Time window used to evaluate the velocity changes. Given in the form Tuple[tw_start, tw_end]. Not required if dv has attribute dv_processing.
freq0 (float, optional) – high-pass frequency of the used bandpass filter. Not required if dv has attribute dv_processing.
freq1 (float, optional) – low-pass frequency of the used bandpass filter. Not required if dv has attribute dv_processing.

Raises:

ValueError – For times that are not covered by the dvs time-series.

Returns:

A gridded resolution matrix.

Return type:

np.ndarray

plot(plot_stations: bool = True, cmap: str = 'seismic_r', ax: Optional[Axis] = None, *args, **kwargs) → Axis[source]#

Simple heatplot of the dv/v grid in Cartesian coordinates.

Parameters:

plot_stations (bool, optional) – plot station coordinates, defaults to True
cmap (str, optional) – Colormap to use, defaults to ‘seismic’
ax (Optional[mpl.axis.Axis], optional) – Axis to plot into, defaults to None

Returns:

Returns the axis handle

Return type:

mpl.axis.Axis

align_dvs_to_grid(dvs: List[DV], utc: UTCDateTime, steps: int, corr_thres: float, save: bool | str = False) → List[DV][source]#

Align dv/v curves to the forward model at a given time.

Parameters:

dvs (List[DV]) – List of dv/v objects to align
utc (UTCDateTime) – Time at which the dv/v curve should be aligned
steps (int) – Number of steps to use for the alignment
corr_thres (float) – Correlation threshold
save (bool | str, optional) – If given, save the aligned dv/v curves to the given directory. Defaults to False.

Returns:

List of dv/v curves that have been aligned

Return type:

List[DV]

seismic.monitor.spatial.align_dv_curves(dv: DV, utc: UTCDateTime, steps: int, value: float = 0.0)[source]#

Align a dv/v curve to a given value at a given time.

Parameters:

dv (DV) – dv/v object to align
utc (UTCDateTime) – Time at which the dv/v curve should be aligned
steps (int) – Number of steps to use for the alignment
value (float, optional) – Value to align the dv/v curve to, defaults to 0.

seismic.monitor.spatial.dv_starts(dv: DV, utc: UTCDateTime, corr_thres: float) → bool[source]#

Check if a dv/v curve is available at a given time and if it has a correlation value above a given threshold.

Parameters:

dv (DV) – dv/v object to check
utc (UTCDateTime) – Time at which the dv/v curve should be aligned
corr_thres (float) – Correlation threshold

Returns:

True if the dv/v curve is available and has a correlation value above the threshold

Return type:

bool

seismic.monitor.wfc#

Module to analyse the waveform coherence.

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Friday, 5th November 2021 08:19:58 am Last Modified: Tuesday, 11th July 2023 04:42:52 pm

class seismic.monitor.wfc.WFC(wfc_dict: dict, stats: CorrStats, wfc_processing: dict)[source]#

Bases: dict

Object to hold the properties of a waveform coherency for a single frequency and lapse time window.

The waveform coherency can be used as a measure of stability of a certain correlation. See Steinmann, et. al. (2021) for details.

Initialise WFC object.

Parameters:

wfc_dict (dict) – Dictionary holding the correlation values for each of the reference traces. The keys follow the syntax reftr_%n where n is the number of the reference trace (starting with 0).
stats (CorrStats) – CorrStats object from the original CorrBulk. Also has to contain the keys holding the information about the bandpass filter and the used lapse time window.

compute_average()[source]#

Computes the average of the waveform coherency over the time axis.

Parameters:: average_reftrcs (bool, optional) – Average over all reference traces. Defaults to True.

save(path: str)[source]#

Saves the object to an npz file.

Parameters:: path (str) – Path to save to.

class seismic.monitor.wfc.WFCBulk(wfcl: List[WFC])[source]#

Bases: object

Object to hold the waveform coherency for one station and subsequently plot it against frequency and lapse time window.

plot(title: str = None, log: bool = False, cmap: str = 'viridis') → Axes[source]#

Create a plot of the waveform coherency against frequency and lapse time window.

Parameters:

title (str, optional) – Title of the plot, defaults to None
log (bool, optional) – plot on logarithmic frequency axis, defaults to False

seismic.monitor.wfc.read_wfc(path: str) → WFC[source]#

Read a WFC object saved as npz. Returns a list if wildcards are used

Parameters:: path (str) – Path to file
Returns:: The WFC object.
Return type:: WFC

seismic.monitor.post_corr_process#

Holds the postprocessing functions for CorrBulks.

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Christoph Sens-Schönefelder Peter Makus (makus@gfz-potsdam.de)

Created: Monday, 14th June 2021 08:50:57 am Last Modified: Wednesday, 4th October 2023 09:29:49 am

seismic.monitor.post_corr_process.corr_mat_clip(A: ndarray, thres: float, axis: int) → ndarray[source]#

Clip the Input Matrix upper and lower bounds to a multiple of its standard deviation. thres determines the factor and axis the axis the std should be computed over

Parameters:

A (np.ndarray) – Input Array to be clipped
thres (float) – factor of the standard deviation to clip by
axis (int) – Axis to compute the std over and, subsequently clip over. Can be None, if you wish to compute floating point rather than a vector. Then, the array will be clipped evenly.

Returns:

The clipped array

Return type:

np.ndarray

seismic.monitor.post_corr_process.corr_mat_smooth(data: ndarray, wsize: int, wtype: str = 'flat', axis: int = 1) → ndarray[source]#

Smoothing of a correlation matrix.

Parameters:

corr_mat (dictionary of the type correlation matrix) – correlation matrix to be smoothed
wsize (int) – Window size
wtype (string) – Window type. It can be one of: [‘flat’, ‘hanning’, ‘hamming’, ‘bartlett’, ‘blackman’] defaults to ‘flat’
axis (int) – Axis along with apply the filter. O: smooth along correlation lag time axis 1: smooth along time axis

Return type:

ndarray

Returns:

X: Filtered matrix

seismic.monitor.post_corr_process.unicode_to_string(input)[source]#: Convert all unicode strings to utf-8 strings

seismic.monitor.post_corr_process.corr_mat_filter(data: ndarray, stats: CorrStats, freqs: Tuple[float, float], order=3) → ndarray[source]#

Filter a correlation matrix.

Filters the correlation matrix corr_mat in the frequency band specified in freqs using a zero phase filter of twice the order given in order.

Parameters:

data (np.ndarray) – correlation matrix
stats (CorrStats) – The stats object corresponding to the CorrBulk object
freqs (array-like of length 2) – lower and upper limits of the pass band in Hertz
order (int) – half the order of the Butterworth filter

Return type:

np.ndarray

Returns:

filtered correlation matrix

seismic.monitor.post_corr_process.corr_mat_resample(data: ndarray, stats: CorrStats, start_times: List[UTCDateTime], end_times=[]) → Tuple[ndarray, CorrStats][source]#

Function to create correlation matrices with constant sampling

When created with recombine_corr_data the correlation matrix contains all available correlation streams but homogeneous sampling is not guaranteed as correlation functions may be missing due to data gaps. This function restructures the correlation matrix by inserting or averaging correlation functions to provide temporally homogeneous sampling. Inserted correlation functions consist of ‘nan’ if gaps are present and averaging is done if more than one correlation function falls in a bin between start_times[i] and end_times[i]. If end_time is an empty list (default) end_times[i] is set to start_times[i] + (start_times[1] - start_times[0])

Parameters:

data (np.ndarray) – 2D matrix holding the correlation data
Stats (Stats) – The stats object belonging to the CorrBulk object.
start_times (list of UTCDateTime objects) – list of starting times for the bins of the new sampling
end_times (list of UTCDateTime objects) – list of starting times for the bins of the new sampling

Return type:

tuple

Returns:

the new data array and altered stats object.

seismic.monitor.post_corr_process.corr_mat_correct_decay(data: ndarray, stats: CorrStats) → ndarray[source]#

Correct for the amplitude decay in a correlation matrix.

Parameters:

data (np.ndarray) – Matrix holding correlation data
stats (CorrStats) – Stats object afiliated to a CorrBulk object

Returns:

correlations with decay corrected amplitudes

Return type:

np.ndarray

seismic.monitor.post_corr_process.corr_mat_envelope(data: ndarray) → ndarray[source]#

Calculate the envelope of a correlation matrix.

The corrlation data of the correlation matrix are replaced by their Hilbert envelopes.

Parameters:: data (np.ndarray) – correlation matrix.
Return type:: np.ndarray
Returns:: hilb: The envelope of the correlation data in the same shape as the input array.

seismic.monitor.post_corr_process.corr_mat_normalize(data: ndarray, stats: CorrStats, starttime: float = None, endtime: float = None, normtype: str = 'energy') → ndarray[source]#

Correct amplitude variations with time in a correlation matrix.

Parameters:

data (np.ndarray) – correlation matrix from CorrBulk
stats (CorrStats) – The stats object from the CorrBulk object.
starttime (float) – beginning of time window in seconds with respect to the zero position
endtime (float) – end time window in seconds with respect to the zero position
normtype (string) – one of the following ‘energy’, ‘max’, ‘absmax’, ‘abssum’ to decide about the way to calculate the normalization.

Return type:

np.ndarray

Returns:

Correlation matrix with normalized ampitudes.

seismic.monitor.post_corr_process.corr_mat_mirror(data: ndarray, stats: CorrStats) → Tuple[ndarray, CorrStats][source]#

Average the causal and acausal parts of a correlation matrix.

Parameters:

data (np.ndarray) – A matrix holding correlation data
stats (CorrStats) – the stats object corresponding to a CorrBulk object.

Returns:

as the input but with avaraged causal and acausal parts of the correlation data

Return type:

Tuple[np.ndarray, CorrStats]

seismic.monitor.post_corr_process.corr_mat_taper(data: ndarray, stats: CorrStats, width: float) → ndarray[source]#

Taper a correlation matrix.

Apply a taper to all traces in the correlation matrix.

Parameters:

data (np.ndarray) – correlation matrix`
width (float) – width to be tapered in seconds (per side)

Return type:

np.ndarray

Returns:

The tapered matrix

..note:: In-place operation

seismic.monitor.post_corr_process.corr_mat_taper_center(data: ndarray, stats: CorrStats, width: float, slope_frac: float = 0.05) → ndarray[source]#

Taper the central part of a correlation matrix.

Parameters:

data (np.ndarray) – Correlation data from CorrBulk object
stats (CorrStats) – Stats object from CorrBulk
width (float) – width of the central window to be tapered in seconds (in total, i.e. not per taper side)
slope_frac (float, optional) – fraction of width used for soothing of edges, defaults to 0.05

Returns:

The tapered matrix

Return type:

np.ndarray

seismic.monitor.post_corr_process.corr_mat_resample_time(data: ndarray, stats: CorrStats, freq: float) → Tuple[ndarray, CorrStats][source]#

Resample the lapse time axis of a correlation matrix. The correlations are automatically filtered with a highpass filter of 0.4*sampling frequency to avoid aliasing.

Parameters:

data (np.ndarray) – correlation matrix
stats (CorrStats) – The stats object associated to the CorrBulk object.
freq (float) – new sampling frequency

Return type:

Tuple[np.ndarray, CorrStats]

Returns:

Resampled data and changed stats object.

..note:: This action is performed in-place.

seismic.monitor.post_corr_process.corr_mat_decimate(data: ndarray, stats: CorrStats, factor: int) → Tuple[ndarray, CorrStats][source]#

Downsample a correlation matrix by an integer sample. A low-pass filter is applied before decimating.

Parameters:

data (np.ndarray) – correlation matrix
stats (CorrStats) – The stats object associated to the CorrBulk object.
factor (int) – The factor to reduce the sampling frequency by

Return type:

Tuple[np.ndarray, CorrStats]

Returns:

Decimated data and changed stats object.

..note:: This action is performed in-place.

seismic.monitor.post_corr_process.corr_mat_resample_or_decimate(data: ndarray, stats: CorrStats, freq: float) → Tuple[ndarray, CorrStats][source]#

Downsample a correlation matrix. Decides automatically whether to decimate or downsampled depneding on the target frequency. A low-pass filter is applied to avoid aliasing.

Parameters:

data (np.ndarray) – correlation matrix
stats (CorrStats) – The stats object associated to the CorrBulk object.
freq (float) – new sampling frequency

Return type:

Tuple[np.ndarray, CorrStats]

Returns:

Decimated data and changed stats object.

..note:: This action is performed in-place.

seismic.monitor.post_corr_process.corr_mat_extract_trace(data: ndarray, stats: CorrStats, method: str = 'mean', percentile: float = 50.0) → ndarray[source]#

Extract a representative trace from a correlation matrix.

Extract a correlation trace from the that best represents the correlation matrix. Method decides about method to extract the trace. The following possibilities are available

mean averages all traces in the matrix
median extracts the median trace
norm_mean averages the traces normalized after normalizing for maxima
similarity_percentile averages the percentile % of traces that
best correlate with the mean of all traces. This will exclude abnormal traces. percentile = 50 will return an average of traces with correlation (with mean trace) above the median.

Parameters:

corr_mat (dictionary) – correlation matrix dictionary
method (string) – method to extract the trace
percentile (float) – only used for method==’similarity_percentile’

Return type:

trace dictionary of type correlation trace

Return tracetrace:

extracted trace

exception seismic.monitor.post_corr_process.Error[source]#: Bases: Exception

exception seismic.monitor.post_corr_process.InputError(msg)[source]#: Bases: Error

seismic.monitor.post_corr_process.corr_mat_stretch(cdata: ndarray, stats: CorrStats, ref_trc: ndarray = None, tw: List[ndarray] = None, stretch_range: float = 0.1, stretch_steps: int = 100, sides: str = 'both', return_sim_mat: bool = False, ref_tr_trim: Optional[Tuple[float, float]] = None, ref_tr_stats: Optional[CorrStats] = None) → dict[source]#

Time stretch estimate through stretch and comparison.

This function estimates stretching of the time axis of traces as it can occur if the propagation velocity changes.

Time stretching is estimated comparing each correlation function stored in the corr_data matrix (one for each row) with stretch_steps stretched versions of reference trace stored in ref_trc. If ref_trc is None the mean of all traces is used. The maximum amount of stretching may be passed in stretch_range. The time axis is multiplied by 1/(1 + stretch). The best match (stretching amount and corresponding correlation value) is calculated on different time windows. If tw = None the stretching is estimated on the whole trace.

Parameters:

corr_mat (dictionary) – correlation matrix dictionary as produced by recombine_corr_data
ref_trc (ndarray) – 1D array containing the reference trace to be stretched and compared to the individual traces in mat
tw (list of ndarray of int) – list of 1D ndarrays holding the indices of samples in the time windows to be use in the time shift estimate. The samples are counted from the zero lag time with the index of the first sample being 0. If tw = None the full time range is used.
stretch_range (scalar) – Maximum amount of relative stretching. Stretching and compression is tested from -stretch_range to stretch_range.
stretch_steps (scalar`) – Number of shifted version to be tested. The increment will be (2 * stretch_range) / stretch_steps
sides (str) – Side of the reference matrix to be used for the stretching estimate (‘both’ | ‘left’ | ‘right’ | ‘single’) single is used for one-sided signals from active sources with zero lag time is on the first sample. Other options assume that the zero lag time is in the center of the traces.

Return type:

Dictionary

Returns:

dv: Dictionary with the following keys

corr: 2d ndarray containing the correlation value for the best: match for each row of mat and for each time window. Its dimension is: :func:(len(tw),mat.shape[1])
value: 2d ndarray containing the stretch amount corresponding to: the best match for each row of mat and for each time window. Stretch is a relative value corresponding to the negative relative velocity change -dv/v. Its dimension is: :func:(len(tw),mat.shape[1])
sim_mat: 3d ndarray containing the similarity matricies that: indicate the correlation coefficient with the reference for the different time windows, different times and different amount of stretching. Its dimension is: (len(tw),mat.shape[1],len(strvec))()
second_axis: It contains the stretch vector used for the velocity: change estimate.
vale_type: It is equal to ‘stretch’ and specify the content of: the returned ‘value’.
method: It is equal to ‘single_ref’ and specify in which “way” the: values have been obtained.

seismic.monitor.post_corr_process.corr_mat_shift(data: ndarray, stats: CorrStats, ref_trc: ndarray = None, tw: List[ndarray] = None, shift_range: float = 10, shift_steps: int = 101, sides: str = 'both', return_sim_mat: bool = False) → dict[source]#

Time shift estimate through shifting and comparison.

This function estimates shifting of the time axis of traces as it can occur if the clocks of digitizers drift.

Time shifts are estimated comparing each correlation function stored in the corr_data matrix (one for each row) with shift_steps shifted versions of reference trace stored in ref_trc. The maximum amount of shifting may be passed in shift_range. The best match (shifting amount and corresponding correlation value) is calculated on different time windows. If tw = None the shifting is estimated on the whole trace.

Parameters:

corr_mat (dictionary) – correlation matrix dictionary as produced by recombine_corr_data
ref_trc (ndarray) – 1D array containing the reference trace to be stretched and compared to the individual traces in mat
tw (list of ndarray of int) – list of 1D ndarrays holding the indices of sampels in the time windows to be use in the time shift estimate. The sampels are counted from the zero lag time with the index of the first sample being 0. If tw = None the full time range is used.
shift_range (scalar) – Maximum amount of shifting in samples. Shifting is tested from -stretch_range to stretch_range.
shift_steps (scalar`) – Number of shifted version to be tested. The increment will be (2 * shift_range) / shift_steps
sides (str) – Side of the reference matrix to be used for the shifting estimate (‘both’ | ‘right’) single is used for one-sided signals from active sources with zero lag time is on the first sample. Other options assume that the zero lag time is in the center of the traces.

Return type:

Dictionary

Returns:

dv: Dictionary with the following keys

corr: 2d ndarray containing the correlation value for the best: match for each row of mat and for each time window. Its dimension is: :func:(len(tw),mat.shape[1])
value: 2d ndarray containing the stretch amount corresponding to: the best match for each row of mat and for each time window. Stretch is a relative value corresponding to the negative relative velocity change -dv/v. Its dimension is: :func:(len(tw),mat.shape[1])
sim_mat: 3d ndarray containing the similarity matricies that: indicate the correlation coefficient with the reference for the different time windows, different times and different amount of stretching. Its dimension is: (len(tw),mat.shape[1],len(strvec))()
second_axis: It contains the stretch vector used for the velocity: change estimate.
vale_type: It is equal to ‘shift’ and specify the content of: the returned ‘value’.
method: It is equal to ‘single_ref’ and specify in which “way” the: values have been obtained.

seismic.monitor.post_corr_process.measure_shift(data: ndarray, stats: CorrStats, ref_trc: Optional[ndarray] = None, tw: Optional[List[List]] = None, shift_range: float = 10, shift_steps: int = 101, sides: str = 'both', return_sim_mat: bool = False) → List[DV][source]#

Time shift estimate through shifting and comparison.

This function estimates shifting of the time axis of traces as it can occur if the clocks of digitizers drift.

Time shifts are estimated comparing each trace (e.g. correlation function stored in the corr_data matrix (one for each row) with shifted versions of reference trace stored in ref_trc. The range of shifting to be tested is given in shift_range in seconds. It is used in a symmetric way from -shift_range to +``shift_range``. Shifting is tested shift_steps times. shift_steps should be an odd number to test zero shifting. The best match (shifting amount and corresponding correlation value) is calculated in specified time windows. Multiple time windows may be specified in tw.

Parameters:

data (ndarray) – correlation matrix or collection of traces. lag time runs along rows and number of traces is data.shape[0].
ref_trc (ndarray) – 1D array containing the reference trace to be shifted and compared to the individual traces in data. If ref_trc``is None, the average of all traces in ``data is used.
tw (list of lists with two floats) – list of lists that contain the start and end times of the time windows in which the similarity is to be estimated. Time are given in seconds of lag time with respect to zero lag time. Using multiple lists of start end end time the best shift can be estimated in multiple time windows. If tw = None one time window is used containing the full time range of lag times in the traces.
shift_range (scalar) – Maximum amount of shifting in second. Shifting is tested from -stretch_range to stretch_range.
shift_steps (scalar`) – Number of shifted version to be tested. The increment will be (2 * shift_range) / shift_steps
sides (str) – Side of the traces to be used for the shifting estimate (‘both’ | ‘single’). single is used for one-sided signals from active sources or if the time window shall not be symmetric. For both the time window will be mirrored about zero lag time, e.g. [start,end] will result in time windows [-end:-start] and [start:end] being used simultaneousy

Return type:

List[Dictionary]

Returns:

dv: List of len(tw) Dictionaries with the following keys

corr: 2d ndarray containing the correlation value for the best: match for each row of mat and for each time window. Its dimension is: :func:(len(tw),mat.shape[1])
value: 2d ndarray containing the stretch amount corresponding to: the best match for each row of mat and for each time window. Stretch is a relative value corresponding to the negative relative velocity change -dv/v. Its dimension is: :func:(len(tw),mat.shape[1])
sim_mat: 3d ndarray containing the similarity matricies that: indicate the correlation coefficient with the reference for the different time windows, different times and different amount of stretching. Its dimension is: (len(tw),mat.shape[1],len(strvec))()
second_axis: It contains the stretch vector used for the velocity: change estimate.
vale_type: It is equal to ‘shift’ and specify the content of: the returned ‘value’.
method: It is equal to ‘single_ref’ and specify in which “way” the: values have been obtained.

seismic.monitor.post_corr_process.apply_shift(data: ndarray, stats: CorrStats, shifts: ndarray) → Tuple[ndarray, CorrStats][source]#

Apply a shift to a correlation matrix.

This function applies a shift given in seconds to each trace in the matrix.

Parameters:

data (ndarray) – correlation matrix or collection of traces. lag time runs along rows and number of traces is data.shape[0].
shifts (ndarray) – shifts in seconds for each trace in ‘data’

seismic.monitor.post_corr_process.apply_stretch(data: ndarray, stats: CorrStats, stretches: ndarray) → Tuple[ndarray, CorrStats][source]#

Apply a stretches to a correlation matrix.

This function applies a stretch given in relative units to each trace in the matrix.

Parameters:

data (ndarray) – correlation matrix or collection of traces. lag time runs along rows and number of traces is data.shape[0].
stretches (ndarray) – stretches in relative units for each trace in ‘data’

seismic.monitor.stretch_mod#

Holds functions for the velocity estimation via stretching.

copyright:
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Tuesday, 15th June 2021 03:42:14 pm Last Modified: Friday, 29th September 2023 10:59:17 am

seismic.monitor.stretch_mod.time_windows_creation(starting_list: list, t_width: List[int]) → ndarray[source]#

Time windows creation.

A matrix containing one time window for each row is created. The starting samples of each one of them are passed in the starting_list parameter. The windows length t_width can be scalar or a list of values. In the latter case both lists starting_list and t_width must have the same length.

Parameters:

starting_list (list) – List of starting points
t_width (int or list of int) – Windows length

Return type:

ndarray

Returns:

tw_mat: 2d ndarray containing the indexes of one time window for each row

seismic.monitor.stretch_mod.compute_wfc(mat: ndarray, tw: ndarray, refcorr: ndarray, sides: str, remove_nans: bool = True) → ndarray[source]#

Computes the waveform coherency (WFC) between the given reference correlation and correlationmatrix.

See Steinmann, et. al. (2021) for details

Parameters:

mat (np.ndarray) – 2D Correlation Matrix (function of corrstart and time lag)
tw (np.ndarray) – Lag time window to use for the computation
refcorr (np.ndarray) – 1D reference Correlation Trace extracted from mat.
sides (str) – Which sides to use. Can be both, right, left, or single.
remove_nans (bool, optional) – Remove nans from CorrMatrix, defaults to True

Raises:

ValueError – unknown option in sides

Returns:

A 1D array with len=N(time windows)

Return type:

np.ndarray

seismic.monitor.stretch_mod.wfc_multi_reftr(corr_data: ndarray, ref_trs: ndarray, tw: ndarray, sides: str, remove_nans: bool = True) → dict[source]#

Computes the waveform coherency (WFC) between the given reference correlation(s) and correlation matrix. If several references are given, it will loop over these

See Steinmann, et. al. (2021) for details

Parameters:

corr_data (np.ndarray) – 2D Correlation Matrix (function of corrstart and time lag)
refcorr (np.ndarray) – 1 or 2D reference Correlation Trace extracted from corr_data. If 2D, it will be interpreted as one refcorr trace per row.
tw (np.ndarray) – Lag time window to use for the computation
sides (str) – Which sides to use. Can be both, right, left, or single.
remove_nans – Remove nans from CorrMatrix, defaults to True

Returns:

A dictionary with one correlation array for each reference trace.

Return type:

dict

seismic.monitor.stretch_mod.velocity_change_estimate(mat: ndarray, tw: ndarray, strrefmat: ndarray, strvec: ndarray, sides: str = 'both', return_sim_mat: bool = False, remove_nans: bool = True) → dict[source]#

Velocity change estimate through stretching and comparison.

Velocity changes are estimated comparing each correlation function stored in the mat matrix (one for each row) with strrefmat.shape[0] stretched versions of a reference trace stored in strrefmat. The stretch amount for each row of strrefmat must be passed in the strvec vector. The best match (stretch amount and corresponding correlation value) is calculated on different time windows (each row of tw is a different one) instead of on the whole trace.

Parameters:

mat (ndarray) – 2d ndarray containing the correlation functions. One for each row.
tw (ndarray of int) – 2d ndarray of time windows to be use in the velocity change estimate.
strrefmat (ndarray) – 2D array containing stretched version of the reference matrix
strvec (ndarray or list) – Stretch amount for each row of strrefmat
sides (string) – Side of the reference matrix to be used for the velocity change estimate (‘both’ | ‘left’ | ‘right’ | ‘single’)
remove_nans (bool) – If True applay nan_to_num() to the given correlation matrix before any other operation.

Return type:

Dictionary

Returns:

dv: Dictionary with the following keys

corr: 2d ndarray containing the correlation value for the best: match for each row of mat and for each time window. Its dimension is: :func:(len(tw),mat.shape[1])
value: 2d ndarray containing the stretch amount corresponding to: the best match for each row of mat and for each time window. Its dimension is: :func:(len(tw),mat.shape[1])
sim_mat: 3d ndarray containing the similarity matricies that: indicate the correlation coefficient with the reference for the different time windows, different times and different amount of stretching. Its dimension is: (len(tw),mat.shape[1],len(strvec))()
second_axis: It contains the stretch vector used for the velocity: change estimate.
vale_type: It is equal to ‘stretch’ and specify the content of: the returned ‘value’.
method: It is equal to ‘single_ref’ and specify in which “way” the: values have been obtained.

seismic.monitor.stretch_mod.time_stretch_estimate(corr_data: ndarray, ref_trc: ndarray = None, tw: ndarray = None, stretch_range: float = 0.1, stretch_steps: int = 100, sides: str = 'both', remove_nans: bool = True, ref_tr_trim: Optional[Tuple[float, float]] = None, ref_tr_stats=None) → dict[source]#

Time stretch estimate through stretch and comparison.

This function estimates stretching of the time axis of traces as it can occur if the propagation velocity changes.

Time stretching is estimated comparing each correlation function stored in the corr_data matrix (one for each row) with stretch_steps stretched versions of reference trace stored in ref_trc. The maximum amount of stretching may be passed in stretch_range. The time axis is multiplied by 1/(1 + stretch). The best match (stretching amount and corresponding correlation value) is calculated on different time windows. If tw = None the stretching is estimated on the whole trace.

Parameters:

corr_data (ndarray) – 2d ndarray containing the correlation functions. One for each row.
ref_trc (ndarray) – 1D array containing the reference trace to be shifted and compared to the individual traces in mat
tw (list of ndarray of int) – list of 1D ndarrays holding the indices of sampels in the time windows to be use in the time shift estimate. The sampels are counted from the zero lag time with the index of the first sample being 0. If tw = None the full time range is used.
stretch_range (scalar) – Maximum amount of relative stretching. Stretching and compression is tested from -stretch_range to stretch_range.
stretch_steps (scalar`) – Number of shifted version to be tested. The increment will be (2 * stretch_range) / stretch_steps
sides (str) – Side of the reference matrix to be used for the stretching estimate (‘both’ | ‘left’ | ‘right’ | ‘single’) single is used for one-sided signals from active sources with zero lag time is on the first sample. Other options assume that the zero lag time is in the center of the traces.
remove_nans (bool) – If True applay nan_to_num() to the given correlation matrix before any other operation.

Return type:

Dictionary

Returns:

dv: Dictionary with the following keys

corr: 2d ndarray containing the correlation value for the best: match for each row of mat and for each time window. Its dimension is: :func:(len(tw),mat.shape[1])
value: 2d ndarray containing the stretch amount corresponding to: the best match for each row of mat and for each time window. Stretch is a relative value corresponding to the negative relative velocity change -dv/v. Its dimension is: :func:(len(tw),mat.shape[1])
sim_mat: 3d ndarray containing the similarity matricies that: indicate the correlation coefficient with the reference for the different time windows, different times and different amount of stretching. Its dimension is: (len(tw),mat.shape[1],len(strvec))()
second_axis: It contains the stretch vector used for the velocity: change estimate.
vale_type: It is equal to ‘stretch’ and specify the content of: the returned ‘value’.
method: It is equal to ‘single_ref’ and specify in which “way” the: values have been obtained.

seismic.monitor.stretch_mod.multi_ref_vchange(corr_data: ndarray, ref_trs: ndarray, tw: ndarray = None, stretch_range: float = 0.1, stretch_steps: int = 100, sides: str = 'both', remove_nans: bool = True, ref_tr_trim: Optional[Tuple[float, float]] = None, ref_tr_stats=None) → dict[source]#

Velocity change estimate with single or multiple reference traces.

This function estimates the velocity change corresponding to each row of the corr_data matrix respect to the reference trace/s passed in ref_trs.

The velocity change is estimated comparing each correlation function stored in the corr_data matrix (one for each row) with stretch_steps stretched versions of reference/s trace stored in ref_trs. The maximum amount of stretching may be passed in stretch_range. The best match (stretching amount and corresponding correlation value) is calculated on different time windows. If tw = None the stretching is estimated on the whole trace.

The output is a dictionary with keys of the form "reftr_%d" % i: One for each reference trace. The corresponding value is aslo a dictionary that has a structure conforming with time_stretch_estimate output.

Parameters:

corr_data (ndarray) – 2d ndarray containing the correlation functions. One for each row.
ref_trc (ndarray) – 1D array containing the reference trace to be shifted and compared to the individual traces in mat
tw (list of ndarray of int) – list of 1D ndarrays holding the indices of sampels in the time windows to be use in the time shift estimate. The sampels are counted from the zero lag time with the index of the first sample being 0. If tw = None the full time range is used.
stretch_range (scalar) – Maximum amount of relative stretching. Stretching and compression is tested from -stretch_range to stretch_range.
stretch_steps (scalar`) – Number of shifted version to be tested. The increment will be (2 * stretch_range) / stretch_steps
single_sided (bool) – if True zero lag time is on the first sample. If False the zero lag time is in the center of the traces.
sides (str) – Side of the reference matrix to be used for the stretching estimate (‘both’ | ‘left’ | ‘right’ | ‘single’) single is used for one-sided signals from active sources with zero lag time is on the first sample. Other options assume that the zero lag time is in the center of the traces.
remove_nans (bool) – If True applay nan_to_num() to the given correlation matrix before any other operation.

Return type:

dictionary

Returns:

multi_ref_panel: It is a dictionary that contains as much dictionaries as reference traces have been used. The key format is "reftr_%d" % i and the corresponding value is also a dictionary with the structure described in time_stretch_estimate

seismic.monitor.stretch_mod.est_shift_from_dt_corr(dt1: ndarray, dt2: ndarray, corr1: ndarray, corr2: ndarray) → Tuple[float, float][source]#

Estimation of a baseline shift between velocity-change measurements preformed with different references.

The use of different reference traces obtaind from different reference periods will result in a shift of the velocity-change curves that ideally characterizes the velocity variation between the two reference periods. Instead of directly measuring this velocity change from the two reference traces it is calulated here as the weighted average of the point differences between the two velocity-change curves weighted by their inverse variance according to Weaver et al. GJI 2011 (On the precision of noise correlation interferometry)

Input vertors must all be of the same lenth.

Pram dt1:: Velocity variation measured for reference A
Pram dt2:: Velocity variation measured for reference B
Pram corr1:: Correlation between velocity corrected trace and reference A
Pram corr2:: Correlation between velocity corrected trace and reference B

seismic.monitor.stretch_mod.estimate_reftr_shifts_from_dt_corr(multi_ref_panel: dict, return_sim_mat: bool = False) → dict[source]#

Combine velocity-change measurements of the same data performed with different references to a single curve.

For a set of velocity-change measurements performed with different references this function estimates the relative offsets between all pairs of the measurements as a weighted average of their difference with the function est_shift_from_dt_corr. A least squares solution in computed that combines the pairwise differences to a consistent set of reference shifts. These shifts should be similar to the velocity variations measured between the reference traces. The consistent set of reference shifts is used to correct i.e. shift the similarity matricies to a common reference. Finally the corrected similarity matrices are averaged resulting in a single matrix that is interpreted as before. The position of the crest is the combined velocity change and the height of the crest is the correlation value.

Parameters:

multi_ref_panel (dictionay) – It is a dictionary with one (key,value) pair for each reference trace. Its structure is described in multi_ref_vchange()
return_sim_mat (bool) – If True the returning dictionary contains also the similarity matrix sim_mat.

Return type:

dict

Returns:

dv: Dictionary with the following keys - corr: 2d ndarray containing the correlation value for the best match for each row of mat and for each time window. Its shape is: (len(tw),mat.shape[1]) - value: 2d ndarray containing the stretch amount corresponding to the best match for each row of mat and for each time window. Stretch is a relative value corresponding to the negative relative velocity change -dv/v. Its shape is: (len(tw),mat.shape[1]) - sim_mat: 3d array containing the similarity matricies that indicate the correlation coefficient with the reference for the different time windows, different times and different amount of stretching. Its shape is: (len(tw),mat.shape[1],len(strvec)) - second_axis: It contains the stretch vector used for the velocity change estimate. - value_type: It is equal to stretch and specify the content of the returned value. - method: It is equal to single_ref and specify in which way the values have been obtained.

seismic.monitor.stretch_mod.multi_ref_vchange_and_align(corr_data: ndarray, ref_trs: ndarray, tw: ndarray = None, stretch_range: float = 0.1, stretch_steps: int = 100, sides: str = 'both', return_sim_mat: bool = False, remove_nans: bool = True, ref_tr_trim: Optional[Tuple[float, float]] = None, ref_tr_stats=None) → dict[source]#

Multi-reference dv estimate and alignment

Parameters:

corr_data (ndarray) – 2d ndarray containing the correlation functions. One for each row.
ref_trc (ndarray) – 1D array containing the reference trace to be shifted and compared to the individual traces in mat
tw (list of ndarray of int) – list of 1D ndarrays holding the indices of sampels in the time windows to be use in the time shift estimate. The sampels are counted from the zero lag time with the index of the first sample being 0. If tw = None the full time range is used.
stretch_range (scalar) – Maximum amount of relative stretching. Stretching and compression is tested from -stretch_range to stretch_range.
stretch_steps (scalar`) – Number of shifted version to be tested. The increment will be (2 * stretch_range) / stretch_steps
single_sided (bool) – if True zero lag time is on the first sample. If False the zero lag time is in the center of the traces.
sides (str) – Side of the reference matrix to be used for the stretching estimate (‘both’ | ‘left’ | ‘right’ | ‘single’) single is used for one-sided signals from active sources with zero lag time is on the first sample. Other options assume that the zero lag time is in the center of the traces.
return_sim_mat (bool) – If True the returning dictionary contains also the similarity matrix sim_mat.
remove_nans (bool) – If True applay nan_to_num() to the given correlation matrix before any other operation.

Return type:

Dictionary

Returns:

dv: Dictionary with the following keys

corr: 2d ndarray containing the correlation value for the best: match for each row of mat and for each time window. Its dimension is: :func:(len(tw),mat.shape[1])
value: 2d ndarray containing the stretch amount corresponding to: the best match for each row of mat and for each time window. Stretch is a relative value corresponding to the negative relative velocity change -dv/v. Its dimension is: :func:(len(tw),mat.shape[1])
sim_mat: 3d ndarray containing the similarity matricies that: indicate the correlation coefficient with the reference for the different time windows, different times and different amount of stretching. Its dimension is: (len(tw),mat.shape[1],len(strvec))()
second_axis: It contains the stretch vector used for the velocity: change estimate.
vale_type: It is equal to ‘stretch’ and specify the content of: the returned ‘value’.
method: It is equal to ‘single_ref’ and specify in which “way” the: values have been obtained.

seismic.monitor.stretch_mod.time_shift_estimate(corr_data: ndarray, ref_trc: ndarray = None, tw: List[ndarray] = None, shift_range: float = 10, shift_steps: int = 100, single_sided: bool = False, return_sim_mat: bool = True, remove_nans: bool = True) → dict[source]#

Time shift estimate through shifting and comparison.

This function is intended to estimate shift of traces as they can occur in noise cross-correlation in case of drifting clocks.

Parameters:

corr_data (ndarray) – 2d ndarray containing the correlation functions. One for each row.
ref_trc (ndarray) – 1D array containing the reference trace to be shifted and compared to the individual traces in mat
tw (list of ndarray of int) – list of 1D ndarrays holding the indices of sampels in the time windows to be use in the time shift estimate. The sampels are counted from the zero lag time with the index of the first sample being 0. If tw = None the full time range is used.
shift_range (scalar) – Maximum amount of time shift in samples (in one direction). Shifting is tested in both directions from -shift_range to shift_range
shift_steps (scalar`) – Number of shifted version to be tested. The increment will be (2 * shift_range) / shift_steps
single_sided – If True the zero lag time of the traces is in the first sample. If False zero lag is assumed to be in the center of the traces and the shifting is evaluated on the causal and acausal parts of the traces separately and averaged. This is done to avoid bias from velocity changes (stretching) in the case of strongly asymmetric traces.
remove_nans (bool) – If True applay nan_to_num() to the given correlation matrix before any other operation.

Return type:

Dictionary

Returns:

dt: Dictionary with the following keys

corr: ndarray 2d ndarray containing the correlation: value for the best match for each row of mat and for each time window. Its dimension is: :func:(len(tw),mat.shape[1])
value: ndarray 2d ndarray containing the amount of: shifting corresponding to the best match for each row of mat and for each time window. Shift is measured in units of the sampling interval. Its dimension is: (len(tw),mat.shape[1])()
sim_mat: 2d ndarray containing the similarity matrix that: indicate the time shift respect to the reference for the selected time windows, for different times and different amount of shift. Its dimension is: (mat.shape[1],shift_steps)()
second_axis: It contains the shift vector used for the velocity: change estimate.
vale_type: It is equal to ‘shift’ and specify the content of: the returned ‘value’.
method: It is equal to ‘time_shift’ and specify in which “way” the: values have been obtained.

seismic.monitor.stretch_mod.time_shift_apply(corr_data: ndarray, shift: ndarray, single_sided: bool = False) → ndarray[source]#

Correct for clock drifts that are given in shifts.

Parameters:

corr_data (np.ndarray) – Matrix with one correlation trace per row
shift (np.ndarray) – 1 or 2D array holding shifts in number of samples compared to a reference that is assumed to be accurate. Hence, the input CorrStream will be shifted by -shift. If 1D: each correlation trace will be shifted by a scalar value. If 2D: Shape has to be equal to corr_data. The clock drift varies in the CorrTrace (i..e, shift is not constant anymore).
single_sided (bool, optional) – corr_data contains only causal side, defaults to False

Returns:

The shifted correlation matrix

Return type:

np.ndarray

seismic.monitor.stretch_mod.time_stretch_apply(corr_data: ndarray, stretch: ndarray, single_sided: bool = False) → ndarray[source]#

Apply time axis stretch to traces.

Stretch the time axis of traces e.g. to compensate a velocity shift in the propagation medium. Such shifts can occur in corrlation traces in case of a drifting clock. This function applies the stretches. To correct for stretching estimated with time_stretch_estimate() you need to apply negative stretching.

Parameters:

corr_data (ndarray) – 2d ndarray containing the correlation functions that are to be shifted. One for each row.
stretch (ndarray) – ndarray with stretch.shape[0] = corr_data.shape[0] containing the stretches relative units.

Return type:

ndarray

Returns:

stretched_mat: stretched version of the input matrix

seismic.monitor.stretch_mod.create_shifted_ref_mat(ref_trc: ndarray, stats: CorrStats, shifts: ndarray) → array[source]#: Create a matrix of shifted versions of a reference trace.

seismic.monitor.stretch_mod.compare_with_modified_reference(data: ndarray, ref_mat: ndarray, indices: ndarray) → ndarray[source]#: Compare a correlation matrix with a modified references.

seismic.monitor.trim#

Trim correlation matrices.

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Tuesday, 15th November 2022 05:15:30 pm Last Modified: Monday, 16th January 2023 11:14:37 am

seismic.monitor.trim.corr_mat_trim(data: ndarray, stats: CorrStats, starttime: float, endtime: float) → Tuple[ndarray, CorrStats][source]#

Trim the correlation matrix to a given period.

Trim the correlation matrix corr_mat to the period from starttime to endtime given in seconds from the zero position, so both can be positive and negative. If starttime and endtime are datetime.datetime objects they are taken as absolute start and endtime.

Parameters:

corr_mat – correlation matrix to be trimmed
starttime (float) – start time in seconds with respect to the zero position. Hence, this parameter describes a lag!
order – end time in seconds with respect to the zero position

Rtype tuple:

Tuple[np.ndarray, CorrStats]

Returns:

trimmed correlation matrix and new stats object

seismic.plot#

seismic.plot.plot_correlation#

Plot Green’s function estimations.

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: GNU Lesser General Public License, Version 3 (https://www.gnu.org/copyleft/lesser.html)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Monday, 19th July 2021 11:37:54 am Last Modified: Monday, 4th March 2024 02:42:42 pm

seismic.plot.plot_correlation.plot_ctr(corr: Trace, tlim: list = None, ax: Axes = None, outputdir: str = None, clean: bool = False) → Axes[source]#

seismic.plot.plot_correlation.plot_cst(cst: Stream, sort_by: str = 'corr_start', timelimits: list = None, ylimits: list = None, scalingfactor: float = 2.0, ax: Axes = None, linewidth: float = 0.25, outputfile: str = None, title: str = None, type: str = 'heatmap', cmap: str = 'seismic', vmin: float = None, vmax: float = None, **kwargs)[source]#

Creates a section or heat plot of all correlations in this stream.

Parameters:

cst (CorrStream) – Input CorrelationStream
sort_by (str, optional) – Which parameter to plot against. Can be either corr_start or distance, defaults to ‘corr_start’.
timelimits (Tuple[float, float], optional) – xlimits (lag time) in seconds, defaults to None
ylimits (Tuple[float, float], optional) – limits for Y-axis (either a datetime.datetime or float in km (if plotted against distance)), defaults to None.
scalingfactor (float, optional) – Which factor to scale the Correlations with. Play around with this if you want to make amplitudes bigger or smaller, defaults to None (automatically chosen).
ax (plt.Axes, optional) – Plot in existing axes? Defaults to None
linewidth (float, optional) – Width of the lines to plot, defaults to 0.25
outputfile (str, optional) – Save the plot? defaults to None
title (str, optional) – Title of the plot, defaults to None
type (str, optional) – can be set to ‘heat’ for a heatmap or ‘section’ for a wiggle type plot. Default to heat
cmap (str, optional) – Decides about colormap if type == ‘heatmap’. Defaults to ‘inferno’.

Returns:

returns the axes object.

Return type:

matplotlib.pyplot.Axes

seismic.plot.plot_correlation.plot_corr_bulk(corr_bulk, timelimits: Optional[Tuple[datetime, datetime]] = None, ylimits: Optional[Tuple[datetime, datetime]] = None, clim: Optional[Tuple[float, float]] = None, plot_colorbar: bool = False, outputfile: Optional[str] = None, title: Optional[str] = None, ax: Optional[Axes] = None) → Axes[source]#

Plots a CorrBulk object.

Parameters:

corr_bulk (CorrBulk) – The CorrBulk to plot.
timelimits (Optional[Tuple[datetime.datetime, datetime.datetime]], optional) – Limits time axis, defaults to None
ylimits (Optional[Tuple[datetime.datetime, datetime.datetime]], optional) – Limits of y-axis, defaults to None
clim (Optional[Tuple[float, float]], optional) – Limits of Colobar, defaults to None
plot_colorbar (bool, optional) – add colorbar to plot, defaults to False
outputfile (Optional[str], optional) – save file to, defaults to None
title (Optional[str], optional) – Title, defaults to None
ax (Optional[plt.Axes], optional) – Axis to plot into, defaults to None

Returns:

The current axis

Return type:

plt.Axes

seismic.plot.plot_correlation.sect_plot_corr_start(cst: Stream, ax: Axes, scalingfactor: float, linewidth: float) → ndarray[source]#

seismic.plot.plot_correlation.heat_plot_corr_start(cst: Stream, ax: Axes, cmap: str, vmin: float, vmax: float)[source]#

seismic.plot.plot_correlation.sect_plot_dist(cst: Stream, ax: Axes, scalingfactor: float, linewidth: float, plot_reference_v: bool = False, ref_v: List[float] = [1, 2, 3]) → ndarray[source]#

seismic.plot.plot_correlation.heat_plot_corr_dist(cst: Stream, ax: Axes, cmap: str, vmin: float, vmax: float, plot_reference_v: bool = False, ref_v: List[float] = [1, 2, 3])[source]#

seismic.plot.plot_dv#

Plot velocity change time series.

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: GNU Lesser General Public License, Version 3 (https://www.gnu.org/copyleft/lesser.html)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Friday, 16th July 2021 02:30:02 pm Last Modified: Monday, 4th March 2024 02:44:41 pm

seismic.plot.plot_dv.plot_fancy_dv(dv, xlim: Tuple[datetime, datetime] = None, ylim: Tuple[float, float] = None, return_ax: bool = False, title: str = None, dateformat: str = '%d %b %y', ax: axes = None, *args, **kwargs)[source]#: Prettier than the technical plot, but less informative

seismic.plot.plot_dv.plot_technical_dv(dv, save_dir='.', figure_file_name=None, mark_time=None, normalize_simmat=False, sim_mat_Clim=[], figsize=(9, 11), dpi=72, ylim: Tuple[float, float] = None, xlim: Tuple[datetime, datetime] = None, title: str = None, plot_scatter: bool = False, return_ax: bool = False, *args, **kwargs) → Tuple[figure, List[axis]][source]#

Plot a seismic.monitor.dv.DV object The produced figure is saved in save_dir that, if necessary, it is created. It is also possible to pass a “special” time value mark_time that will be represented in the dv/v and corr plot as a vertical line; It can be a string (i.e. YYYY-MM-DD) or directly a datetime object. if the dv dictionary also contains a ‘comb_mseedid’ keyword, its value (i.e. MUST be a string) will be reported in the title. In case of the chosen filename exist in the save_dir directory, a prefix _<n> with n:0..+Inf, is added. The aspect of the plot may change depending on the matplotlib version. The recommended one is matplotlib 1.1.1

Parameters:

dv (dict) – velocity change estimate dictionary as output by multi_ref_vchange_and_align and successively “extended” to contain also the timing in the form {‘time’: time_vect} where time_vect is a ndarray of datetime objects.
save_dir (string) – Directory where to save the produced figure. It is created if it doesn’t exist.
figure_file_name (string) – filename to use for saving the figure. If None the figure is displayed in interactive mode.
mark_time (string or datetime object) – It is a “special” time location that will be represented in the dv/v and corr plot as a vertical line.
normalize_simmat (Bool) – if True the simmat will be normalized to a maximum of 1. Defaults to False
sim_mat_Clim (2 element array_like) – if non-empty it set the color scale limits of the similarity matrix image
ylim (Tuple[float, float], optional) – Limits for the stretch axis. Defaults to None
return_ax (bool, optional) – Return plt.figure and list of axes. Defaults to False. This overwrites any choice to save the figure.

Returns:

If return_ax is set to True it returns fig and axes.

seismic.plot.plot_dv.plot_dv(style: str, *args, **kwargs)[source]#

seismic.plot.plot_multidv#

Plot velocity change time series.

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: GNU Lesser General Public License, Version 3 (https://www.gnu.org/copyleft/lesser.html)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Tuesday, 5th October 2021 11:50:22 am Last Modified: Thursday, 21st October 2021 02:44:05 pm

seismic.plot.plot_multidv.plot_multiple_dv(indir: str, only_mean: bool = False, title: str = None, outfile: str = None, fmt: str = None, dpi: int = 300, legend: bool = False, plot_median: bool = False, ylim: Tuple[float, float] = None, xlim: Tuple[datetime, datetime] = None, statfilter: List[str] = None, plt_abs: bool = False)[source]#

Plots several Velocity variations in one single plot

Parameters:

indir (str) – directory in which the .npz DV files are located.
only_mean (bool, optional) – plot only the average of one station instead of all of its components. The averages have to be computed before. Defaults to false.
title (str, optional) – Title for the figure, defaults to None
outfile (str, optional) – File to save figure to, defaults to None
fmt (str, optional) – Figure format, defaults to None
dpi (int, optional) – DPI for pixel-based formats, defaults to 300
legend (bool, optional) – Should a legend be plotted, defaults to False
plot_median (bool, optional.) – Plot the median of all dataset as a black line. Defaults to False.
ylim (Tuple[float, float], optional) – y limit, defaults to None
xlim (Tuple[datetime, datetime], optional) – x limit, defaults to None
statfilter (List[str]) – Only plot data from the station combinations with the following codes. Given in the form [‘net0-net0.sta0-stat1.ch0-ch1’, …]. Note that wildcards are allowed. Defaults to None.
plt_abs (bool, optional) – plot absolute values. Defaults to False

seismic.plot.plot_utils#

Useful little plotting functions

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Monday, 17th May 2021 12:25:54 pm Last Modified: Monday, 4th March 2024 02:43:54 pm

seismic.plot.plot_utils.set_mpl_params()[source]#

seismic.plot.plot_utils.remove_all(ax: Axes = None, top=False, bottom=False, left=False, right=False, xticks='none', yticks='none')[source]#: Removes all frames and ticks.

seismic.plot.plot_utils.remove_topright(ax: Axes = None)[source]#: Removes top and right border and ticks from input axes.

seismic.plot.plot_spectrum#

Plot spectrograms.

Plotting specral time series.

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE Version 1.2 (https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Wednesday, 21st June 2023 04:54:20 pm Last Modified: Wednesday, 7th August 2024 02:39:58 pm

seismic.plot.plot_spectrum.plot_spct_series(S: ndarray, f: ndarray, t: ndarray, norm: str = None, norm_method: str = None, flim: Tuple[int, int] = None, tlim: Tuple[datetime, datetime] = None, cmap: str = 'batlow', vmin=None, vmax=None, log_scale: bool = False) → Axes[source]#

Plots a spectral series.

Parameters:

S (np.ndarray) – A spectral series with dim=2.
f (np.ndarray) – Vector containing frequencies (Hz)
t (np.ndarray) – Vector containing times as UTCDateTimes
norm (str, optional) – Normalise the spectrum either on the time axis with norm=’t’ or on the frequency axis with norm=’f’, defaults to None.
norm_method – Normation method to use. Either ‘linalg’ (i.e., length of vector), ‘mean’, or ‘median’.
title (str, optional) – Plot title, defaults to None
outfile (str, optional) – location to save the plot to, defaults to None
fmt (str, optional) – Format to save figure, defaults to ‘pdf’
flim (Tuple[int, int]) – Limit Frequency axis and Normalisation to the values in the given window
tlim (Tuple[datetime, datetime]) – Limit time axis to the values in the given window

seismic.plot.plot_wfc#

Plot the waveform coherence

copyright:: Module to plot waveform coherencies.

The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Monday, 8th November 2021 02:46:15 pm Last Modified: Tuesday, 11th July 2023 01:46:20 pm

seismic.plot.plot_wfc.plot_wfc_bulk(time_lag: ndarray, cfreq: ndarray, wfc: ndarray, title: str, log: bool, cmap: str)[source]#

seismic.utils#

Collection of tools

seismic.utils.io#

Load correlations produced by MIIC

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: GNU Lesser General Public License, Version 3 (https://www.gnu.org/copyleft/lesser.html)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Thursday, 8th July 2021 10:37:21 am Last Modified: Friday, 10th February 2023 04:09:21 pm

seismic.utils.io.flatten(x)[source]#

Return the flattened version of the input array x

This funciton works with all Iterable that can be nested in an irregular fashion.

Param:: Iterable to be flattened
Return type:: list
Returns:: x as a flattened list

seismic.utils.io.flatten_recarray(x)[source]#: Flatten a recarray

seismic.utils.io.mat_to_corrtrace(mat: dict) → CorrTrace[source]#

Reads in a correlation dictionary created with the old (i.e., python2) miic and translates it into a CorrTrace object.

Parameters:: mat (dict) – dictionary as produced by scipy.io.loadmat()
Returns:: A Correlation Trace
Return type:: CorrTrace

seismic.utils.io.corrmat_to_corrbulk(path: str) → CorrBulk[source]#

Reads in a correlation .mat file created with the old (i.e., python2) miic and translates it into a CorrTrace object.

Parameters:: path (str) – Path to file
Returns:: A Correlation Bulk object
Return type:: CorrBulk
Warning:: This reads both stacks and subdivisions into the same object.

seismic.utils.fetch_func_from_str#

Import and Return a function

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Friday, 9th April 2021 04:07:06 pm Last Modified: Monday, 16th January 2023 11:13:58 am

seismic.utils.fetch_func_from_str.func_from_str(string: str) → Callable[source]#

Given a full string in the form ‘module.submodule.function’. This function imports said function and returns it, so that it can be called. The actual function can be located in any module.

Parameters:: string (str) – string in the form module.submodule.function
Returns:: The imported function ready to be executed
Return type:: function

seismic.utils.miic_utils#

Set of smaller functions.

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE v. 1.2

(https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Monday, 29th March 2021 12:54:05 pm Last Modified: Wednesday, 19th June 2024 03:13:18 pm

seismic.utils.miic_utils.trace_calc_az_baz_dist(stats1: Stats, stats2: Stats) → Tuple[float, float, float][source]#

Return azimuth, back azimhut and distance between tr1 and tr2 This funtions calculates the azimut, back azimut and distance between tr1 and tr2 if both have geo information in their stats dictonary. Required fields are: tr.stats.sac.stla and tr.stats.sac.stlo

Parameters:

stats1 (Stats) – First trace to account
stats2 (Stats) – Second trace to account

Return type:

float

Returns:

az: Azimuth angle between tr1 and tr2

Return type:

float

Returns:

baz: Back-azimuth angle between tr1 and tr2

Return type:

float

Returns:

dist: Distance between tr1 and tr2 in m

seismic.utils.miic_utils.filter_stat_dist(inv1: Inventory, inv2: Inventory, thres: float) → bool[source]#

Very simple function to check whether to stations are closer than thres to each other.

Parameters:

inv1 (Inventory) – Inventory of station 1
inv2 (Inventory) – Inventory of station 2
thres (float) – Threshold distance in m

Returns:

True if closer (or equal) than thres, False if not.

Return type:

bool

seismic.utils.miic_utils.inv_calc_az_baz_dist(inv1: Inventory, inv2: Inventory) → Tuple[float, float, float][source]#

Return azimuth, back azimuth and distance between stat1 and stat2

Parameters:

tr1 (Inventory) – First trace to account
tr2 (Inventory) – Second trace to account

Return type:

float

Returns:

az: Azimuth angle between stat1 and stat2

Return type:

float

Returns:

baz: Back-azimuth angle between stat2 and stat2

Return type:

float

Returns:

dist: Distance between stat1 and stat2 in m

seismic.utils.miic_utils.resample_or_decimate(data: Trace, sampling_rate_new: int, filter=True) → Stream[source]#

Decimates the data if the desired new sampling rate allows to do so. Else the signal will be interpolated (a lot slower).

Parameters:

data (Stream) – Stream to be resampled.
sampling_rate_new (int) – The desired new sampling rate

Returns:

The resampled stream

Return type:

Stream

seismic.utils.miic_utils.trim_stream_delta(st: Stream, start: float, end: float, *args, **kwargs) → Stream[source]#

Cut all traces to starttime+start and endtime-end. args and kwargs will be passed to trim()

Parameters:

st (Stream) – Input Stream
start (float) – Delta to add to old starttime (in seconds)
end (float) – Delta to remove from old endtime (in seconds)

Returns:

The trimmed stream

Return type:

Stream

Note

This operation is performed in place on the actual data arrays. The raw data will no longer be accessible afterwards. To keep your original data, use copy() to create a copy of your stream object.

seismic.utils.miic_utils.trim_trace_delta(tr: Trace, start: float, end: float, *args, **kwargs) → Trace[source]#

Cut all traces to starttime+start and endtime-end. args and kwargs will be passed to trim().

Parameters:

st (Trace) – Input Trace
start (float) – Delta to add to old starttime (in seconds)
end (float) – Delta to remove from old endtime (in seconds)

Returns:

The trimmed trace

Return type:

Trace

Note

This operation is performed in place on the actual data arrays. The raw data will no longer be accessible afterwards. To keep your original data, use copy() to create a copy of your stream object.

seismic.utils.miic_utils.save_header_to_np_array(stats: Stats) → dict[source]#

Converts an obspy header to a format that allows it to be saved in an npz file (i.e., several gzipped npy files)

Parameters:: stats (Stats) – input Header
Returns:: Dictionary, whose keys are the names of the arrays and the values the arrays themselves. Can be fed into np.savez as **kwargs
Return type:: dict

seismic.utils.miic_utils.load_header_from_np_array(array_dict: dict) → dict[source]#

Takes the dictionary-like return-value of np.load and converts the corresponding keywords into an obspy header.

Parameters:: array_dict (dict) – Return value of np.load
Returns:: The obspy header as dictionary
Return type:: dict

seismic.utils.miic_utils.convert_utc_to_timestamp(utcdt: UTCDateTime) → ndarray[source]#

Converts obspy.core.utcdatetime.UTCDateTime objects to floats.

Parameters:: utcdt (UTCDateTimeorList[UTCDateTime]) – The input times, either a list of utcdatetimes or one utcdatetime
Returns:: A numpy array of timestamps
Return type:: np.ndarray

seismic.utils.miic_utils.convert_timestamp_to_utcdt(timestamp: ndarray) → List[UTCDateTime][source]#

Converts a numpy array holding timestamps (i.e., floats) to a list of UTCDateTime objects

Parameters:: timestamp (np.ndarray) – numpy array holding timestamps
Returns:: a list of UTCDateTime objects
Return type:: List[UTCDateTime]

seismic.utils.miic_utils.get_valid_traces(st: Stream)[source]#

Return only valid traces of a stream.

Remove traces that are 100% masked from a stream. This happens when a masked trace is trimmed within a gap. The function works in place.

Parameters:: st (obspy.Stream) – stream to work on

seismic.utils.miic_utils.discard_short_traces(st: Stream, length: float)[source]#

Discard all traces from stream that are shorter than length.

Parameters:

st (Stream) – inputer obspy Stream
length (float) – Maxixmum Length that should be discarded (in seconds).

Note

Action is performed in place.

seismic.utils.miic_utils.nan_moving_av(data: ndarray, win_half_len: int, axis: int = -1) → ndarray[source]#

Returns a filtered version of data, disregarding the nans. Moving mean window length is win_half_len*2+1.

Parameters:

data (np.ndarray) – Array to be filtered
win_half_len (int) – Half length of the boxcar filter (len = halflen*2+1)
axis (int, optional) – Axis to filter along, defaults to -1

Returns:

The filtered array

Return type:

np.ndarray

seismic.utils.miic_utils.stream_require_dtype(st: Stream, dtype: type) → Stream[source]#

Often it might make sense to change the data type of seismic data before saving or broadcasting (e.g., to save memory). This function allows to do so

Parameters:

st (Stream) – input Stream
dtype (type) – desired datatype, e.g. np.float32

Returns:

Stream with data in new dtype

Return type:

Stream

Note

This operation is performed in place.

seismic.utils.miic_utils.correct_polarity(st: Stream, inv: Inventory)[source]#

Sometimes the polarity of a seismometer’s vertical component is reversed. This simple functions flips the polarity if this should be the case.

Note

This function acts on the data in-place.

Parameters:

st (Stream) – input stream. If there is no Z-component, nothing will happen.
inv (Inventory) – Inventory holding the orientation information.

seismic.utils.miic_utils.nan_helper(y: ndarray) → Tuple[ndarray, ndarray][source]#

Helper to handle indices and logical indices of NaNs.

Parameters:: y (1d np.ndarray) – numpy array with possible NaNs
Returns:: Tuple of two numpy arrays: The first one holds a logical mask indicating the positions of nans, the second one is a function that can be called to translate the first output to indices.
Return type:: Tuple[np.ndarray, np.ndarray]

seismic.utils.miic_utils.gap_handler(st: Stream, max_interpolation_length: int = 10, retain_len: float = 0, taper_len: float = 10) → Stream[source]#

This function tries to interpolate gaps in a stream. If the gaps are too large and the data snippets too short, the data will be discarded. Around remaining gaps, the trace will be tapered.

Parameters:

st (Stream) – Stream that might contain gaps
max_interpolation_length (int, optional) – maximum length in number of samples to interpolate, defaults to 10
retain_len (int, optional) – length in seconds that a data snippet has to have to be retained, defaults to 0
taper_len (float, optional) – Length of the taper in seconds, defaults to 10

Returns:

The sanitised stream

Return type:

Stream

seismic.utils.miic_utils.interpolate_gaps(A: ndarray, max_gap_len: int = -1) → ndarray[source]#

perform a linear interpolation where A is filled with nans.

Parameters:

A (np.ndarray) – np.ndarray, containing nans
max_gap_len (int, optional) – maximum length to fill, defaults to -1

Returns:

The array with interpolated values

Return type:

np.ndarray

seismic.utils.miic_utils.interpolate_gaps_st(st: Stream, max_gap_len: int = -1) → Stream[source]#

perform a linear interpolation where A is filled with nans.

Parameters:

st (Stream) – Stream containing nans
max_gap_len (int, optional) – maximum length to fill, defaults to -1

Returns:

The stream with interpolated values

Return type:

Stream

seismic.utils.miic_utils.sort_combinations_alphabetically(netcomb: str, stacomb: str, loccomb: str, chacomb: str) → Tuple[str, str, str, str][source]#

Sort the combinations of network, station, location and channel alphabetically to avoid ambiguities.

Parameters:

netcomb (str) – Network combination
stacomb (str) – Station combination
loccomb (str) – Location combination
chacomb (str) – Channel combination

Returns:

The sorted combinations

Return type:

Tuple[str, str, str, str]

seismic.utils.raw_analysis#

Analysing raw waveform data.

Module for waveform data analysis. Contains spectrogram computation.

copyright:: The SeisMIC development team (makus@gfz-potsdam.de).
license:: EUROPEAN UNION PUBLIC LICENCE Version 1.2 (https://joinup.ec.europa.eu/collection/eupl/eupl-text-eupl-12)
author:: Peter Makus (makus@gfz-potsdam.de)

Created: Wednesday, 21st June 2023 12:22:00 pm Last Modified: Thursday, 31st August 2023 04:40:30 pm

seismic.utils.raw_analysis.spct_series_welch(streams: Iterator[Stream], window_length: int, freqmax: float, remove_response: bool = True)[source]#

Computes a spectral time series. Each point in time is computed using the welch method. Windows overlap by half the windolength. The input stream can contain one or several traces from the same station. Frequency axis is logarithmic.

Note

MPI support since version 0.4.2. Then the result is only available on rank 0. All other ranks return None. (to save RAM)

Parameters:

st (Stream) – Input Stream with data from one station.
window_length (int or float) – window length in seconds for each datapoint in time
freqmax (float) – maximum frequency to be considered

Returns:

Arrays containing a frequency and time series and the spectral series.

Return type:

np.ndarray

seismic.utils.raw_analysis.preprocess(tr: Trace, freqmax: float, remove_response: bool)[source]#: Some very basic preprocessing on the string in order to plot the spectral series. Does the following steps: 1. Remove station response 2. Detrend 3. Decimate if sampling rate>50 4. highpass filter with corner period of 300s. :param st: The input Stream, should only contain Traces from one station. :type st: ~obspy.core.Stream :return: The output stream and station inventory object :rtype: ~obspy.core.Stream and ~obspy.core.Inventory