pyglimer.utils#

pyglimer.utils.create_geom#

Automatic creation of raysum geom files

copyright:: The PyGLImER 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, 14th May 2020 10:23:03 Last Modified: Monday, 30th May 2022 02:58:04 pm

pyglimer.utils.create_geom.create_geom(N: int, bazv: ndarray, raypv: ndarray, shift_max: int, filename: str, shape='cross')[source]#

Creates geometry files for Raysum.

Parameters:

N (int) – Number of stations. Has to be uneven if shape=cross. Else has to be N=M**2. Were M is a natural number
bazv (np.ndarray(1D)) – 1D array containing the backzimuths per station (deg).
raypv (np.ndarray(1D)) – 1D array containing the slownesses in s/m per backzimuth.
shift_max (int) – Maximum horizontal shift in m.
filename (str) – Name of the output file.
shape (str) – shape of the array

Raises:

ValueError – For Even Ns.

Return type:

None.

pyglimer.utils.createvmodel#

Create a 3D velocity model using Litho1.0

copyright:: The PyGLImER 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, 01st May 2020 12:11:03 Last Modified: Thursday, 21st October 2021 03:48:11 pm

pyglimer.utils.createvmodel.load_avvmodel() → AverageVelModel[source]#

Creates a model over the average P and S-wave velocities in the upper crust. These are used by the P-SV-SH rotation algorithm. Model data is extracted from the Litho1.0 model (location provided above). The package is distributed with a readily compiled model.

Litho1.0 must be installed and location must be set correct for a complete compilation! However, function will first look in RAM and then in data for a pickle file.

Compiling takes up to one hour!

Return type:: Model containing average velocities for upper crust.

pyglimer.utils.createvmodel.load_gyps(save=False, latb: Optional[tuple] = None, lonb: Optional[tuple] = None, flatten=True) → ComplexModel[source]#

Compiles the GyPSuM 3D-velocity object from included GyPSuM text files

Parameters:

save (Bool, optional) – Pickle the 3D velocity model after compiling it for the first time. This will allow for faster access to the model. Saving the model takes about 800 MB disk space. The default is False, as it lead to unstabilities with joblib.
latb (Tuple, optional) – Creates a submodel from the full model. In form (minlat, maxlat).
lonb (Tuple, optional) – (minlon, maxlon)

Returns:

Object that can be queried for velocities.

Return type:

ComplexModel object

pyglimer.utils.even2Dpoints#

pyglimer.utils.even2Dpoints.distance(x, y, x0: float, y0: float)[source]#

Compute euclidean 2D distance between points

Parameters:

x (arraylike) – x coordinates
y (arraylike) – y coordinates
x0 (float) – x coordinate to measure distance to
y0 (float) – y coordinatte to measure distance to

Returns:

distances

Return type:

float or arraylike

Notes

Author:: Lucas Sawade (lsawade@princeton.edu)
Last Modified:: 2021.03.16 16.30

pyglimer.utils.even2Dpoints.even2Dpoints(n, width, height, radius, seed=None)[source]#

Creates evenly distributed points in 2D space.

Parameters:

n (int) – number of points to be generated
width (floatt) – width of the 2D space
height (float) – height of the 2D space
radius (float) – min distance between the points

Returns:

tuple – x,y coordinates in 2D space
Does not really need a unit test, since its testing its distances
itself, and the distance function is really not a complicated thing…

pyglimer.utils.geo_utils#

Geographical utilities.

copyright:: The PyGLImER development team (makus@gfz-potsdam.de).
license:: GNU Lesser General Public License, Version 3 (https://www.gnu.org/copyleft/lesser.html)
author:: Lucas Sawade (lsawade@princeton.edu)

Last Update: November 2019

pyglimer.utils.geo_utils.cart2geo(x: ndarray, y: ndarray, z: ndarray) → Tuple[ndarray, ndarray, ndarray][source]#

Computes geographical coordinates from cartesian coordinates

Parameters:

arrays) (x (numpy.array corresponding to y/z) –
(numpy.array) (z) –
(numpy.array) –

Return type:

r, latitude, longitude tuple

pyglimer.utils.geo_utils.epi2euc(epi: float) → float[source]#: Converts epicentral distance in to a euclidean distance along the corresponding chord.

pyglimer.utils.geo_utils.euc2epi(euc: float) → float[source]#: Converts euclidean distance to epicentral distance

pyglimer.utils.geo_utils.gctrack(lat, lon, dist: float = 1.0, constantdist: bool = True) → Tuple[ndarray, ndarray, ndarray][source]#

Given waypoints and a point distance, this function computes evenly spaced points along the great circle and the given waypoints.

Warning

This function is not very accurate, but it is fast. You may want to resort to more accurate functions if that is of importance for your work. Here this function is only used for the purpose of making cross sections.

Warning

``constantdist`` The reason for this option is the fact that two waypoints are most likely not a multiple of dist apart from each other. A previous implementation included an interpolation between, but the interpolation of a parametric curve in 3D space is not as simple as I had hoped. The following explanation is for N segments defined by N+1 waypoints. constantdist sets the spacing constant, but the waypoint is most likely not exactly hit making the distance between the point before the (n+1)-th waypoint smaller than the the rest. This will be the same in all segments between n and (n+1)th waypoints. The alternative is to create linspace between the waypoints with N = round(total_dist/dist) + 1 points including n and (n+1) waypoints. This is essentially the more accurate way of defining the segments, but each segment has a different spacing that is not exactly dist.

Parameters:

lat (np.ndarray) – waypoint latitudes
lon (np.ndarray) – waypoint longitudes
dist (float) – distance in degrees
constantdist (bool) – Not really important for the enduser, but there are two ways of computing a gctrack along a curve. for 2 points only the most accurate way is to set constant_dist to False.

Return type:

tuple(np.ndarray, np.ndarray)

Notes

Authors:: Lucas Sawade (lsawade@princeton.edu)
Last Modified:: 2021.10.10 02.17

pyglimer.utils.geo_utils.geo2cart(r: float, latitude: ndarray, longitude: ndarray) → Tuple[ndarray, ndarray, ndarray][source]#

Computes cartesian coordinates from geographical coordinates

Parameters:

arrays) (r (float or numpy.array corresponding to lat/lon) –
location –
(numpy.array) (longitude) –
(numpy.array) –

Return type:

xyz coordinates coordinates of same size as latitude and longitude arrays

pyglimer.utils.geo_utils.geodiff(lat, lon)[source]#: Computes Azimuths and distances between geographical points.

pyglimer.utils.geo_utils.geodist(lat, lon)[source]#: Computes Azimuths and distances between geographical points.

pyglimer.utils.geo_utils.reckon(lat: float, lon: float, distance: float, bearing: float) → Tuple[float, float][source]#

Computes new latitude and longitude from bearing and distance.

Parameters:

lat (in degrees) –
lon (in degrees) –
bearing (in degrees) –
distance (in degrees) –

Return type:

lat, lon

lat1 = math.radians(52.20472) # Current lat point converted to radians lon1 = math.radians(0.14056) # Current long point converted to radians bearing = np.pi/2 # 90 degrees # lat2 52.20444 - the lat result I’m hoping for # lon2 0.36056 - the long result I’m hoping for.

pyglimer.utils.log#

copyright:: The PyGLImER 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, 13th September 2021 10:38:53 am Last Modified: Friday, 7th January 2022 11:05:39 am

pyglimer.utils.log.create_mpi_logger(logger: Logger, rank: int) → Logger[source]#

Creates a very similar logger to the input logger, but with name and filehandler dependent on mpi rank.

Parameters:: logger (logging.Logger) – The original logger to crete a rank-dependent version on.
Returns:: The rank dependent logger (different name and different file)
Return type:: logging.Logger

pyglimer.utils.log.start_logger_if_necessary(name: str, logfile: str, loglvl: int) → Logger[source]#

Initialise a logger.

Parameters:

name (str) – The logger’s name
logfile (str) – File to log to
loglvl (int) – Log level

Returns:

the logger

Return type:

logging.Logger

pyglimer.utils.nextpowof2#

copyright:: The PyGLImER 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, 18th October 2019 08:11:03 pm Last Modified: Thursday, 25th March 2021 03:56:56 pm

pyglimer.utils.nextpowof2.nextPowerOf2(n: int)[source]#: just returns the next higher power of two from n

pyglimer.utils.Ra2b#

pyglimer.utils.Ra2b.Ra2b(a, b)[source]#

Gets rotation matrix for R3 vectors that rotates a -> b. This is the linear algebra version of rodriquez formula. Theory: https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula

Parameters:

a (np.ndarray) – vector to rotate towards
b (np.ndarray) – vector to rotate from

Returns:

3x3 rotation matrix

Return type:

np.ndarray

Notes

Author:: Lucas Sawade (lsawade@princeton.edu)
Last Modified:: 2021.04.17 00.30

pyglimer.utils.roundhalf#

Rounds to next 0.5

copyright:: The PyGLImER 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: Saturday, 02nd May 2020 01:48:03 pm Last Modified: Thursday, 25th March 2021 04:00:02 pm

pyglimer.utils.roundhalf.roundhalf(number) → float[source]#

Rounds to next half of integer

Parameters:: number (float/int) – number to be rounded.
Returns:: Closest half of integer.
Return type:: float

pyglimer.utils.signalproc#

copyright:: The PyGLImER 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: Sunday, 20th October 2019 10:31:03 am Last Modified: Wednesday, 13th April 2022 02:55:06 pm

pyglimer.utils.signalproc.convf(u: ndarray, v: ndarray, nf: int, dt: float) → ndarray[source]#

Convolution conducted in the frequency domain.

Parameters:

u (np.array) – Array 1.
v (np.array) – Array 2.
nf (INTEGER) – Array length in frequency domain (use next power of 2).
dt (FLOAT) – Sampling Interval [s].

Returns:

c – Convolution of u with v.

Return type:

np.array

pyglimer.utils.signalproc.corrf(u: ndarray, v: ndarray, nf: int) → ndarray[source]#

Cross-correlation in frequency domain, calculates the x=crosscorr(u,v). Hence, v is the flipped vector.

Parameters:

u (np.array) – flipped vector.
v (np.array) – Vector that the correlation is measured to.
nf (INTEGER) – Array length in frequency domain (use next power of 2).

Returns:

x – Correlation vector.

Return type:

np.array

pyglimer.utils.signalproc.filter(s: ndarray, F: ndarray, dt: float) → ndarray[source]#

Convolves a filter with a signal (given in time domain).

Parameters:

s (np.array) – Signal given in time domain.
F (np.array) – Filter’s amplitude response (i.e., frequency domain).
dt (FLOAT) – Sampling interval [s].

Returns:

s_f – Filtered signal.

Return type:

np.array

pyglimer.utils.signalproc.gaussian(N: int, dt: float, width: float) → ndarray[source]#

Create a zero-phase Gaussian function (i.e., a low-pass filter). In particular meant to be convolved with the impulse response that is the output of an iterative deconvolution

Parameters:

N (INTEGER) – Length of the desired array.
dt (FLOAT) – Sampling interval [s].
width (FLOAT) – Width parameter of the Gaussian function.

Returns:

G – Gaussian window.

Return type:

np.array

pyglimer.utils.signalproc.noise(N: int, A: float) → ndarray[source]#: create random noise :param N: signal length :type N: int :param A: maximal amplitude :type A: float

pyglimer.utils.signalproc.resample_or_decimate(data: Stream, 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).

Note:

The stream has to be filtered a priori to avoid aliasing.

Parameters:

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

Returns:

The resampled stream

Return type:

Stream

pyglimer.utils.signalproc.ricker(sigma: float, N2: int, dt: float) → tuple[source]#

create a zero-phase Ricker / Mexican hat wavelet

Parameters:

sigma (float) – standard deviation
dt (float) – sampling Interval
N2 (int) – number of samples/2 (halflength). Full length = 2*N2+1

pyglimer.utils.signalproc.sshift(s: ndarray, N2: int, dt: float, shift: float) → ndarray[source]#

shift a signal by a given time-shift in the frequency domain

Parameters:

s (Arraylike) – signal
N2 (int) – length of the signal in f-domain (usually equals the next pof 2)
dt (float) – sampling interval
shift (float) – the time shift in seconds

Return type:

Timeshifted signal

pyglimer.utils.SphericalNN#

Spherical nearest neighbour idea taken from ObsPy originally, but heavily modified to include interpolation using shephard’s and modified shephard’s method of inverse distance weighting.

class pyglimer.utils.SphericalNN.SphericalNN(lat, lon)[source]#

Bases: object

Spherical nearest neighbour queries using scipy’s fast kd-tree implementation.

data#

cartesian point data array [x,y,z]

Type:: numpy.ndarray

kd_tree#

a KDTree used to query data

Type:: scipy.spatial.cKDTree

query(qlat, qlon)[source]#: Query a set of latitudes and longitudes

query_pairs(maximum_distance)[source]#: Find pairs of points that are within a certain distance of each other

interp(data, qlat, qlon)#: Use the kdtree to interpolate data corresponding to the points of the Kdtree onto a new set of points using nearest neighbor interpolation or weighted nearest neighbor interpolation (default).

Notes

Authors:: Lucas Sawade (lsawade@princeton.edu)
Last Modified:: 2021.04.21 20.00

Initialize class

Parameters:

lat (numpy.ndarray) – latitudes
lon (numpy.ndarray) – longitudes

interpolator(qlat, qlon, maximum_distance=None, no_weighting=False, k: Optional[int] = None, p: float = 2.0)[source]#

Spherical interpolation function using the SphericalNN object. Returns an interpolator that can be used for interpolating the same set of locations based on the KDTree. The only input the interpolator takes are the data corresponding to the points in the KDTree.

Parameters:

qlat (numpy.ndarray) – query latitudes
qlon (numpy.ndarray) – query longitude
maximum_distance (float, optional) – max distace for the interpolation in degree angle. Default None. If the mindistance to any points is larger than maximum_distance the interpolated value is set to np.nan.
no_weighting (bool, optional) – Whether or not the function uses a weightied nearest neighbor interpolation
k (int, optional) – Define maximum number of neighbors to be used for the weighted interpolation. Not used if no_weighting = True. Default None
p (float, optional) – Exponent to compute the inverse distance weights. Note that in the limit p->inf is just a nearest neighbor interpolation. Default is 2

Notes

In the future, I may add a variable weighting function for the weighted interpolation.

Please refer to https://en.wikipedia.org/wiki/Inverse_distance_weighting for the interpolation weighting.

query(qlat, qlon)[source]#

Query latitude and longitude values from kdtree

Parameters:

qlat (np.ndarray or float) – query latitude
qlon (np.ndarray or float) – query longitude

Returns:

(distance, indeces) to closest point in kdtree

Return type:

tuple

query_pairs(maximum_distance=180.0)[source]#

Query pairs within the kdtree

Parameters:: maximum_distance (float) – Maximum query distance in deg
Returns:: Set of pairs (i, j) where i < j
Return type:: set or ndarray

sparse_distance_matrix(other: Optional[SphericalNN] = None, maximum_distance=180.0, sparse: bool = False, km: bool = False)[source]#

Computes the sparse distance matrix between two kdtree. if no other kdtree is provided, this kdtree is used

Parameters:

other (SphericalNN) – other SphericalNN tree
maximum_distance (float) – Maximum query distance in deg
sparse (bool) – Flag to output sparse array for memory. Default False
km (bool) – flag to output distances in kilometers. Default False

Returns:

Set of pairs (i, j) where i < j

Return type:

set or ndarray

static spherical2cartesian(lat, lon)[source]#: Converts a list of Station objects to an array of shape(len(list), 3) containing x/y/z in meters.

pyglimer.utils.statistics#

Toolbox for statistical tools to quantify errors in receiver functions.

copyright:: The PyGLImER 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, 27th September 2021 03:14:17 pm Last Modified: Thursday, 21st October 2021 10:23:01 am

pyglimer.utils.statistics.stack_case_resampling(data: ndarray, b: int = 1000) → List[ndarray][source]#

Monte Carlo Case resampling algorithm that can be used with single station RF stacks (aka bootstrap).

Parameters:

data (np.ndarray) – Moveout corrected receiver function data that is used for the bootstrap. Given in form of a 2D matrix -> 1 line = 1 RF.
b (int, optional) – Number of iterations that the bootstrap is supposed to run. Defaults to 1000

Returns:

A list with each of the boot strapped stacks / averages.

Return type:

List[np.ndarray]

pyglimer.utils.utils#

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

Created: Tue May 26 2019 13:31:30 Last Modified: Monday, 26th September 2022 10:45:48 am

pyglimer.utils.utils.check_UTC_overlap(start: List[UTCDateTime], end: List[UTCDateTime]) → List[bool][source]#

Checks a list of starttimes and endtimes for overlap

Parameters:

start (List[UTCDateTime]) – List of UTCDatetime starttimes
end (List[UTCDateTime]) – List of UTCDateTime endtimes

Returns:

List of booleans. If True no overlap, if False overlap

Return type:

List[bool]

pyglimer.utils.utils.chunkdict(d: dict, chunksize) → list[source]#

Chunks upd dictionary with values into list of dictionaries where each dictionary has chunksize or less entries.

Parameters:

d (dict) – dictionary of lists
chunksize (int) – number of elements per list in dictionary

Returns:

returns a list of dictionaries subsets of lists in the original dictionary

Return type:

list

pyglimer.utils.utils.chunks(lst: list, n: int) → list[source]#: Yield successive n-sized chunks from lst. Useful for multi-threading

pyglimer.utils.utils.cos_taper(tr: Trace, taper_len: float, taper_at_masked: bool, side: str = 'both') → 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

Returns:

Tapered Trace

Return type:

Trace

Note

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

pyglimer.utils.utils.cos_taper_st(st: Stream, taper_len: float, taper_at_masked: bool, side: str = 'both') → 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

Returns:

Tapered Stream

Return type:

Stream

Note

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

pyglimer.utils.utils.create_bulk_str(networks: Union[str, List[str]], stations: Union[str, List[str]], location: str, channel: Union[str, List[str]], t0: UTCDateTime, t1: UTCDateTime) → List[tuple][source]#

Function to generate the input for the obspy functions: get_stations_bulk() and get_waveforms_bulk().

Parameters:

networks (str or List[str]) – The requested networks, can be str or List.
stations (str or List[str]) – The requested stations, can be str or List.
location (str) – Location string
channel (str) – Channel string
t0 (UTCDateTime or str or List[UTCDateTime]) – starttimes
t1 (UTCDateTime or str or List[UTCDateTime]) – endtimes

Raises:

ValueError – For invalid input types or inputs

Returns:

The List to be used as input with the aforementioned functions.

Return type:

List[tuple]

Note

All parameters accept wildcards.

pyglimer.utils.utils.download_full_inventory(statloc: str, fdsn_client: list)[source]#

This utility loops through statloc and redownloads the whole response (i.e., all channels and all times) for every station. Thus, overwriting the old xml file.

Parameters:

statloc (str) – Folder in which the old stationxmls are saved
fdsn_client (list) – List of FDSN providers that should be queried.

pyglimer.utils.utils.dt_string(dt: float) → str[source]#: Returns Time elapsed string depending on how much time has passed. After a certain amount of seconds it returns minutes, and after a certain amount of minutes it returns the elapsed time in hours.

pyglimer.utils.utils.get_multiple_fdsn_clients(clients: List[str]) → Tuple[Client][source]#

Returns a tuple holding all the queried fdsn providers. Also finds all available fdsn providers if input is None. Will also sort the clients so the big providers (IRIS and ORFEUS) come last.

Just a modified version of a code in the obspy mass downloader.

Parameters:: clients (List[str] or str or None) – List of strings, each describing one client. Put None if you want to use all available.
Returns:: Tuple of the fdsn Client objects.
Return type:: Tuple[Client]

pyglimer.utils.utils.save_raw(saved: dict, st: Stream, rawloc: str, inv: Inventory, saveh5: bool)[source]#

Save the raw waveform data in the desired format. The point of this function is mainly that the waveforms will be saved with the correct associations and at the correct locations.

Parameters:

saved (dict) – Dictionary holding information about the original streams to identify them afterwards.
st (Stream) – obspy stream holding all data (from various stations)
rawloc (str) – Parental directory (with phase) to save the files in.
inv (Inventory) – The inventory holding all the station information
saveasdf (bool) – If True the data will be saved in asdf format.

pyglimer.utils.utils.save_raw_mseed(evt: Event, slst: Stream, rawloc: str, net: str, stat: str)[source]#

Saves input stream in the correct location in mini seed format.

Parameters:

evt (Event) – event that is associated to the recording.
slst (Stream) – The selected stream. Should contain all components of one recording of a teleseismic event.
rawloc (str) – Location (with phase) to save the miniseed files in (will also create subfolders).
net (str) – Network code
stat (str) – Station Code

pyglimer.utils.utils.save_stream(st: Stream, rawloc: str, saveh5: bool, saved: dict, inv: Inventory, network: Optional[str] = None, station: Optional[str] = None)[source]#

Saves a stream to either mseed or HDF5 format depending on flag, and the event-channel dictionary saved which with created in a parent function.

Parameters:

st (Stream) – stream to save
rawloc (str) – path to database raw waveforms databases
saveh5 (bool) – whether to save hdf5 or mseed
saved (dict) – event-channel dictionary that contains the bulk requests made
inv (Inventory) – inventory
network (str | None, optional) – network string, by default None
station (str | None, optional) – station string, by default None