relmt.utils module

relmt.utils.approx_time_lookup(time1, time2, include_at=inf)

Create lookup table for time1 approximately matching time2

All events in time1s must be present in time2, i.e. len(time1) <= len(time2).

Parameters:

• time1 (list[str] | list[float]) – Times either as date_string understood by datetime.datetime.fromisoformat(), or as float in arbitrary units

• time2 (list[str] | list[float]) – Times either as date_string understood by datetime.datetime.fromisoformat(), or as float in arbitrary units

• include_at (float) – Maximum time difference (seconds) to consider a match.

Returns:

dict – Lookup table time1 -> time2

Raises:

LookupError: – When time1 is longer than time2

relmt.utils.cartesian_distance(x1, y1, z1, x2, y2, z2)

Cartesian distance between points (x1, y1, z1) and (x2, y2, z2)

Parameters:

• x1 (float | ndarray)

• y1 (float | ndarray)

• z1 (float | ndarray)

• x2 (float | ndarray)

• y2 (float | ndarray)

• z2 (float | ndarray)

Return type:

float | ndarray

relmt.utils.ccorf3_all(gmat, batch_size=200000)

Triplet-wise S-wave correlations of seismogram matrix Generalized and optimized ccorf3 for all column triples.

Parameters:

• gmat (ndarray) – (samples, events) Input matrix. Each event is assumed to be normalized

• batch_size (int) – Process up to this many triples per batch to control memory.

Returns:

cc_all ((C(event,3), 3) float64) – (events, events, events) S-wave cross correlation coefficient per event triplet

relmt.utils.collect_takeoff(phase_dict, event_index, stations=None)

Take-off azimuth and plunge angles for an event

Parameters:

• phase_dict (dict[slice(<class ‘str’>, <class ‘relmt.core.Phase’>, None)]) – Lookup table for phase arrivals

• event_index (int) – Index of the event to consider

• stations (set[str] | list[str] | None) – Limit output to certain stations (e.g. those with actual amplitude observations)

Returns:

• azimuths (np.ndarray) – Take-off azimuths (degree)

• plunges (np.ndarray) – Take-off plunges (degree)

• stas (np.ndarray) – Corrsponding station names

• phases (np.ndarray) – Corresponding phase types

relmt.utils.common_event_phase_lowpass(passbands, quantile)

Find common lowpass frequency across events for each station and phase

Parameters:

• passbands (dict[str, dict[int, tuple[float, float]]]) – Mapping from waveform ID to mapping of event number to high- and lowpass filter corener (Hz)

• quantile (float) – Quantile of the lowpass frequencies across events to use as the common lowpass frequency

Returns:

dict[str, dict[int, tuple[float, float]]] – New passbands dictionary with a common lowpass frequency for each phase of each event.

relmt.utils.concat_components(arr)

Rearange waveform array so that components are concatenated

Parameters:

arr (ndarray) – Waveform array of shape (events, components, samples)

Returns:

ndarray – (events, components * samples) Waveform matrix

relmt.utils.corner_frequency(magnitude, phase, stress_drop, vs)

Return a rought estimate of the corner frequency based on (Madariaga, 1976):

\[f_c = (16/7 \sigma k_s^3 v_S^3 / M_0 )^{1/3}\]

Parameters:

• magnitude (float) – Magnitude \(M_0\)

• phase (str) – P (\(k_s\) = 0.32) or S (\(k_s\) = 0.21)

• stress_drop (float) – Stress drop \(\sigma\) in Pa (use ~5e6)

• vs (float) – S wave velocity \(v_S\) in m/s

Returns:

float – Corner frequency in Hz

relmt.utils.element_redundancy(pairs, ignore=None)

Count of the number of contributing elements per pair

Parameters:

• pairs (ndarray) – (n, 2) array of event index pairs

• ignore (list[int] | None) – Do not count elements that contain these event indices

Returns:

ndarray – Number of element occurrences contributing to each pair

relmt.utils.event_indices(amps)

Index arrays of unique event numbers

For each event, find the constituting pair or triplet wise amplitude observations

Parameters:

amps (list[P_Amplitude_Ratio] | list[S_Amplitude_Ratios]) – List of amplitude observations, either P or S

Returns:

dict[int, ndarray] – Mapping from event number to observation indices

relmt.utils.fftw_data_window(minimum_time, sampling_rate, components)

Set ‘data_window’ based on sampling rate and number of components

Parameters:

• minimum_time (float) – Minimum length of single-channel time window

• sampling_rate (float) – Sampling rate of the seismic waveform (Hertz)

• components (str) – One-character component names ordered as in the waveform array, as one string (e.g. ‘ZNE’)

Returns:

float – Length of single-channel time window. The concated components will be fftw optimized

relmt.utils.fisher_average(ccarr, axis=-1)

Average cross correlation coefficients applying Fisher z-transform

Parameters:

• ccarr (ndarray) – 1, 2, or 3 dimensional array holding values bound in interval [-1, 1]

• axis (int) – Array axis along which to operate

Returns:

ndarray – Averaged array, reduced by one dimension

relmt.utils.interpolate_phase_dict(phase_dict, event_dict, station_dict, obs_min=1)

Interpolate phase arrivals, ray azimuth and plunge

Arrival time is interpolated via an average velocity along the ray path.

Ray plunge is interpolated linearly between neighboring points in depth

Azimuth is calculated assuming a horizontally straight ray.

Parameters:

• phase_dict (dict[str, Phase]) – All seismic phases

• event_dict (dict[int, Event]) – The seismic event catalog

• station_dict (dict[str, Station]) – Station table

• obs_min (int) – Minimum number of observations to make interpolation

Returns:

dict[str, Phase] – New phase dictionary containing the interpolated phases

relmt.utils.item_count(array)

Count of the items in the list

Parameters:

array (Iterable)

Return type:

ndarray

relmt.utils.mt_clusters(mt_dict, method='kagan', distance_matrix=None, max_distance=30, min_ev=1, link_method='average')

Cluster moment tensors

Parameters:

• mt_dict (dict[int, MT]) – list of moment tensors to cluster

• method (str) – Method to compare two events. Choices are: - ‘ccp’ Correlation coefficient of the P radiation - ‘ccs’ Correlation coefficient of the S radiation - ‘scalar’ Normaized scalar product of the tenors - ‘kagan’ Kagan angle of the DC components

• distance_matrix (ndarray | None) – Ignore ‘method’, but give a pre-computed distance matrix. Can Accalerate the process in finding max_dist and ev_min. Second return argument from a previous run is suited. See scipy.spacial.distance.pdist() for format.

• max_distance (float) – maximum distance for pairs to include in a cluster

• min_ev (int) – Minimum number of events per cluster. Unclustered events will receive 0 label.

• link_method (str) – Method of linking events passed on to scipy.cluster.hirearchy.linkage()

Returns:

• labels – Labels of the event clusters, where 0 is the label for unclustered events

• distance_matrix – Compressed pairwise distance matrix

• representative – The representative element for each label

relmt.utils.nearest_neighbors(n_neighbors, events, event_dict)

Parameters:

• n_neighbors (int)

• events (ndarray)

• event_dict (dict)

Return type:

set[tuple[int, int]]

relmt.utils.next_fftw_size(samples, devisor=3, odd=True)

The next integer for which FFTW is optimized

Parameters:

• samlpes – Number of target samples

• devisor (int) – The result must be devisable by this number

• odd (bool) – The result divided by devisor must be odd

• samples (int)

Return type:

int

relmt.utils.pair_redundancy(triplets, ignore=None)

Count of the number of contributing pairs per triplet

Parameters:

• triplets (ndarray) – (n, 3) array of event index triplets

• ignore (list[int] | None) – Do not count pairs that contain these event indices

Returns:

ndarray – Number of pairs contributing to each triplet

relmt.utils.pc_index(mtx, phase)

Return principal component sorting of seismogram matrix

For P-phases, sort by zero-th left hand singular vector For S-phases, sort angle in the plane spanned by the zero-th and first left hand singular vectors

Parameters:

• mtx (ndarray) – Waveform matrix of shape (events, components * samples)

• phase (str) – P or S phase

Returns:

ndarray – Index array

Raises:

ValueError: – If phase is not P or S

relmt.utils.phase_dict_azimuth(phase_dict, event_dict, station_dict, overwrite=False)

Fill phase dictionary with azimuth using trigonometry

Parameters:

• phase_dict (dict[str, Phase]) – All seismic phases

• event_dict (dict[int, Event]) – The seismic event catalog

• station_dict (dict[str, Station]) – Station table

• overwrite (bool) – Overwrite existing azimuth values (False: only replace NaN values)

Returns:

dict[str, Phase] – New phase dictionary containing computed azimuths

relmt.utils.phase_dict_hash_plunge(phase_dict, event_dict, station_dict, vmodel, overwrite=False, strict=False, nquerry=100, nray=1000, station_depth=True)

Fill phase dictionary with plunge values from HASH.

Uses the python implementation of the HASH ray-tracer (Skoumal, Hardebeck & Shearer, 2024; Hardebeck & Shearer, 2002) by courtesey of USGS.

Parameters:

• phase_dict (dict[str, Phase]) – All seismic phases

• event_dict (dict[int, Event]) – The seismic event catalog

• station_dict (dict[str, Station]) – Station table

• vmodel (ndarray) – (layers, 3) array of depth (m), P- and S-wave velocities (m/s)

• overwrite (bool) – Overwrite existing plunge values (False: only replace NaN values)

• strict (bool) – Raise KeyError when an event or station is missing. If False, check phase_dict for missing events and stations.

• nquerry (int) – Number of distance and depth querry points for plunge lookup table

• nray – Number of trail rays (should be larger than nquerry)

• station_depth (bool) – If True, consider station depth by cropping velocity model above to the depth of the station. Produces NaN plunges for events above the station. False places all stations at the shallowest depth of the velocity model, consistent with SKHASH behavior.

Returns:

dict[str, Phase] – New phase dictionary containing computed plunges

relmt.utils.reshape_ccvec(ccvec, ns, combinations=None)

Reshape 1D vector of cc coefficeints c_ijk to (n, n, n) 3D matrix

Parameters:

• ccvec (ndarray)

• ns (int)

Return type:

ndarray

relmt.utils.select_events(arr, select, events)

Return waveforms of specific evetn numbers

arr:

(events, ...) or waveform array or matrix

select:

List of event numbers to select

events:

List of event numbers in the array

Returns:

ndarray – (selected_events, ...) Waveform array or matrix

Parameters:

• arr (ndarray)

• select (list[int])

• events (list[int])

relmt.utils.signed_log(x)

Signed log transform: sign(x) log(1 + | x |)

relmt.utils.signed_log_inverse(y)

Inverse of the signed log transform: sign(x) log(1 + | x |)

relmt.utils.source_duration(magnitude)

Approximate duration of an earthquake with given magnitude

\[t = 3^{M-5}\]

Parameters:

magnitude (float) – Magnitude \(M\) of the event

Returns:

float – Approximate source duration \(t\) (second)

relmt.utils.station_gap(station_dict, event_dict)

Azimuthal gap from center of event clud that would remain if station was removed

Parameters:

• station_dict (dict[str, Station]) – Station table

• event_dict (dict[int, Event]) – The seismic event catalog

Returns:

dict[str, float] – Mapping from station name to azimuthal gap (degree)

relmt.utils.subset_list(lst, indices)

Return subset of event dictionary

Parameters:

• event_list – Event list

• indices (Iterable[int]) – Indices of event_list to include in dictionary

• lst (list)

Returns:

list – Subset list

relmt.utils.valid_combinations(events, pairs, phase)

Indices to the event combinations that are in pairs

Parameters:

• pairs (set[tuple[int, int]]) – Event numbers (a, b) that may be combined pairwise. We assume sorted pairs: a < b

• events (list[int]) – List of event numbers

• phase (str) – Phase type, either ‘P’ or ‘S’

Returns:

ndarray – * Array of shape (combinations, 2) (if phase=’P’) or (combinations, 3) * (if phase=’S’) with the indices to the valid event combinations

relmt.utils.xyzarray(table)

Spatial coordinates as an array

Parameters:

table (dict[str, Station] | dict[int, Event] | Event | Station) – Event, event dictionary, station, or station dictionary

Returns:

ndarray – (rows, 3) or (3,) north-east-down coordinates