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

start()

Starts the downloads.

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

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()

Starts the downloads.

The catalog of seismic events is iterated over and, for each event, all the waveforms from the receivers associated with stations_config will be downloaded, if available. The waveforms are detrended, demeaned, removed of the instrument response, and resampled. Finally, they are saved in savedir/data/event, where ‘event’ is the origin time (in the obspy timestamp format) of the considered earthquake.

The algorithm keeps track of the progress made. This allows one to stop the downloads and start from where it was left off at any time.

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, plotting=False)

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().

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
coordsndarray of shape (n, 4)

Coordinates (lat1, lon1, lat2, lon2) of the station pairs for which a dispersion curve has been calculated

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:

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]
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:

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.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)

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.

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