Teleseismic Earthquakes (`seislib.eq`)

To retrieve earthquake-based phase-velocity measurements, SeisLib implements a two-station method [1]. The two-station method builds on the assumption that a fundamental-mode surface wave that has traveled large distances (~20°) from the epicenter of a strong earthquake (magnitudes $\gtrsim 5.5$) can be decomposed into a sum of monochromatics, or plane waves. This assumption is particularly convenient if we consider two receivers approximately lying on the same great-circle path as the epicenter. As the wavefront propagates from the first receiver to the second, its amplitude gets modified (by inter-station attenuation and/or site effects) and its phase gets shifted according to inter-station distance $\Delta$ and average velocity of propagation c. If attenuation is small, the phase shift can be expressed, in the frequency domain, as $\phi_2(\omega) - \phi_1(\omega) = \frac{\omega \Delta}{c(\omega)}$ [2], where the subscripts denote the receiver and $\phi$ the phase.

It is understood that the observed phase delay is invariant under $2\pi$ translations, giving rise to an ambiguity in phase velocity

\[c(\omega) = \frac{\omega \Delta}{\phi_2(\omega) - \phi_1(\omega) + 2n\pi},\]

where n is integer. In other words, the delay in the arrival time of the wavefront at the second receiver can be explained by different average inter-station phase velocities. Such ambiguity needs to be solved algorithmically, and SeisLib does it through the method extract_dispcurve() of seislib.eq.eq_velocity.TwoStationMethod.

References

1: Meier et al. (2004), One-dimensional models of shear wave velocity for the eastern Mediterranean obtained from the inversion of Rayleigh wave phase velocities and tectonic implications, GJI
2: Magrini et al. 2020, Arrival-angle effects on two-receiver measurements of phase velocity, GJI

Download of Teleseismic-Earthquake Recordings

The below class allows for automated downloads of seismograms recording teleseismic surface waves. The data will be saved on disk in the format required by seislib.eq.eq_velocity.EQVelocity and seislib.eq.eq_velocity.TwoStationMethod to calculate inter-station dispersion curves based on a two-station method.

class seislib.eq._eq_downloader.EQDownloader(savedir, inventory_name, inv_provider='iris', user=None, password=None, ev_provider='iris', vmin=1.5, vmax=5.5, sampling_rate=1, units='disp', prefilter=(0.001, 0.005, 0.05, 0.4), attach_response=False, stations_config={}, events_config={}, sleep_time_after_exception=30, verbose=True)

Bases: object

Downloads surface-wave data to be used via the two-station method

Parameters

savedirstr

Absolute path to the directory where the data will be saved. Inside this directory, a folder named ‘data’ will be created. Inside ‘data’ the seismograms will be organized in subdirectories (one for each seismic event) in the format net.sta.loc.cha.sac.

inventory_namestr

Name of the inventory (obspy.core.inventory.inventory.Inventory) associated with the downloads. This will be saved in xml format in the savedir

inv_providerstr

Provider of station metadata, passed to obspy.clients.fdsn.Client. Default is iris

user, passwordstr, optional

User name and password for accessing restricted data. Passed to obspy.clients.fdsn.Client

ev_providerstr

Provider of earthquake metadata, passed to obspy.clients.fdsn.Client. Default is iris

vmin, vmaxfloat

Minimum and maximum velocity of the surface waves. These values are used to establish the starting and ending time of the seismograms to be downloaded with respect to the origin time of given earthquake. Relatively loose limits are to be preferred. Default values are 1.5 and 5.5 (km/s), respectively

sampling_rateint, float

Final sampling rate of the waveforms, in Hz. Default is 1 Hz

units{‘DISP’, ‘VEL’, ‘ACC’}

Physical units of the final waveforms. Can be either ‘DISP’, ‘VEL’, or ‘ACC’ (case insensitive). Default is ‘DISP’.

prefilter(4, ) tuple of floats

Bandapass filter applied to the waveforms before removing the intrument response. Default is (0.001, 0.005, 0.05, 0.4)

attach_responsebool

If True, the details of the instrument response are attached to the header of the waveforms during their downloads. It can slow down the downloads, but it will make errors due to response information less likely. Default is False

stations_configdict

Python dictionary containing optional keys and values passed to obspy.clients.fdsn.Client.get_stations.

events_configdict

Python dictionary containing optional keys and values passed to obspy.clients.fdsn.Client.get_events. Default values are:

distmin = 2223.9 (km, corresponding to 20°)
distmax = 15567.3 (km, corresponding to 140°)
depthmax = 50 (km)
magmin = 6
magmax = 8.5
starttime = UTC(2000, 1, 1)
endtime = UTC(2021, 1, 1)

Warning

The obspy key words ‘minmag’, ‘maxmag’, and ‘maxdepth’ have been renamed as ‘magmin’, ‘magmax’, and ‘depthmax’, for better clarity in the code. Passing to events_config one of the (obspy) key words among ‘minmag’, ‘maxmag’, and ‘maxdepth’ will result in an error.

Note

The above listed default values include ‘distmin’ and ‘distmax’ (in km). These refer to the minimum and maximum distance of the epicenter from a given receiver. Note that these two key words correspond to the obspy’s ‘minradius’ and ‘maxradius’, with the only difference that the latters are expressed in degrees. ‘distmin’ and ‘distmax’ have been introduced only for achieving a higher consistency in the use of physical units throughout the code. The user should be aware that, although the use of obspy’s ‘minradius’ and ‘maxradius’ will not result in any error, their use is suggested against. In fact, if ‘distmin’ and ‘distmax’ are not specified, their default value will be used in the downloads. And if this default is more “restrictive” than ‘minradius’ and ‘maxradius’, ‘minradius’ and ‘maxradius’ will simpy be ignored in the downloads.

sleep_time_after_exceptionfloat

Time to wait (in s) after an obspy.clients.fdsn.header.FDSNException. The excpetion could be due to temporal issues of the client, so waiting a little bit before the next download could be useful to get things back to normality.

verbosebool

If True, information on the downloads will be printed in console

Examples

The following will download seismic data from all the seismic networks and stations managed by iris within the region specified by minlatitude, maxlatitude, minlongitude, maxlongitude. All such parameters should be passed to stations_config:

stations_config=dict(
        channel='LH*,BH*,HH*',
        includerestricted=False,
        maxlatitude=12,
        minlatitude=-18,
        minlongitude=90,
        maxlongitude=140)    

The above also specifies the preference of downloading 3-component seismograms (ZNE), as implied by the asterisk sign. For a given receiver recording an earthquake, either of the LH, BH, or HH channels will be downloaded, with priority given from left to right, i.e., LH>BH>HH: if LH is not available, BH will be downloaded; if BH is not available, HH will be downloaded.

In this example, we use all earthquakes characterized by a magnitude between 6 and 8.5, by a maximum depth of 50km, and generated between the 1.1.2000 and 1.1.2021. The distance between each earthquake and a given receiver should be more than 20° (expressed in km) and less than 140° (in km):

from obspy import UTCDateTime as UTC
from obspy.geodetics import degrees2kilometers

events_config=dict(
        starttime=UTC(2000, 1, 1),
        endtime=UTC(2021, 1, 1),
        depthmax=50,
        magmin=6, 
        magmax=8.5,
        distmin=degrees2kilometers(20),
        distmax=degrees2kilometers(140),
            )      

We initialize the EQDownloader instance, and then start it:

from seislib.eq import EQDownloader

downloader = EQDownloader(savedir='/path/to/directory',
                          inv_provider='iris',
                          ev_provider='iris',
                          inventory_name='iris.xml',
                          sampling_rate=1,
                          prefilter=(0.001, 0.005, 0.1, 0.4),
                          vmin=1.5,
                          vmax=5.5,
                          units='disp',
                          attach_response=False,
                          stations_config=stations_config,
                          events_config=events_config,
                          verbose=True)        
downloader.start()

Attributes

savedirstr: Absolute path of the directory that will be created
staclientobspy.clients.fdsn.Client: Provider of station metadata and waveforms
evclientobspy.clients.fdsn.Client: Provider of events metadata
channelslist: List of candidate channels for the downloads
components{3, 1}: Either 3 (ZNE) or 1 (Z)
vmin, vmaxfloat: Minimum and maximum expected velocity of the surface waves
sampling_ratefloat: Sampling rate of the final seismograms
prefiltertuple: Filter applied to the waveform before removal of the instrument response
unitsstr: Physical units of the final seismogram
attach_responsebool: Whether or not the response information is attached to the waveforms during the download
verbosebool: Whether or not progress information are displayed in the console
sleep_time_after_exceptionfloat: The downloads are stopped for the specified time (in s) after a FDSN exception
distmin, distmaxfloat: Minimum and maximum distance of a station from the epicenter (km)
depthmaxfloat: Maximum depth of the earthquakes (km)
magmin, magmaxfloat: Minimum and maximum magnitude
starttime, endtimeobspy.core.utcdatetime.UTCDateTime: Starttime and entime of the events catalogue
stations_config, events_configdict: Dictionary object with information relevant to the download of the stations and events metadata
inventoryobspy.core.inventory.inventory.Inventory: Stations metadata
_exceptionscollections.defaultdict(int): Dictionary object with information on the occurred exceptions
_downloadedint: Number of waveforms downloaded
_no_stationsint: Number of stations to download
_events_doneint: Number of events for which the downloads have been finished

Methods

`active_channels`(station)	Channels available for the given station among those to download
`adjust_channels`(st)	If the stream contains the Z12 channels, these are rotated towards ZNE.
`build_inventory`(**kwargs)	Builds an obspy inventory containing stations information
`collect_waveforms`(network, station, ...)	Downloads obspy stream
`compile_header_and_save`(st, savedir, stla, ...)	Compiles the header of the obspy stream (sac format) and writes to disk
`event_coordinates_and_time`(event)	Fetch event coordinates (lat, lon, depth) and origin time.
`fetch_catalog`(t1, t2, **kwargs)	Fetch catalog of seismic events
`get_event_info`(event)	Fetch event information
`handle_multiple_locations`(st, station_info)	Automatic selection of location.
`inventory_iterator`(inventory[, reverse])	Generator function to iterate over an obspy inventory
`plot_stations`([ax, show, oceans_color, ...])	Creates a maps of seismic receivers available for download
`prepare_data`(st)	Demean, detrend, tapering, removing response and resampling
`preprocessing`(st, station, baz)	Preprocessing of the obspy stream
`select_components`(st, baz)	Handles the absence of some components in the final stream
`station_coordinates`(station_info)	Fetch station coordinates (lat, lon, elevation).
`station_was_active`(station, time)	Wheater or not the seismic station was active at the given time

start

active_channels(station)

Channels available for the given station among those to download

Parameters

stationobspy.core.inventory.station.Station

Returns

List

adjust_channels(st)

If the stream contains the Z12 channels, these are rotated towards ZNE.

Parameters

stobspy.core.stream.Stream

Returns

obspy.core.stream.Stream if the rotation is successful, else None

build_inventory(**kwargs)

Builds an obspy inventory containing stations information

The strategy is to first download inventory information at station level. Then, for each station, the instrument response is downloaded for each channel. This may increase the downloading time, but prevents possible “timed out errors”.

Parameters

**kwargsdict, optional: Additional key word arguments passed to the get_stations method of obspy.clients.fdsn.client

Returns

invobspy.core.inventory.inventory.Inventory

collect_waveforms(network, station, channels, starttime, endtime)

Downloads obspy stream

Parameters

network, stationstr: network and station codes
channelsiterable of str: iterable containing the channels codes suited to the download. Each channel will be tried to be used in the downloads. The first successfull attempt determines the returned waveforms. (Higher priority is given to decreasing indexes)
starttime, endtimeobspy.core.utcdatetime.UTCDateTime: start and end time of the stream.

Returns

obspy.core.stream.Stream if the download is successful, else None

compile_header_and_save(st, savedir, stla, stlo, stel, evla, evlo, evdp, otime, mag, dist, az, baz)

Compiles the header of the obspy stream (sac format) and writes to disk

Parameters

streamobspy.core.stream.Stream
stla, stlo, stelfloat: Latitude, longitude, elevation of the seismic station
evla, evlo, evdpfloat: Latitude, longitude, depth of the earthquake
otimeobspy.core.utcdatetime.UTCDateTime: Origin time of the earthquake
mag, distfloat: Event magnitude and distance (km) of the event from the receiver
az, bazfloat: Azimuth and back azimuth of the epicenter with respect to the receiver

classmethod event_coordinates_and_time(event)

Fetch event coordinates (lat, lon, depth) and origin time.

Parameters

eventobspy.core.event.event.Event

Returns

(2,) Tuple: (latitude, longitude, depth) and origin time in obspy.UTCDateTime.timestamp format

fetch_catalog(t1, t2, **kwargs)

Fetch catalog of seismic events

Parameters

t1, t2obspy.core.utcdatetime.UTCDateTime: Starttime and endtime of the catalog
**kwargsdict, optional: Additional key-word arguments passed to obspy.clients.fdsn.client.Client.get_events

Returns

obspy.core.event.catalog.Catalog: Catalog of seismic events

classmethod get_event_info(event)

Fetch event information

Parameters

eventobspy.core.event.event.Event

Returns

(2, ) Tuple: (origin time, magnitude), (latitude, longitude, depth)

handle_multiple_locations(st, station_info)

Automatic selection of location.

Parameters

stobspy.core.stream.Stream
station_infoobspy.core.inventory.station.Station

Returns

stobspy.core.stream.Stream

classmethod inventory_iterator(inventory, reverse=False)

Generator function to iterate over an obspy inventory

Parameters

inventoryobspy.core.inventory.inventory.Inventory
reversebool: If True, the inventory will be iterated over from bottom to top

Yields

(2,) tuple: Yields network and station information at each iteration

plot_stations(ax=None, show=True, oceans_color='water', lands_color='land', edgecolor='k', projection='Mercator', resolution='110m', color_by_network=True, legend_dict={}, **kwargs)

Creates a maps of seismic receivers available for download

Parameters

stationsdict

Dictionary object containing stations information. This should structured so that each key corresponds to a station code ($network_code.$station_code) and each value is a tuple containing latitude and longitude of the station.

For example:

{ net1.sta1 : (lat1, lon1),
  net1.sta2 : (lat2, lon2),
  net2.sta3 : (lat3, lon3)
  }

axcartopy.mpl.geoaxes.GeoAxesSubplot

If not None, the receivers are plotted on the GeoAxesSubplot instance. Otherwise, a new figure and GeoAxesSubplot instance is created

showbool

If True, the plot is shown. Otherwise, a GeoAxesSubplot instance is returned. Default is True

oceans_color, lands_colorstr

Color of oceans and lands. The arguments are ignored if ax is not None. Otherwise, they are passed to cartopy.feature.NaturalEarthFeature (to the argument ‘facecolor’). Defaults are ‘water’ and ‘land’

edgecolorstr

Color of the boundaries between, e.g., lakes and land. The argument is ignored if ax is not None. Otherwise, it is passed to cartopy.feature.NaturalEarthFeature (to the argument ‘edgecolor’). Default is ‘k’ (black)

projectionstr

Name of the geographic projection used to create the GeoAxesSubplot. (Visit the cartopy website for a list of valid projection names.) If ax is not None, projection is ignored. Default is ‘Mercator’

resolution{‘10m’, ‘50m’, ‘110m’}

Resolution of the Earth features displayed in the figure. Passed to cartopy.feature.NaturalEarthFeature. Valid arguments are ‘110m’, ‘50m’, ‘10m’. Default is ‘110m’

color_by_networkbool

If True, each seismic network will have a different color in the resulting map, and a legend will be displayed. Otherwise, all stations will have the same color. Default is True

legend_dictdict, optional

Dictionary of keyword arguments passed to matplotlib.pyplot.legend

**kwargsdict, optional

Additional keyword arguments passed to matplotlib.pyplot.scatter

Returns

If show is True, None, else ax, i.e. the GeoAxesSubplot

prepare_data(st)

Demean, detrend, tapering, removing response and resampling

Parameters

stobspy.core.stream.Stream

Returns

obspy.core.stream.Stream if the processing is successful, else None

preprocessing(st, station, baz)

Preprocessing of the obspy stream

The function calls sequentially the methods handle_multiple_locations(), prepare_data(), and adjust_channels().

Parameters

stobspy.core.stream.Stream
stationobspy.core.inventory.station.Station
bazfloat: Back azimuth (degrees) of the epicenter with respect to the receiver

Returns

obspy.core.stream.Stream if the preprocessing is successful, else None

select_components(st, baz)

Handles the absence of some components in the final stream

Parameters

stobspy.Stream
bazint, float (in degrees): Back-azimuth used for the rotation NE->RT

Returns

stobspy.core.stream.Stream: If the expected three-component stream lacks of one or both the horizontal components, the vertical component is returned. If the stream lacks the vertical component but has the two horizontal ones, it returns the horizontal components rotated towards the back azimuth. If only one horizontal component is available, it returns None

start()

classmethod station_coordinates(station_info)

Fetch station coordinates (lat, lon, elevation).

Parameters

station_infoobspy.core.inventory.station.Station

Returns

(3,) Tuple of floats: latitude, longitude, elevation

classmethod station_was_active(station, time)

Wheater or not the seismic station was active at the given time

Parameters

stationobspy.core.inventory.station.Station
timeobspy.core.utcdatetime.UTCDateTime

Returns

bool: Whether or not the receiver was active at the specified time

Inter-Station Dispersion Curves

The below classes provide automated routines to calculate Rayleigh and Love inter-station dispersion curves based on teleseismic earthquakes, through the two-station method.

class seislib.eq._eq_velocity.EQVelocity(src, savedir=None, component='Z', verbose=True)

Bases: object

Class to obtain surface-wave phase velocities from teleseismic earthquakes, using the two-station method.

Parameters

srcstr: Absolute path to the directory containing the files associated with the earthquake recordings. The directory should contain, for each event, a sub-directory named after the origin time of the earthquake in timestamp format (see obspy.UTCDateTime.timestamp). Inside each sub-directory, the seismograms associated with the respective event recorded at all the available receivers should be stored into sacfiles. These should be named after the trace id, i.e., net.sta.loc.cha.sac (see obspy.Trace.id). Furthermore, inside each event directory there should be either an xml containing the event coordinates (latitude, longitude), or these coordinates should be stored in the sac header of the related seismograms.

Note

To download the data in the proper format, it is suggested to use seislib.eq.EQDownloader.
savedirstr, optional: Absolute path to the directory where the results are saved. If not provided, its value will be set to the parent directory of src (i.e., $src/..). All the results will be saved into the directory $savedir/eq_velocity/$component (see the component parameter)
component{‘Z’, ‘R’, ‘T’}: Either ‘Z’, ‘R’, (corresponding to vertically and radially polarized Rayleigh waves, respectively), or ‘T’ (for Love waves)
verbosebool: Whether or not information on progress is printed in the console

Examples

The following will calculate the dispersion curves from teleseismic earthquakes recorded at a set of stations. These are stored in the directory src=’/path/to/data’. The data have been downloaded using seislib.eq.eq_downloader.EQDownloader. The dispersion curves are extracted for vertically-polarized Rayleigh waves (on the vertical - Z - component). We first define the EQVelocity instance and prepare the data for the subsequent analysis:

>>> from seislib.eq import EQVelocity
>>> eq = EQVelocity(src, component='Z')
>>> eq.prepare_data(azimuth_tolerance=7, 
                    distmin=None, 
                    distmax=2500, 
                    min_no_events=8)
>>> print(eq)
RAYLEIGH-WAVE VELOCITY FROM TELESEISMIC EARTHQUAKES (TWO-STATION METHOD)
        =========================================================================
        EVENTS: X
        RECEIVERS: N.A.
        TRIPLETS (EPICENTER-RECEIVERS) AVAILABLE: N.A.
        DISPERSION CURVES RETRIEVED: Y
        =========================================================================
        SOURCE DIR: /path/to/data
        SAVE DIR: /path/to/eq_velocity/Z

The above will find, among our seismograms, station pairs approximately lying (azimuthal tolerance of $\pm7^{\circ}$) on the same great-circle path as the epicenters. We set a maximum inter-station distance of 2500 km, and all pairs of receivers for which less than 8 events are available for measuring the dispersion curves are discarded. Printing eq will give information on the data available. (In the above output, X and Y are integers depending on your data.)

Note

prepare_data() will extract the geographic coordinates of each seismic receivers from the header of the corresponding sac files, those of the seismic events, and store the information on the triplets of epicenters-receivers that will be used to obtain the dispersion curves.

The geographic coordinates of the seismic receivers are saved at $eq.savedir/stations.pickle, and are stored in a dictionary object where each key corresponds to a station code ($network_code.$station_code) and each value is a tuple containing latitude and longitude of the station. For example:

{ 
net1.sta1 : (lat1, lon1),
net1.sta2 : (lat2, lon2),
net2.sta3 : (lat3, lon3)
    }

The geographic coordinates of the events are saved at $eq.savedir/events.pickle, and have a similar structure to the above: each key corresponds to an event origin time (in the obspy.UTCDateTime.timestamp format), and each value is (lat, lon, mag) of the earthquake, where mag is magnitude.

The triplets of epicenter-receivers are saved at $eq.savedir/triplets.pickle, and are stored in a dictionary object where each key is a tuple corresponding to a station pair and each value is a list of events that are aligned on the same great circle path as the receivers (within the limits defined by the input params of prepare_data()). For example:

{ 
(net1.sta1, net2.sta2) : ['1222701570.84',
                          '1237378863.3', 
                          '1237486660.74',
                          '1238981562.62',
                          '1239825695.33',
                          etc.],

(net3.sta3, net4.sta4) : ['1222701570.84',
                          '1237378863.3',
                          '1237486660.74',
                          '1238981562.62',
                          '1239825695.33',
                          etc.]
    }

Note

Each event origin time corresponds to a sub-directory of the source data (src). Since the min_no_events passed to prepare_data() is 8, the length of the lists of events associated with each receiver pair will be at least 8.

To visualize the location of the epicenters that will be used to calculate inter-station dispersion curves, run, for example,:

eq.plot_events(color='r', 
               min_size=5, 
               max_size=250, 
               marker='*', 
               legend_dict=dict(fontsize=12))

We can now extract the dispersion curves in a automated fashion. The following will allow us to extract dispersion measurements at 75 different surface-wave periods, linearly spaced in the range 15-150 s. All values of phase velocity outside the velocity range 2.5-5 km/s will be discarded, and will only periods for which a ratio between inter-station distance and wavelength > 1 will be considered. (The wavelength is inferred from the reference curve.):

eq.extract_dispcurves(refcurve, 
                      periodmin=15, 
                      periodmax=150, 
                      no_periods=75, 
                      cmin=2.5, 
                      cmax=5, 
                      min_no_wavelengths=1, 
                      plotting=True)

Attributes

srcstr: Absolute path to the directory containing the files associated with the earthquake-based recordings
savedirstr: Absolute path to the directory where the results are saved
component{‘Z’, ‘R’, or ‘T’}
eventslist: List of available events to calculate the dispersion curves
verbosebool: Whether or not information on progress is printed in the console
stationsdict: Coordinates of the station pairs that can be employed to retrieve dispersion curves (can be accessed after calling the method prepare_data)
tripletsdict: Triplets of epicenters-receivers that can be employed to retrieve dispersion curves (can be accessed after calling the method prepare_data)

Methods

`delete_unused_files`()	Deletes every file in the data directory which is not useful for extracting dispersion curves (i.e., those waveform-files that are not included in triplets dict).
`extract_dispcurves`(refcurve[, periodmin, ...])	Automatic extraction of the dispersion curves for all available pairs of receivers.
`get_coords_and_triplets`(events[, ...])	Retrieves stations, events information, and the triplets of epicenter-receivers to be used to calculate the phase velocities
`get_events_used`()	Retrieves the events id for which triplets of epicenter-receivers are available to extract dispersion measurements
`interpolate_dispcurves`(period)	Interpolates the dispersion curves found at $self.savedir/dispcurves at the specified period(s).
`lie_on_same_gc`(stla1, stlo1, stla2, stlo2, ...)	Boolean function.
`plot_events`([ax, show, oceans_color, ...])	Creates a map of epicenters
`plot_stations`([ax, show, oceans_color, ...])	Maps the seismic receivers for which data are available
`prepare_data`([azimuth_tolerance, distmin, ...])	Saves to disk the geographic coordinates of the seismic receivers and of the seismic events, along with the triplets of epicenters-receivers to be used for retrieving the dispersion curves.
`prepare_input_tomography`(savedir, period[, ...])	Prepares a txt file for each specified period, to be used for calculating phase-velocity maps using `seislib.tomography.SeismicTomography`.

delete_unused_files(): Deletes every file in the data directory which is not useful for extracting dispersion curves (i.e., those waveform-files that are not included in triplets dict).

Warning

Use it with caution. It will not be possible to restore the deleted files.

extract_dispcurves(refcurve, periodmin=15, periodmax=150, no_periods=75, cmin=2.5, cmax=5, min_no_wavelengths=1.5, approach='freq', prob_min=0.25, prior_sigma_10s=0.7, prior_sigma_200s=0.3, manual_picking=False, shuffle=False, plotting=False, show=True)

Automatic extraction of the dispersion curves for all available pairs of receivers.

The results are saved to $self.savedir/dispcurves in the npy format, and consist of ndarrays of shape (n, 2), where the 1st column is period and the 2nd phase velocity (in m/s).

The routine iterates over all the available pair of receivers for which there are epicenters aligned on the same great circle path as the receivers (see the EQVelocity.prepare_data method); for each such pair of stations, (i) it first extracts dispersion measurements from all the event available, and then (ii) merges the dispersion measurements to obtain a “probability” density distribution of the thus retrieved dispersion measurements, which is function of period and phase velocity. (iii) Finally, the dispersion curve is extracted from the regions of greater “probability”. All this is done under the hood calling the measure_dispersion() and extract_dispcurve() of TwoStationMethod.

Parameters

refcurvendarray of shape (n, 2): Reference curve used to extract the dispersion curves. The 1st column should be period, the 2nd velocity (in either km/s or m/s). The reference curve is automatically converted to km/s, the physical unit employed in the subsequent analysis.
periodmin, periodmaxfloat: Minimum and maximum period analysed by the algorithm (default are 15 and 150 s). The resulting dispersion curves will be limited to this period range
no_periodsint: Number of periods between periodmin and periodmax (included) used in the subsequent analysis. The resulting periods will be equally spaced (linearly) from each other. Default is 75
cmin, cmaxfloat: Estimated velocity range spanned by the dispersion curves (default values are 2.5 and 5 km/s). The resulting dispersion curves will be limited to this velocity range
min_no_wavelengthsfloat: Ratio between the estimated wavelength $\lambda = T c_{ref}$ of the surface-wave at a given period T and interstation distance $\Delta$. Whenever this ratio is > min_no_wavelength, the period in question is not used to retrieve a dispersion measurement. Values < 1 are suggested against. Default is 1.5
approach{‘time’, ‘freq’}: Passed to TwoStationMethod.measure_dispersion(). It tells if the dispersion measurements are extracted in the frequency domain (‘freq’) or in the time domain (‘time’). Default is ‘freq’
prob_minfloat: Minimum acceptable “probability” in the density of dispersion measurements, at a given period, below which the dispersion curve is not picked. Larger values are more restrictive. Good values are between ~0.2 and ~0.35. Default is 0.25. See TwoStationMethod.extract_dispcurve().
prior_sigma_10s, prior_sigma_200sfloat: Standard deviations of the Gaussian functions built around the reference model (refcurve) at the periods of 10 and 200 s to calculate the prior probability of the dispersion measurements. At each analysed period, the standard deviation is interpolated (and eventually linearly extrapolated) based on these two values. Smaller values give more “weight” to the reference curve in the picking of the phase velocities. Defaults are 0.7 and 0.3. See TwoStationMethod.extract_dispcurve().
shufflebool: If True, the dispersion curves are calculated following a random order. This is indicated for large data sets that need to be processed through multiple processes to decrease the computational time.
plottingbool: If True, a figure is created for each retrieved dispersion curve. This is automatically displayed and saved in $self.savedir/figures
manual_pickingbool: If True, the user is required to pick the dispersion curve manually. The picking is carried out through an interactive plot.

get_coords_and_triplets(events, azimuth_tolerance=5, distmin=None, distmax=None)

Retrieves stations, events information, and the triplets of epicenter-receivers to be used to calculate the phase velocities

Parameters

eventslist: List of events (i.e., sub-directories in the data folder src)
azimuth_tolerancefloat: Maximum allowed deviation from the great circle path in degrees
distmin, distmaxfloat, optional: Minimum and maximum allowed inter-station distance (in km). Default is None

Returns

stationsdict

Each key corresponds to a station code ($network_code.$station_code) and each value is a tuple containing latitude and longitude of the station. For example:

{ net1.sta1 : (lat1, lon1),
  net1.sta2 : (lat2, lon2),
  net2.sta3 : (lat3, lon3)
  }

events_infodict

Each key corresponds to an event origin time and each value is a tuple containing latitude, longitude, and magnitude of the event:

{ '1222701570.84' : (lat1, lon1, mag1),
    '1237486660.74' : (lat2, lon2, mag2),
    '1239825695.33' : (lat3, lon3, mag3)
    }                 

tripletsdict

Each key is a tuple corresponding to a station pair and each value is a list of events that are aligned on the same great circle path as the receivers (within the limits defined by the input params). For example:

{ (net1.sta1, net2.sta2) : ['1222701570.84',
                            '1237378863.3', 
                            '1237486660.74',
                            '1238981562.62',
                            '1239825695.33'],

  (net3.sta3, net4.sta4) : ['1222701570.84',
                            '1237378863.3', 
                            '1237486660.74',
                            '1238981562.62',
                            '1239825695.33']
  }

Note that each event in the list corresponds to a sub-directory of the source data src.

get_events_used()

Retrieves the events id for which triplets of epicenter-receivers are available to extract dispersion measurements

Returns

events_useddict: Dictionary object where each key corresponds to an event (origin time in obspy.UTCDateTime.timestamp format, i.e., the name of the respective directory in self.src), and the associated values include all the station codes that exploit that event to extract a dispersion measurement

interpolate_dispcurves(period)

Interpolates the dispersion curves found at $self.savedir/dispcurves at the specified period(s). (No extrapolation is made.)

Parameters

periodfloat, array-like: Period (or periods) at which the dispersion curves will be interpolated

Returns

codesndarray of shape (n, 2): Codes associated with the station pairs for which a dispersion curve has been calculated
coordsndarray of shape (n, 4): Coordinates (lat1, lon1, lat2, lon2) of the station pairs corresponding to the station codes
measurementsndarray of shape (n, p): Phase velocity calculated for station pair contained in coords at the wanted period(s). p is the number of periods

Note

measurements could contain nans

Examples

Assume you calculated dispersion curves from your data, which are located in /path/to/data, and you had initialized your class:EQVelocity instance, to calculate inter-station dispersion curves, as:

from seislib.eq import EQVelocity

eq = EQVelocity(src=SRC, component='Z')
eq.prepare_data(azimuth_tolerance=7,
                distmin=None, 
                distmax=3000, 
                min_no_events=8)

Even if the computation of the dispersion curves is not yet finished, you can use the above instance to interpolate all the available dispersion curves as follows:

period = [20, 30, 40, 50]
codes, coords, measurements = eq.interpolate_dispcurves(period)  

classmethod lie_on_same_gc(stla1, stlo1, stla2, stlo2, evla, evlo, azimuth_tolerance=5, distmin=None, distmax=None)

Boolean function. If the station pair and the epicenter lie on the same great circle path, it returns True.

Parameters

stla1, stlo1float: Latitude and longitude of station 1
stla2, stlo2float: Latitude and longitude of station 2
evla, evlofloat: Latitude and longitude of the epicenter
azimuth_tolerancefloat: Maximum deviation from the great circle path in degrees
distmin, distmaxfloat, optional: Minimum and maximum allowed inter-station distance (in km). Default is None

Returns

bool

plot_events(ax=None, show=True, oceans_color='water', lands_color='land', edgecolor='k', projection='Mercator', resolution='110m', min_size=5, max_size=150, legend_markers=4, legend_dict={}, **kwargs)

Creates a map of epicenters

Parameters

lat, lonndarray of shape (n,): Latitude and longitude of the epicenters
magndarray of shape(n,), optional: If not given, the size of each on the map will be constant
axcartopy.mpl.geoaxes.GeoAxesSubplot: If not None, the receivers are plotted on the GeoAxesSubplot instance. Otherwise, a new figure and GeoAxesSubplot instance is created
showbool: If True, the plot is shown. Otherwise, a GeoAxesSubplot instance is returned. Default is True
oceans_color, lands_colorstr: Color of oceans and lands. The arguments are ignored if ax is not None. Otherwise, they are passed to cartopy.feature.NaturalEarthFeature (to the argument ‘facecolor’). Defaults are ‘water’ and ‘land’
edgecolorstr: Color of the boundaries between, e.g., lakes and land. The argument is ignored if ax is not None. Otherwise, it is passed to cartopy.feature.NaturalEarthFeature (to the argument ‘edgecolor’). Default is ‘k’ (black)
projectionstr: Name of the geographic projection used to create the GeoAxesSubplot. (Visit the cartopy website for a list of valid projection names.) If ax is not None, projection is ignored. Default is ‘Mercator’
resolution{‘10m’, ‘50m’, ‘110m’}: Resolution of the Earth features displayed in the figure. Passed to cartopy.feature.NaturalEarthFeature. Valid arguments are ‘110m’, ‘50m’, ‘10m’. Default is ‘110m’. Ignored if ax is not None.
min_size, max_sizefloat: Minimum and maximum size of the epicenters on the map. These are used to interpolate all magnitudes associated with each event, so as to scale them appropriately on the map. (The final “sizes” are passed to the argument s of matplotlib.pyplot.scatter)
legend_markersint: Number of markers displayed in the legend. Ignored if s (size of the markers in matplotlib.pyplot.scatter) is passed
legend_dictdict, optional: Keyword arguments passed to matplotlib.pyplot.legend
**kwargsdict, optional: Additional keyword arguments passed to matplotlib.pyplot.scatter

Returns

If show is True, None, else ax, i.e. the GeoAxesSubplot

plot_stations(ax=None, show=True, oceans_color='water', lands_color='land', edgecolor='k', projection='Mercator', resolution='110m', color_by_network=True, legend_dict={}, **kwargs)

Maps the seismic receivers for which data are available

Parameters

axcartopy.mpl.geoaxes.GeoAxesSubplot: If not None, the receivers are plotted on the GeoAxesSubplot instance. Otherwise, a new figure and GeoAxesSubplot instance is created
showbool: If True, the plot is shown. Otherwise, a GeoAxesSubplot instance is returned. Default is True
oceans_color, lands_colorstr: Color of oceans and lands. The arguments are ignored if ax is not None. Otherwise, they are passed to cartopy.feature.NaturalEarthFeature (to the argument ‘facecolor’). Defaults are ‘water’ and ‘land’
edgecolorstr: Color of the boundaries between, e.g., lakes and land. The argument is ignored if ax is not None. Otherwise, it is passed to cartopy.feature.NaturalEarthFeature (to the argument ‘edgecolor’). Default is ‘k’ (black)
projectionstr: Name of the geographic projection used to create the GeoAxesSubplot. (Visit the cartopy website for a list of valid projection names.) If ax is not None, projection is ignored. Default is ‘Mercator’
resolution{‘10m’, ‘50m’, ‘110m’}: Resolution of the Earth features displayed in the figure. Passed to cartopy.feature.NaturalEarthFeature. Valid arguments are ‘110m’, ‘50m’, ‘10m’. Default is ‘110m’. Ignored if ax is not None.
color_by_networkbool: If True, each seismic network will have a different color in the resulting map, and a legend will be displayed. Otherwise, all stations will have the same color. Default is True
legend_dictdict, optional: Dictionary of keyword arguments passed to matplotlib.pyplot.legend
**kwargsdict, optional: Additional keyword arguments passed to matplotlib.pyplot.scatter

Returns

If show is True, None, else ax, i.e. the GeoAxesSubplot

prepare_data(azimuth_tolerance=5, distmin=None, distmax=None, min_no_events=5, recompute=False, delete_unused_files=False)

Saves to disk the geographic coordinates of the seismic receivers and of the seismic events, along with the triplets of epicenters-receivers to be used for retrieving the dispersion curves.

Parameters

azimuth_tolerancefloat: Maximum allowed deviation from the great circle path in degrees. All triplets of epicenter-receivers for which the receivers are not aligned within the tolerance indicated are rejected. Larger values will identify more triplets to be used in the following analysis. But if this value is too large, the assumptions behind the two- station method [1] may not be met. Suggested values are between 3 and 8. Default is 5.
distmin, distmax: float, optional: Minimum and maximum allowed inter-station distance (in km). Default is None
min_no_eventsint: Minimum number of events available for a given station pair to be considered in the calculation of the phase velocities.
recomputebool: If True, the station coordinates and triplets will be removed from disk and recalculated. Otherwise (default), if they are present, they will be loaded into memory, avoiding any computation. This parameter should be set to True whenever one wants to change the other parameters of this function, which control the selection of the epicenter-receivers triplets
delete_unused_filesbool: If True, every waveform-file that is not contained in the triplets object (i.e., those that are not used to extract dispersion curves in the subsequent analysis) will be permanently deleted from the system.

Notes

The geographic coordinates of the seismic receivers are saved at $self.savedir/stations.pickle, and are stored in a dictionary object where each key corresponds to a station code ($network_code.$station_code) and each value is a tuple containing latitude and longitude of the station. For example:

{ net1.sta1 : (lat1, lon1),
  net1.sta2 : (lat2, lon2),
  net2.sta3 : (lat3, lon3)
  }

The geographic coordinates of the events are saved at $self.savedir/events.pickle, and have a similar structure to the above: each key corresponds to an event origin time (in obspy.UTCDateTime.timestamp format), and each value is (lat, lon, mag) of the epicenter, where mag is the magnitude of the event.

The triplets of epicenter-receivers are saved at $self.savedir/triplets.pickle, and are stored in a dictionary object where each key is a tuple corresponding to a station pair and each value is a list of events that are aligned on the same great circle path as the receivers (within the limits defined by the input params). For example:

{ (net1.sta1, net2.sta2) : ['1222701570.84',
                            '1237378863.3', 
                            '1237486660.74',
                            '1238981562.62',
                            '1239825695.33'],

  (net3.sta3, net4.sta4) : ['1222701570.84',
                            '1237378863.3', 
                            '1237486660.74',
                            '1238981562.62',
                            '1239825695.33']
  }

Note that each event in the list corresponds to a sub-directory of the source data src.

References

1: Magrini et al. 2020, Arrival-angle effects on two-receiver measurements of phase velocity, GJI

prepare_input_tomography(savedir, period, min_no_wavelengths=1.5, outfile='input_%.2fs.txt')

Prepares a txt file for each specified period, to be used for calculating phase-velocity maps using seislib.tomography.SeismicTomography.

Parameters

savedirstr: Absolute path to the directory where the file(s) is (are) saved. If savedir does not exist, it will be created
periodfloat, array-like: Period (or periods) at which the dispersion curves will be interpolated using interpolate_dispcurves()
min_no_wavelengthsfloat: Ratio between the estimated wavelength $\lambda = T c_{ref}$ of the surface-wave at a given period T and interstation distance $\Delta$. Whenever this ratio is > min_no_wavelength, the period in question is not used to retrieve a dispersion measurement. Values < 1 are suggested against. Default is 1.5
outfilestr: Format for the file names. It must include either %s or %.Xf (where X is integer), since this will be replaced by each period analysed (one for file)

Examples

Assume you calculated dispersion curves from your data, which are located in /path/to/data, and you had initialized your class:EQVelocity instance, to calculate inter-station dispersion curves, as:

from seislib.eq import EQVelocity
eq = EQVelocity(src=SRC, component='Z')
eq.prepare_data(azimuth_tolerance=7,
                distmin=None, 
                distmax=3000, 
                min_no_events=8)

You can, even if the computation of the dispersion curves is not finished, use the above instance to extract the information needed to produce phase-velocity maps as follows:

savedir = /arbitrary/path/to/savedir
periods = [20, 30, 40, 50, 75, 100]
eq.prepare_input_tomography(savedir=savedir,
                            period=periods)

The above will save one txt file for each specified period.

class seislib.eq._eq_velocity.TwoStationMethod(refcurve, periodmin=15, periodmax=150, no_periods=75, cmin=2.5, cmax=5, ttol=0.3, min_no_wavelengths=1.5, approach='freq', gamma_f=14, gamma_w=None, distances=None)

Bases: object

Class providing methods to process and extract dispersion curves from earthquake recordings generated by strong teleseismic earthquakes.

Parameters

refcurvendarray of shape (n, 2)

Reference curve used to extract the dispersion curves. The 1st column should be period, the 2nd velocity (in either km/s or m/s). The reference curve is automatically converted to km/s, the physical units employed in the subsequent analysis.

periodmin, periodmaxfloat

Minimum and maximum period analysed by the algorithm (default are 15 and 150 s). The resulting dispersion curves will be limited to this period range

no_periodsint

Number of periods between periodmin and periodmax (included) used in the subsequent analysis. The resulting periods will be equally spaced (linearly) from each other. Default is 75

cmin, cmaxfloat

Estimated velocity range spanned by the dispersion curves (default values are 2.5 and 5 km/s). The resulting dispersion curves will be limited to this velocity range

ttolfloat

Tolerance, with respect to the reference curve, used to taper the seismograms around the expected arrival time of the surface wave (at a given period). In practice, at a given period, everything outside of the time range given by tmin and tmax (see below) is set to zero through a taper. tmin and tmax are defined as:

tmin = dist / (ref_vel + ref_vel*ttol)
tmax = dist / (ref_vel - ref_vel*ttol)

where dist is inter-station distance. Default is 0.3, i.e., 30% of the reference velocity

min_no_wavelengthsfloat

Ratio between the estimated wavelength $\lambda = T c_{ref}$ of the surface-wave at a given period T and interstation distance $\Delta$. Whenever this ratio is > min_no_wavelength, the period in question is not used to retrieve a dispersion measurement. Values < 1 are suggested against. Default is 1.5

approach{‘freq’, ‘time’}

Passed to TwoStationMethod.measure_dispersion(). It tells if the dispersion measurements are extracted in the frequency domain (‘freq’) or in the time domain (‘time’). Default is ‘freq’

gamma_ffloat

Controls the width of the bandpass filters, at a given period, used to isolate the fundamental mode in the seismogram [1].

gamma_w, distancesndarray of shape (m,), optional

Control the width of tapers used to taper, at a given period, the cross correlations in the frequency domain (these two parameters are ignored if approach is ‘time’). distances should be in km. If not given, they will be automatically set to:

gamma_w = np.linspace(5, 50, 100)
distances = np.linspace(100, 3000, 100)

These two arrays are used as interpolators to calculate gamma based on the inter-station distance. gamma is the parameter that actually controls the width of the tapers [1], and is defined as:

gamma = np.interp(dist, distances, gamma_w)       

References

1(1,2): Soomro et al. (2016), Phase velocities of Rayleigh and Love waves in central and northern Europe from automated, broad-band, interstation measurements, GJI

Examples

In the following we will calculate a dispersion curve based on several earthquakes recorded at two receivers (sta1, sta2) on the vertical components (Z). The seismograms (two per earthquake) are stored in as many directories as the number of earthquakes available in the analysis. These directories are named after the event origin time (in the obspy.UTCDateTime.timestamp format) and located at /path/to/data. The seismic data must be in the .sac format, with the information on the epicentral distance (dist) compiled in the header. It is assumed that receivers and epicenter lie on the same great-circle path.

We will retrieve dispersion measurements for 75 periods linearly spaced between 15 and 150s. We assume that in this period range surface waves travel in the velocity range 2.5 - 5 km/s. We need to pass a reference curve consisting of 2 columns (1st: period, 2nd: phase velocity in either km/s or m/s), that will be used in the subsequent analysis. We will only employ periods for which the ratio between expected wavelength (based on the reference curve) and inter-station distance is not smaller than 1.5.

First, we initialize the TwoStationMethod class as:

from seislib.eq import TwoStationMethod

tsm = TwoStationMethod(refcurve=refcurve,
                       periodmin=15,
                       periodmax=150,
                       no_periods=75,
                       cmin=2.5, 
                       cmax=5, 
                       min_no_wavelengths=1.5)

Then, we use this class obtain dispersion measurements for all events available in the data directory. (We are assuming that we only have folders corresponding to events that are aligned on the same great circle path as the epicenter. For a more automatic processing of station pairs and events, see EQVelocity.prepare_data().) We will store all the dispersion measurements in the folder /path/to/savedir/dispersion, in the npy format:

import os
import numpy as np
from obspy import read

src = /path/to/data
dispersion_dir = /path/to/savedir/dispersion

events = os.listdir(src)
for origin_time in events:
    st1 = read(os.path.join(src, origin_time, '%s*'%sta1))
    st2 = read(os.path.join(src, origin_time, '%s*'%sta2))    
    tsm.preprocess(st1, st2, float(origin_time))
    dispersion = tsm.measure_dispersion()
    outfile = os.path.join(dispersion_dir, 
                           '%s_%s_%s.npy'%(sta1, sta2, origin_time))
    np.save(outfile, dispersion)

Now that all the dispersion measurements have been extracted, we can calculate a dispersion curve based on them. (dist_km is inter-station distance in km). The results will displayed in the console:

dispcurve = tsm.extract_dispcurve(refcurve=refcurve,
                                  src=dispersion_dir,
                                  dist_km=dist_km,
                                  plotting=True,
                                  sta1=sta1,
                                  sta2=sta2)

Attributes

periodsndarray of shape (n,): Array of periods to extract the dispersion measurements from
ref_velndarray of shape (n,): Array of reference velocities corresponding to periods
ref_modelscipy.interpolate.interpolate.interp1d: Interpolating function. Input: period. Returns: phase velocity
min_no_wavelegnthsfloat: Minimum number of wavelengths between two receivers, at a given period, below which the dispersion measurements are not extracted
cmin, cmaxfloat: Estimated velocity range (in km/s) spanned by the dispersion curves
ttolfloat: Tolerance (percentage), with respect to the reference curve, used to taper the seismograms around the expected arrival time of the surface wave (at a given period).
gamma_ffloat, int: Controls the width of the bandpass filters, at a given period, used to isolate the fundamental mode in the seismogram.
gamma_w, distancesndarrays of shape (m,), optional: Control the width of tapers used to taper, at a given period, the cross correlations in the frequency domain
approachstr: It indicates if the dispersion measurements are extracted in the frequency domain (‘freq’) or in the time domain (‘time’)
otimefloat: Origin time of the earthquake. Available after calling TwoStationMethod.preprocess()
dist1, dist2float: Epicentral distances at the two receivers. Available after calling TwoStationMethod.preprocess()
distfloat: Inter-station distance, calculated as dist2-dist1. Available after calling TwoStationMethod.preprocess()
dtfloat: Time sampling of the seismograms. Available after calling TwoStationMethod.preprocess()
fsfloat: Sampling rate of the seismograms. Available after calling TwoStationMethod.preprocess()
data1, data2ndarray of shape (n,): Preprocessed seismograms. Available after calling TwoStationMethod.preprocess()
timesndarray of shape (n,): Array of times corresponding to the two seismograms, starting at the origin time. Available after calling TwoStationMethod.preprocess()
periods_maskedndarray of shape (n,): Subset of periods, where the wavelength is such that there is a number of wavelenghts between the two receivers > min_no_wavelength. Available after calling TwoStationMethod.preprocess()

Methods

`adapt_startime_to_t`(data, dt, starttime, ...)	Zero-pads data in order to make the start time coincide with the event origin time
`adapt_times`(tr1, tr2)	Zero-pads and slices the two traces so that their startting and ending times coincide
`build_taper`(center_idx, taper_size, data_size)	Defines the taper used for tapering the seismograms
`convert_to_kms`(dispcurve)	Converts from m/s to km/s the dispersion curve (if necessary).
`extract_dispcurve`(refcurve, src, dist_km[, ...])	Extraction of the dispersion curve from a given set of dispersion measurements.
`freq_domain_dispersion`([taper_p])	Extract dispersion measurements using a frequency-domain approach.
`frequency_w_and_alphaf`(period)	Retrieves frequency, angular frequency, and alpha_f at a given period
`measure_dispersion`([taper_p])	Extract dispersion measurements using either a time-domain or a frequency-domain approach, depending on the parameter approach passed to `TwoStationMethod`.
`preprocess`(st1, st2, otime[, fs])	Read sac infos and prepare to pv extraction.
`taper_from_times`(tmin, tmax[, taper_p])	Retrieves the tapers to be applied to the seismograms
`time_domain_dispersion`([taper_p])	Extract dispersion measurements using a time-domain approach.
`times_for_taper`(dist, ref_vel, tolerance)	Identifies the starttime and endtime of the tapers to be applied to the seismograms

adapt_startime_to_t(data, dt, starttime, origintime, power_2_sized=True)

Zero-pads data in order to make the start time coincide with the event origin time

Parameters

datandarray
dtfloat: Time sampling interval
starttime, origintimeobspy.UTCDateTime.timestamp: Starttime and origintime of the waveform
power_2_sizedbool: If True, the returned data size will be a power of 2. This allows to take the Fourier transform more efficiently

Returns

numpy.ndarray: Padded data

adapt_times(tr1, tr2)

Zero-pads and slices the two traces so that their startting and ending times coincide

Parameters

tr1, tr2obspy.Trace

Returns

tr1, tr2obspy.Trace: Modified copy of the original input

build_taper(center_idx, taper_size, data_size, taper_type=<function tukey>, alpha=0.1)

Defines the taper used for tapering the seismograms

Parameters

center_idxint: Index of corresponding to the center of the taper
taper_sizeint: Lenght of the taper. The taper will be centred onto taper_idx and extend around the center by taper_size/2
data_size :: Length of the data, which corresponds to the final size of the taper
taper_typefunc: Taper function. Default is scipy.signal.windows.tukey
alphafloat: Shape parameter of the taper window.

Returns

tapernumpy.ndarray

classmethod convert_to_kms(dispcurve)

Converts from m/s to km/s the dispersion curve (if necessary).

Parameters

dispcurvendarray of shape (n, 2): The 1st column (typically frequency or period) is ignored. The 2nd column should be velocity. If the first value of velocity divided by 10 is greater than 1, the 2nd column is divided by 1000. Otherwise, the dispersion curve is left unchanged.

Returns

dispcurvendarray of shape (n, 2): Dispersion curve eventually converted to km/s

classmethod extract_dispcurve(refcurve, src, dist_km, prob_min=0.25, smoothing=0.3, min_derivative=- 0.01, prior_sigma_10s=0.7, prior_sigma_200s=0.3, plotting=False, savefig=None, sta1=None, sta2=None, manual_picking=False, show=True)

Extraction of the dispersion curve from a given set of dispersion measurements. These (one per earthquake) should be stored in a directory as npy files.

The algorithm will (i) first gather all the dispersion measurments in a unique ndarray of shape (n, 2), where the 1st column is period and the 2nd phase velocity. From this ensemble, (ii) it will then create a 2-D image representing the density of measurements, via a method similar to a kernel-density estimate (KDE). This image (of shape (k, l) and henceforth referred to as P_obs), is normalized at each period by the maximum at the same period. (iii) A second image, call it P_ref, is obtained from the reference curve: the weights given to each pixel at a given period are defined by a Gaussian function centred onto the reference velocity at the same period, so that they decrease with distance from the reference velocity. As in P_obs, the image has a shape of (k, l) and is normalized at each period by the maximum. (iv) Element- wise multiplication of P_ref and P_obs yields a third image, call it P_cond, which can be interpreted as the probability to measure a dispersion curve at given period and velocity given the dispersion observations and the a-priori knowledge of the reference curve. (v) The dispersion curve is picked starting from the longest periods, where the phase ambiguity is less pronounced, and is driven by two quality parameters: min_derivative (so as to avoid strong velocity decreases with period) and prob_min (areas in P_cond characterized by values smaller than prob_min are not considered).

Parameters

refcurvendarray of shape (n, 2)

Reference curve used to extract the dispersion curves. The 1st column should be period, the 2nd velocity (in either km/s or m/s). The reference curve is automatically converted to km/s, the physical unit employed in the subsequent analysis.

srcstr

Absolute path to the directory containing the dispersion measurements

dist_kmfloat

Inter-station distance (in km)

prob_minfloat

“probability” in P_cond, at a given period, below which the dispersion curve is not picked. Larger values are more restrictive. Good values are between ~0.2 and ~0.35. Default is 0.25

prior_sigma_10s, prior_sigma_200sfloat

Standard deviations of the Gaussians built around the reference model (refcurve) at the periods of 10 and 200 s to calculate P_ref. At each analysed period, the standard deviation is interpolated (and eventually linearly extrapolated) based on these two values. Smaller values give more “weight” to the reference curve in the picking of the phase velocities. Defaults are 0.7 and 0.3

smoothingfloat

Final smoothing applied to the dispersion curve. The smoothing is performed using scipy.signal.savgol_filter. This parameters allows for defining the window_length passed to the SciPy function, via:

window_length = int(len(dispcurve) * smoothing)

Default is 0.3

min_derivativefloat

Minimum derivative of phase-velocity as a function of period. If exceeded, a quality selection will truncate the dispersion curve. To avoid the presence of gaps in the final dispersion curves, if they are present due to this truncation, the branch of dispersion curve associated with the largest probability (calculated from P_cond) is returned, the other(s) are discarded.

Default is -0.01. Smaller (negative) values are less restrictive.

plottingbool

If True, a figure is created for each retrieved dispersion curve

savefigstr, optional

Absolute path to the output image. Ignored if plotting is False

sta1, sta2str, optional

Station codes of the two receivers used to calculate the dispersion curve. If plotting is True, they are displayed in title of the title of the resulting figure. Otherwise they are ignored

manual_pickingbool

If True, the user is required to pick each dispersion curve manually. The picking is carried out through an interactive plot

showbool

If False, the figure is closed before being shown. This is only relevant when the user has the interactive plot enabled: in that case, the figure will not be displayed. Default is True

Returns

dispcurvendarray of shape (m, 2): Dispersion curve, where he 1st column is period and the 2nd is phase velocity in km/s

Raises

DispersionCurveException: If the calculation of the dispersion curve fails.

freq_domain_dispersion(taper_p=0.2)

Extract dispersion measurements using a frequency-domain approach.

Parameters

taper_pfloat: Taper percentage

Returns

solutionsndarray of shape (n, 2): The 1st column is period, the 2nd is phase velocity. Note that, in general, at a given period more than one phase-velocity measurement will be present in solutions, because of the $2\pi$ phase ambiguity [1]. All phase velocities outside the velocity range cmin-cmax are discarded (see TwoStationMethod).

References

1: Magrini et al. 2020, Arrival-angle effects on two-receiver measurements of phase velocity, GJI

frequency_w_and_alphaf(period)

Retrieves frequency, angular frequency, and alpha_f at a given period

Parameters

periodfloat

Returns

frequencyfloat: Inverse of period
wfloat: Angular frequency, i.e. $2\pi \times$ frequency
alpha_ffloat: Controls the width of the bandpass filters, at a given period, used to isolate the fundamental mode in the seismogram [1]. It is a function of gamma_f (See TwoStationMethod).

References

1: Soomro et al. (2016), Phase velocities of Rayleigh and Love waves in central and northern Europe from automated, broad-band, interstation measurements, GJI

measure_dispersion(taper_p=0.2)

Extract dispersion measurements using either a time-domain or a frequency-domain approach, depending on the parameter approach passed to TwoStationMethod. The function will call either time_domain_dispersion() or freq_domain_dispersion(), respectively.

Parameters

taper_pfloat: Taper percentage
solutionsndarray of shape (n, 2): The 1st column is period, the 2nd is phase velocity. Note that, in general, at a given period more than one phase-velocity measurement will be present in solutions, because of the $2\pi$ phase ambiguity [1]. All phase velocities outside the velocity range cmin-cmax are discarded (see TwoStationMethod).

References

1: Magrini et al. 2020, Arrival-angle effects on two-receiver measurements of phase velocity, GJI

preprocess(st1, st2, otime, fs=2)

Read sac infos and prepare to pv extraction. The sac headers must be compiled with traces and event information

Parameters

st1, st2obspy.core.stream.Stream
otimeobspy.UTCDateTime.timestamp: Event origin time
fsint: Sampling rate
taper_pfloat: Percentage taper to apply to the streams.

taper_from_times(tmin, tmax, taper_p=0.2)

Retrieves the tapers to be applied to the seismograms

Parameters

tmin, tmaxfloat: Starting and ending time of the taper
taper_pfloat: Shape parameter of the taper. Default is 0.2

Returns

taperndarray of shape (n,)

time_domain_dispersion(taper_p=0.2)

Extract dispersion measurements using a time-domain approach.

Parameters

taper_pfloat: Taper percentage

Returns

solutionsndarray of shape (n, 2): The 1st column is period, the 2nd is phase velocity. Note that, in general, at a given period more than one phase-velocity measurement will be present in solutions, because of the $2\pi$ phase ambiguity [1]. All phase velocities outside the velocity range cmin-cmax are discarded (see TwoStationMethod).

References

1: Magrini et al. 2020, Arrival-angle effects on two-receiver measurements of phase velocity, GJI

times_for_taper(dist, ref_vel, tolerance)

Identifies the starttime and endtime of the tapers to be applied to the seismograms

Parameters

distfloat: Inter-station distance (in km)
ref_velfloat: Reference velocity (at a given period in km/s)
tolerancefloat: Tolerance, with respect to ref_vel, used to identify the starting and ending time of the taper

Returns

tmin, tmaxfloat: Starting and ending time of the taper

Teleseismic Earthquakes (seislib.eq)

References

Download of Teleseismic-Earthquake Recordings

Inter-Station Dispersion Curves

Teleseismic Earthquakes (`seislib.eq`)