pysep.declust#

Event and station declustering and weighting algorithms for adjoint tomography.

Provides routines for downsizing an event catalog by performing a declustering algorithm to remove nearby stations, and also implements a weighting scheme to promote unique source receiver pairs.

authors:

adjTomo Dev Team (2023) Amanda McPherson (2021) Carl Tape (2018)

Module Contents#

Classes#

Declust

Declustering class in charge of declustering and source receiver weighting

class pysep.declust.Declust(cat, inv=None, data_avail=None, min_lat=None, max_lat=None, min_lon=None, max_lon=None)[source]#

Declustering class in charge of declustering and source receiver weighting

User-input parameters to determine algorithm behavior

Parameters:
  • cat (obspy.core.catalog.Catalog) – Catalog of events to consider. Events must include origin information latitude and longitude

  • inv (obspy.core.inventory.Inventory) – Inventory of stations to consider

  • data_avail (dict) – If None, Declust assumes that all events in cat were recorded by all stations in inv. This is typically not the case however, so this dict allows the user to tell Declust about data availability. Keys of data_avail must match resource IDs, and values must be lists of station names (NN.SSSS)

  • min_lat (float) – optional, minimum latitude for bounding box defining the region of interest. If not given, will use the minimum latitude in the catalog of events

  • max_lat (float) – optional, maximum latitude for bounding box defining the region of interest. If not given, will use the maximum latitude in the catalog of events

  • min_lon (float) – optional, minimum longitude for bounding box defining the region of interest. If not given, will use the minimum longitude in the catalog of events

  • max_lon (float) – optional, maximum longitude for bounding box defining the region of interest. If not given, will use the maximum longitude in the catalog of events

update_metadata(cat=None, inv=None)[source]#

Get metadata like location and event depth and magnitude from a given Cat and Inv and set as internal attributes. Needs to be as separate function as threshold_events will cut down the internal catalog representation so this will need to be re-run

threshold_catalog(zedges=None, min_mags=None, min_data=None)[source]#

Kick out events that fall below a given magnitude range or a given data availability range. Allow this to be done for various depth ranges or for the entire volume at once.

Note

Updates internal cat Catalog object and metadata in place

Parameters:
  • zedges (list of float) – depth [km] slices to partition domain into when thresholding data

  • min_mags (int or list) – a list of minimum magnitude thresholds for each depth slice. If zedges is None, should be a list of length==1, which provides minimum magnitude for entire catalog. Elif zedges is given, should be a list of len(zedges)-1, which defines minimum magnitude for each depth bin. For example if zedges=[0, 35, 400], then one example is min_mags=[4, 6]. Meaning between 0-34km the minimum magnitude is 4, and between 35-400km the minimum magnitude is 6.

  • min_data (int or list) – an integer or list of length len(zedges)-1 that defines the minimum number of stations on for a given event origin time, which allows user to prioritize stations with available data

calculate_srcrcv_weights(cat=None, inv=None, write='weights.txt', plot=False, show=False, save='srcrcvwght.png')[source]#

Calculate event and station specific weights based on event geographic weights, and event-dependent station geographic weights which take into account data availability for each event.

Uses the internal cat and inv attributes of the class. If declustering was run to thin out the catalog, run update_metadata first to update internal attributes which are used for defining weights.

Parameters:
  • cat (obspy.core.catalog.Catalog) – Optional, catalog of events to consider. If none given, will use whatever internal catalog is available

  • inv (obspy.core.inventory.Inventory) – Inventory of stations to consider

  • plot (bool) – plot source weights and an average of station weights on a single figure.

  • write (str) – filename used to write station weights to text file. if None, will not write

  • show (bool) – show figure in GUI

  • save (str) – if given, file id for the name of the output figure to save if not given, will not save figure

Return type:

np.array

Returns:

a 2D array where each row corresponds to a station and each column corresponds to an event. The value for a given row and column is the weight of that station w.r.t all other available stations. Stations that were not available (not on) are given a weight of 0

_get_weights(lons, lats, norm=None, plot=False, save='reference_distance_scan.png')[source]#

Given a set of coordinates (either stations or earthquakes), calculate a reference distance and a set of weights for each coordinate. Reference distance is chosen as one-third the largest possible value for all possible reference distances

Parameters:
  • lons (np.array) – array of longitude values

  • lats (np.array) – array of latitude values

  • plot (bool) – plot a simple scatterplot showing the weights for each station

  • norm (str) – how to normalize the weights - None: don’t normalize, provide raw weights - ‘max’: normalize by the maximum weight - ‘len’: normalize by the length of the array - ‘avg’: normalize by the mean weight value

Return type:

np.array

Returns:

relative, normalized weights for each lon/lat pair

static _covert_dist_to_weight(dists, ref_dist)[source]#

Calculate distance weights from a matrix of distances and a reference distance

Parameters:
  • dists (np.array) – array of inter-source distances

  • ref_dist (float) – user-defined reference distance in units km. larger values for reference distances increase the sensitivity to inter-station distances. lower values tend to reduce scatter.

decluster_events(cat=None, inv=None, choice='cartesian', zedges=None, min_mags=None, nkeep=1, select_by='magnitude', **kwargs)[source]#

Main logic function for choosing how to decluster events. Allow for both cartesian and polar binning of the domain.

See _decluster_events_cartesian and _decluster_events_polar for specific input parameters to control declustering.

Parameters:
  • cat (obspy.core.catalog.Catalog) – Optional, catalog of events to consider. If none given, will use whatever internal catalog is available

  • inv (obspy.core.inventory.Inventory) – Inventory of stations to consider

  • choice (str) – choice of domain partitioning, can be one of: - cartesian: grid the domain as a cube with nx by ny cells - polar: grid the domain with polar coordinates and ntheta bins

  • zedges (list of float) – depth [km] slices to partition domain into. Each slice will be given equal weighting w.r.t to all other slices, independent of slice size. e.g., allows upweighting crustal events

  • min_mags (list) – a list of minimum magnitude thresholds for each depth slice. If zedges is None, should be a list of length==1, which provides minimum magnitude for entire catalog. Elif zedges is given, should be a list of len(zedges)-1, which defines minimum magnitude for each depth bin. For example if zedges=[0, 35, 400], then one example is min_mags=[4, 6]. Meaning between 0-34km the minimum magnitude is 4, and between 35-400km the minimum magnitude is 6.

  • nkeep (int or list of int) – number of events to keep per cell. If zedges is None, then this must be an integer which defines a blanket value to apply. If zedges is given, then this must be a list of length zedges - 1, defining the number of events to keep per cell, per depth slice. See min_mags definition for example.

  • select_by (str) – determine how to prioritize events in the cell - magnitude (default): largest magnitudes prioritized - magnitude_r: smallest magnitudes prioritized - depth: shallower depths prioritized - depth_r: deeper depths prioritized - data: prioritize events which have the most data availability

_decluster_events_cartesian(nx=10, ny=10, zedges=None, nkeep=1, select_by='magnitude_r', plot=False, plot_dir='./', **kwargs)[source]#

Decluster event catalog by partitioning the 3D domain in the X, Y and Z directions, and then selecting a given number of events in each cell.

Parameters:
  • nx (int) – Number of X/longitude cells to partition domain into

  • ny (int) – Number of Y/latitude cells to partition domain into

  • zedges (list of float) – depth [km] slices to partition domain into. Each slice will be given equal weighting w.r.t to all other slices, independent of slice size. e.g., allows upweighting crustal events

  • nkeep (int or list of int) – number of events to keep per cell. If zedges is None, then this must be an integer which defines a blanket value to apply. If zedges is given, then this must be a list of length zedges - 1, defining the number of events to keep per cell, per depth slice. See min_mags definition for example.

  • select_by (str) – determine how to prioritize events in the cell - magnitude (default): largest magnitudes prioritized - magnitude_r: smallest magnitudes prioritized - depth: shallower depths prioritized - depth_r: deeper depths prioritized - data: less data availability prioritized - data_r: more data availability prioritized

  • plot (bool) – create a before and after catalog scatter plot to compare which events were kept/removed. Plots within the cwd

  • plot_dir (str) – directory to save figures to. file names will be generated automatically

Return type:

obspy.core.catalog.Catalog

Returns:

a declustered event catalog

_decluster_events_polar(ntheta=16, zedges=None, nkeep=1, select_by='magnitude_r', plot=False, plot_dir='./', **kwargs)[source]#

Run the declustering agorithm but partition the domain in polar. That is, divide each depth slice into a pie with ntheta partitions and keep events based on events within each slice of the pie. Option to cut each slice of pie by radius (distance from center of domain) and put additional constraints (e.g., more distant events require larger magnitude).

Parameters:
  • ntheta (int) – Number of theta bins to break a polar search into. Used to break up 360 degrees, so e.g., `ntheta`==17 will return bins of size 22.5 degrees ([0, 22.5, 45., 67.5 …. 360.])

  • zedges (list of float) – depth [km] slices to partition domain into. Each slice will be given equal weighting w.r.t to all other slices, independent of slice size. e.g., allows upweighting crustal events

  • nkeep (int or list of int) – number of events to keep per cell. If zedges is None, then this must be an integer which defines a blanket value to apply. If zedges is given, then this must be a list of length zedges - 1, defining the number of events to keep per cell, per depth slice. See min_mags definition for example.

  • select_by (str) – determine how to prioritize events in the cell - magnitude (default): largest magnitudes prioritized - magnitude_r: smallest magnitudes prioritized - depth: shallower depths prioritized - depth_r: deeper depths prioritized - data: less data availability prioritized - data_r: more data availability prioritized

  • plot (bool) – create a before and after catalog scatter plot to compare which events were kept/removed. Plots within the cwd

  • plot_dir (str) – directory to save figures to. file names will be generated automatically

Return type:

obspy.core.catalog.Catalog

Returns:

a declustered event catalog

plot(cat=None, inv=None, color_by='depth', connect_data_avail=False, vmin=None, vmax=None, title=None, cmap='inferno_r', show=True, save=None, equal_scale=False, **kwargs)[source]#

Geranalized plot function used to plot an event catalog and station inventory.

Parameters:
  • cat (obspy.core.catalog.Catalog) – Catalog of events to consider. Events must include origin information latitude and longitude

  • inv (obspy.core.inventory.Inventory) – Inventory of stations to consider

  • color_by (str) –

    how to color the event markers, available are - ‘depth’: color by the event depth - ‘data’: color by data availability for given event - ‘custom’: used by internal plotting routines to provide custom

    color array to the plot

  • connect_data_avail (bool) – connect sources and receivers with a thin line based on data availability

  • vmin (float) – min value for the colorbar, defaults to smallest value in array defined by color_by

  • vmax (float:) – maximum value for colorbar, defaults to largest value in the array defined by color_by

  • title (str) – custom user-input title for the figure, otherwise defaults to useful information about the catalog and inventory

  • cmap (str) – matplotlib colormap to use for array defined by color_by

  • show (bool) – show figure in GUI

  • save (str) – if given, file id for the name of the output figure to save if not given, will not save figure

  • equal_scale (bool) – set the scale of lat and lon equal, False by default

_plot_cartesian(cat, inv, xedges, yedges, title, save)[source]#

Convenience function to plot catalogs with edges

_plot_polar(cat, inv, evrad, theta_array, mid_lon, mid_lat, title, save)[source]#

Convenience function to plot catalogs with radial bin lines