Seismic Tomography (seislib.tomography)

This module allows for inverting inter-station (or epicenter-station) velocity measurements for surface-wave velocity maps. SeisLib implements a least-square inversion scheme based on ray theory [1], i.e., it assumes that the surface waves travel, from a given point on the Earth’s surface to another, without deviating from the great-circle path connecting them.

By default, to discretize the Earth’s surface, seislib.tomography._tomography.SeismicTomography implements an equal-area grid (see seislib.tomography._grid.EqualAreaGrid). Alternatively, a regular grid can be employed (see seislib.tomography._grid.EqualAreaGrid). (The latter is particularly suited to tomographic applications at local scale, where the use of equal-area parameterizations does not have clear advantages.) The grid of choice can be refined, in the areas characterized by a relatively high density of available rays, an arbitrary of times. The refinement is carried out by splitting a given block of the parameterization into four sub-blocks.

1

Boschi & Dziewonski 1999, High-and low-resolution images of the Earth’s mantle: Implications of different approaches to tomographic modeling, JGR

Parameterization

By default, SeisLib discretizes the Earth’s surface by means of equal-area grids. These prevent from artificially increasing the resolution of the resulting tomographic maps at latitudes different than zero (the effect is more prominent nearby the poles), and should therefore be preferred to Cartesian grids when investigating relatively large areas. SeisLib also allows for adaptive parameterizations, with finer resolution in the areas characterized by relatively high density of measurements. If we consider a given block intersected by more than a certain number of inter-station great-circle paths, the finer resolution is achieved by splitting it in four sub-blocks, at the midpoint along both latitude an longitude. The operation can be performed an arbitrary number of times.

class seislib.tomography._grid.EqualAreaGrid(cell_size, latmin=None, lonmin=None, latmax=None, lonmax=None, verbose=True)

Bases: seislib.tomography._grid._Grid

Class that allows for creating an equal-area grid covering the whole Earth’s surface.

Parameters

cell_sizefloat

Size of each side of the equal-area grid

latmin, lonmin, latmax, lonmaxfloat, optional

Boundaries (in degrees) of the grid

verbosebool

If True, information about the grid will be displayed. Default is True

Examples

Let’s first define an equal-area grid of \(10^{\circ} \times 10^{\circ}\). By default, this is created on the global scale.

>>> from seislib.tomography import EqualAreaGrid
>>> grid = EqualAreaGrid(cell_size=10, verbose=True)
-------------------------------------
Optimal grid found in 10 iterations
-------------------------------------
-------------------------------------
GRID PARAMETERS
Lonmin - Lonmax : -180.000 - 180.000
Latmin - Latmax : -90.000 - 90.000
Number of cells : 412
Grid cells of 10.006° : 412
-------------------------------------

>>> print(grid.mesh)
[[ 80.2098   90.     -180.      -60.    ]
[  80.2098   90.      -60.       60.    ]
[  80.2098   90.       60.      180.    ]
..., 
[ -90.      -80.2098 -180.      -60.    ]
[ -90.      -80.2098  -60.       60.    ]
[ -90.      -80.2098   60.      180.    ]]

We can now restrict the above parameterization to an arbitrary region, for example:

>>> grid.set_boundaries(latmin=0, 
...                     lonmin=0,
...                     latmax=10,
...                     lonmax=10,
...                     inplace=True)
*** GRID UPDATED ***
-------------------------------------
GRID PARAMETERS
Lonmin - Lonmax : -10.000 - 20.000
Latmin - Latmax : -10.065 - 10.065
Number of cells : 6
Grid cells of 10.006° : 6
-------------------------------------

>>> print(grid.mesh)
[[ 0.      10.0645 -10.       0.    ]
[  0.      10.0645   0.      10.    ]
[  0.      10.0645  10.      20.    ]
[-10.0645   0.     -10.       0.    ]
[-10.0645   0.       0.      10.    ]
[-10.0645   0.      10.      20.    ]]

Hint

The same result can be obtained by passing the boundaries of the region of interest directly in the initialization of the class instance, e.g.:

grid = EqualAreaGrid(cell_size=10, 
                     latmin=0, 
                     lonmin=0,
                     latmax=10,
                     lonmax=10,
                     verbose=True)

We can refine any block, say at the 0th and 1st index, simply by:

>>> grid.refine_mesh([0, 1], inplace=True)
*** GRID UPDATED ***
-------------------------------------
GRID PARAMETERS
Lonmin - Lonmax : -10.000 - 20.000
Latmin - Latmax : -10.065 - 10.065
Number of cells : 12
Grid cells of 10.006° : 4
Grid cells of 5.003° : 8
-------------------------------------

>>> print(grid.mesh)
[[ 5.03225  10.0645  -10.       -5.     ]
[  5.03225  10.0645   -5.        0.     ]
[  0.        5.03225 -10.       -5.     ]
[  0.        5.03225  -5.        0.     ]
[  5.03225  10.0645    0.        5.     ]
[  5.03225  10.0645    5.       10.     ]
[  0.        5.03225   0.        5.     ]
[  0.        5.03225   5.       10.     ]
[  0.       10.0645   10.       20.     ]
[-10.0645    0.      -10.        0.     ]
[-10.0645    0.        0.       10.     ]
[-10.0645    0.       10.       20.     ]]

Note that the size of the two blocks defined at the first two rows of the above have been halved.

To display the grid, use the plot() method >>> grid.plot(projection=’Mercator’)

Attributes

verbosebool

If True, information about the grid will be displayed.

refinedint

Number of times the grid has been “refined”

lon_spanndarray of shape (n,)

Longitudinal span of each block in the n latitudinal bands defining the grid

meshndarray of shape (m, 4)

Blocks bounded by parallel1, parallel2, meridian1, meridian2

latmin, lonmin, latmax, lonmaxfloat

Minimum and maximum latitudes and longitudes of the blocks defining the grid

Methods

best_grid_parameters(cell_side)	Finds the spatial parameterization that most closely approximates the size of the grid cells required by the user.
global_mesh(ncells, lon_span)	Builds an equal-area global mesh given the number of cells and longitude span as a function of latitude.
grid_parameters(nrings)	Computes the grid parameters (number of cells, cells area, cells side, longitude span as a function of latitude) [1].
index_lon_lat(lon, lat)	Returns the mesh index corresponding with the coordinates (lat, lon)
indexes_in_polygon(poly)	Returns the mesh indexes whose midpoints (lon, lat) fall outside the specified rectangular region.
indexes_in_region(latmin, latmax, lonmin, lonmax)	Returns the mesh indexes whose midpoints (lon, lat) fall inside the specified rectangular region.
indexes_out_region(latmin, latmax, lonmin, ...)	Returns the mesh indexes whose midpoints (lon, lat) fall outside the specified rectangular region.
midpoints_lon_lat()	Returns the midpoints (lon, lat) of each grid's block
parallels_first_pixel([mesh])	Generator function yielding the indexes at which a change in the parallel coordinates is found.
plot([ax, mesh, projection, meridian_min, ...])	Plots the (adaptive) equal-area grid
refine_mesh(ipixels[, mesh, inplace])	Halves the size of the specified pixels
select_cells(indexes[, inplace])	Mesh indexing
set_boundaries(latmin, latmax, lonmin, lonmax)	Restricts the mesh to the required boundaries
split_pixel(newmesh, i_newpixel, parallel1, ...)	Called by refine_mesh() to create new pixels
update_grid_params(mesh[, refined])	Updates the grid boundaries stored in the EqualAreaGrid instance and the number of times that the grid has been refined.

best_grid_parameters(cell_side)

Finds the spatial parameterization that most closely approximates the size of the grid cells required by the user. It exploits a 1-D grid search

Parameters

cell_sidefloat

Side’s size of the desidered grid cells

Returns

grid parameters associated with the best parameterization. (See grid_parameters())

global_mesh(ncells, lon_span)

Builds an equal-area global mesh given the number of cells and longitude span as a function of latitude. [1]

Parameters

ncellsint

Number of grid cells

lon_spanndarray

Longitudinal span the grid cells in each latitudinal ring

Returns

gridndarray (n, 4)

Array containing n pixels bounded by parallel1, parallel2, meridian1, meridian2. The grid is rounded to the 4rd decimal digit, for improved numerical stability in most applications

References

1

Malkin 2019, A new equal-area isolatitudinal grid on a spherical surface, The Astronomical Journal

grid_parameters(nrings)

Computes the grid parameters (number of cells, cells area, cells side, longitude span as a function of latitude) [1].

Parameters

nringsint

Number of latitudinal rings used to subdivide the Earth.

Returns

ncellsint

Number of grid cells

cell_sidefloat

Latitudinal extent of the blocks in degrees, (corresponding to the sqrt of the area)

lon_spanndarray

Longitudinal span the grid cells in each latitudinal band

References

1(1,2)

Malkin 2019, A new equal-area isolatitudinal grid on a spherical surface, The Astronomical Journal

parallels_first_pixel(mesh=None)

Generator function yielding the indexes at which a change in the parallel coordinates is found.

Parameters

meshndarray (n, 4)

Array containing n pixels bounded by parallel1, parallel2, meridian1, meridian2. If None, the mesh stored in the EqualAreaGrid instance (self.mesh) is used. Default is None.

Yields

iint

class seislib.tomography._grid.RegularGrid(cell_size, latmin, lonmin, latmax, lonmax, verbose=True)

Bases: seislib.tomography._grid._Grid

Class that allows allows for creating a regular grid in the format required by seislib.tomography.tomography.SeismicTomography. This class is particularly suited to tomographic applications at local scale, where the use of equal-area parameterizations does not have clear advantages.

Parameters

cell_sizefloat, (2,) tuple

Size of each side of the regular grid. If a (2,) tuple is passed, this will be interpreted as (dlon, dlat), where dlon and dlat are the longitudinal and latitudinal steps characterizing the grid

latmin, lonmin, latmax, lonmaxfloat

Boundaries (in degrees) of the grid

verbosebool

If True, information about the grid will be displayed. Default is True

Examples

Let’s first define a regular grid of \(0.1^{\circ} \times 0.1^{\circ}\). We will restrict the study area to \(9 \leq lon \leq 11\) and \(40 \leq lat \leq 42\).

>>> from seislib.tomography import RegularGrid
>>> grid = RegularGrid(cell_size=0.1, 
                       latmin=40,
                       latmax=42,
                       lonmin=9,
                       lonmax=11,
                       verbose=True)
-------------------------------------
GRID PARAMETERS
Lonmin - Lonmax : 10.000 - 11.000
Latmin - Latmax : 40.000 - 41.000
Number of cells : 100
Grid cells of 0.100° : 100
-------------------------------------

>>> print(grid.mesh)
[[41.9 42.  10.9 11. ]
 [41.9 42.  10.8 10.9]
 [41.9 42.  10.7 10.8]
 ...
 [40.  40.1  9.2  9.3]
 [40.  40.1  9.1  9.2]
 [40.  40.1  9.   9.1]]

We can refine any block, say at the 0th and 1st index, simply by:

>>> grid.refine_mesh([0, 1], inplace=True)
*** GRID UPDATED ***
-------------------------------------
GRID PARAMETERS
Lonmin - Lonmax : 9.000 - 11.000
Latmin - Latmax : 40.000 - 42.000
Number of cells : 406
Grid cells of 0.100° : 398
Grid cells of 0.050° : 8
-------------------------------------

>>> print(grid.mesh)
[[41.95 42.   10.9  10.95]
 [41.95 42.   10.95 11.  ]
 [41.9  41.95 10.9  10.95]
 ...
 [40.   40.1   9.2   9.3 ]
 [40.   40.1   9.1   9.2 ]
 [40.   40.1   9.    9.1 ]]

Note that the size of the two blocks defined at the first two rows of the above have been halved.

To display the grid, use the plot() method >>> grid.plot(projection=’Mercator’)

Attributes

verbosebool

If True, information about the grid will be displayed.

refinedint

Number of times the grid has been “refined”

meshndarray of shape (m, 4)

Blocks bounded by parallel1, parallel2, meridian1, meridian2

latmin, lonmin, latmax, lonmaxfloat

Minimum and maximum latitudes and longitudes of the blocks defining the grid

Methods

create_mesh(dlon, dlat[, latmin, latmax, ...])	Creates the regular grid.
index_lon_lat(lon, lat)	Returns the mesh index corresponding with the coordinates (lat, lon)
indexes_in_polygon(poly)	Returns the mesh indexes whose midpoints (lon, lat) fall outside the specified rectangular region.
indexes_in_region(latmin, latmax, lonmin, lonmax)	Returns the mesh indexes whose midpoints (lon, lat) fall inside the specified rectangular region.
indexes_out_region(latmin, latmax, lonmin, ...)	Returns the mesh indexes whose midpoints (lon, lat) fall outside the specified rectangular region.
midpoints_lon_lat()	Returns the midpoints (lon, lat) of each grid's block
plot([ax, mesh, projection, meridian_min, ...])	Plots the (adaptive) equal-area grid
refine_mesh(ipixels[, mesh, inplace])	Halves the size of the specified pixels
select_cells(indexes[, inplace])	Mesh indexing
set_boundaries(latmin, latmax, lonmin, lonmax)	Restricts the mesh to the required boundaries
split_pixel(newmesh, i_newpixel, parallel1, ...)	Called by refine_mesh() to create new pixels
update_grid_params(mesh[, refined])	Updates the grid boundaries stored in the EqualAreaGrid instance and the number of times that the grid has been refined.

create_mesh(dlon, dlat, latmin=None, latmax=None, lonmin=None, lonmax=None)

Creates the regular grid.

Parameters

dlon, dlatfloat

Longitudinal and latitudinal steps characterizing the grid

latmin, lonmin, latmax, lonmaxfloat

Boundaries (in degrees) of the grid

Returns

gridndarray (n, 4)

Array containing n pixels bounded by parallel1, parallel2, meridian1, meridian2. The grid is rounded to the 4rd decimal digit, for improved numerical stability in most applications

Least-Squares Imaging (Ray Theory)

To map lateral variations in surface-wave velocity, SeisLib implements a least-square inversion scheme based on ray theory. The travel time along the great-circle path can be written \(t = \int_{path}{s(\phi(l), \theta(l)) dl}\), where \(\phi\) and \(\theta\) denote longitude and latitude, and s the sought Earth’s slowness.

Let us consider a discrete parameterization of the Earth’s surface, and assume each block (or grid cell) of such parameterization has constant slowness. The above integral expression can then be reformulated in the discrete form

\[s = \frac{1}{L} \sum_{n}{s_n l_n},\]

where L is the length of the great-circle path and l the distance traveled by the surface wave through the \(n\)th block. The above equation represents the forward calculation that allows for retrieving the average velocity of propagation between two points on the Earth’s surface, provided that the (discrete) spatial variations in velocity (or slowness) are known. If we now define the \(m \times n\) matrix such that \(A_{ij} = \frac{l_j}{L_i}\), where \(L_i\) is the length of the great circle associated with \(i\)th observation, we can switch to matrix notation and write

\[\bf A \cdot x = d,\]

where \(\bf d\) is an m-vector whose \(k\)th element corresponds to the measured slowness, and \(\bf x\) the sought n-vector whose \(k\)th element corresponds to the model coefficient \(s_k\). Matrix \(\bf A\) can be computed numerically in a relatively simple fashion. In real-world seismological applications, however, the above system of equations is often strongly overdetermined, i.e. the number of data points is much larger than the number of model parameters (\(m \gg n\)). This implies that, although \(\bf A\) is known, it is not invertible.

SeisLib solves the above inverse problem in a regularized least-squares sense [1], i.e.

\[{\bf x} = {\bf x}_0 + \left( {\bf A}^T \cdot {\bf A} + \mu^2 {\bf R}^T \cdot {\bf R} \right)^{-1} \cdot {\bf A}^T \cdot ({\bf d - A} \cdot {\bf x}_0),\]

where the roughness operator \(\bf R\) is dependent on the parameterization, \({\bf x}_0\) is a reference model, and the scalar weight \(\mu\) should be chosen via L-curve analysis [2].

References

1

Aster et al. 2018, Parameter estimation and inverse problems

2

Hansen 2001, The L-curve and its use in the numerical treatment of inverse problems

class seislib.tomography._tomography.SeismicTomography(cell_size=4, latmin=None, lonmin=None, latmax=None, lonmax=None, verbose=True, regular_grid=False)

Bases: object

Class to obtain velocity maps using a linearized inversion based on the ray theory (infinite frequency approximation and wave propagation along the great-circle path connecting two points on the Earth’s surface).

Parameters

cell_sizeint

Size of each side of the equal-area grid

latmin, lonmin, latmax, lonmaxfloat, optional

Boundaries (in degrees) of the grid

regular_gridbool

If False (default), the study area is discretized using an equal-area parameterization. Otherwise, a regular grid is employed.

verbosebool

If True, information about the grid and data will be displayed in console. (Default is True)

Examples

The following will calculate a phase-velocity map at a given period, say 10 s, based on inter-station measurements of surface-wave velocity. In practice, our data consist of a ndarray of shape (n, 5), where n is the number of inter-station measurements of phase velocity (extracted, for example, via seislib.an.an_velocity.AmbientNoiseVelocity), and the five columns consist of lat1 (°), lon1 (°), lat2 (°), lon2 (°), velocity (m/s), respectively. (-180<=lon<180, -90<=lat<90). This matrix has been saved to /path/to/data.txt

We will discretize the study area using an equal-area parameterization, characterized by blocks of :math:`2^{circ} times 2^{circ}. We will then iteratively refine the parameterization up to a maximum number of 2 times, to reach a maximum resolution of 0.5° in the areas of the map characterized by a relatively large number of measurements. (This refinement can be carried out an arbitrary number of times.)

First, we need to initialize the SeismicTomography instance and load our data into memory:

>>> from seislib.tomography import SeismicTomography
>>> tomo = SeismicTomography(cell_size=2, regular_grid=False, verbose=True)
>>> tomo.add_data(src='/path/to/data')
-------------------------------------
Optimal grid found in 46 iterations
-------------------------------------
-------------------------------------
GRID PARAMETERS
Lonmin - Lonmax : -180.000 - 180.000
Latmin - Latmax : -90.000 - 90.000
Number of cells : 10312
Grid cells of 2.000° : 10312
-------------------------------------
DATA PARAMETERS
Lonmin - Lonmax data : -124.178 - -69.292
Latmin - Latmax data : 23.686 - 48.470
Number of measurements : 1064
Source : /path/to/data
-------------------------------------        

Hint

we could have directly passed the data matrix without loading it from disk. If data is your ndarray of shape (m, 5), you can pass it to tomo by:

tomo.add_data(data=your_matrix)

Hint

you can add, sequentially, how many data sets you wish. The data information will be automatically updated

Hint

to display general information on the data, type print(tomo). To display general information on the parameterization, type print(tomo.grid). (See also the documentation on seislib.tomography.grid.EqualAreaGrid)

Hint

if you are interested in working at local scale, where the use of equal-area parameterizations does not have clear advantages, consider setting regular_grid=True when initializing the instance of SeismicTomography, together with the boundaries of the region you are interested in (lonmin, lonmax, latmin, latmax). For example:

from seislib.tomography import SeismicTomography

tomo = SeismicTomography(cell_size=0.05, 
                         regular_grid=True, 
                         latmin=40,
                         latmax=41,
                         lonmin=10,
                         lonmax=12,
                         verbose=True)

Now we can restrict the boundaries of the (global) equal-area parameterization to the minimum and maximum latitude and longitude spanned by our data:

>>> tomo.grid.set_boundaries(latmin=tomo.latmin_data,
                             latmax=tomo.latmax_data,
                             lonmin=tomo.lonmin_data,
                             lonmax=tomo.lonmax_data)
*** GRID UPDATED ***
-------------------------------------
GRID PARAMETERS
Lonmin - Lonmax : -126.234 - -66.614
Latmin - Latmax : 22.006 - 50.005
Number of cells : 324
Grid cells of 2.000° : 324
-------------------------------------

Having done so, everything is ready to calculate the coefficients of the \(\bf A\) matrix (i.e., of the data kernel), and to refine the parameterization up to two times in the areas characterized by a relatively high density of measurements. We will define such regions (i.e., model parameters) as those where there are at least 150 inter-station great-circle paths intersecting them. In doing so, we will remove the grid cells of 2° that are not intersected by at least one great-circle path (see the argument keep_empty_cells).

>>> tomo.compile_coefficients(keep_empty_cells=False)
>>> tomo.refine_parameterization(hitcounts=150, 
...                              keep_empty_cells=True)    
>>> tomo.refine_parameterization(hitcounts=150, 
...                              keep_empty_cells=True)    
*** GRID UPDATED ***
-------------------------------------
GRID PARAMETERS
Lonmin - Lonmax : -126.142 - -66.614
Latmin - Latmax : 22.006 - 50.005
Number of cells : 202
Grid cells of 2.000° : 202
-------------------------------------
*** GRID UPDATED ***
-------------------------------------
GRID PARAMETERS
Lonmin - Lonmax : -126.142 - -66.614
Latmin - Latmax : 22.006 - 50.005
Number of cells : 808
Grid cells of 2.000° : 0
Grid cells of 1.000° : 808
-------------------------------------
*** GRID UPDATED ***
-------------------------------------
GRID PARAMETERS
Lonmin - Lonmax : -126.142 - -66.614
Latmin - Latmax : 22.006 - 50.005
Number of cells : 2284
Grid cells of 2.000° : 0
Grid cells of 1.000° : 316
Grid cells of 0.500° : 1968
-------------------------------------

To obtain the velocity map, we now need to carry out the inversion. We will apply a roughness-damping regularization, using a roughness coefficient equal to 3e-3 (to select a proper roughness damping, check lcurve()). We then plot the retrieved velocity:

c = 1 / tomo.solve(rdamp=3e-3)
tomo.plot_map(tomo.grid.mesh, c)

(Note that the solve method returns slowness, hence we took the inverse of the solution.) A checkerboard test, to verify the ability of our data to resolve the lateral differences in the study area, can simply be performed and visualized by:

restest = tomo.checkerboard_test(kx=6,
                                 ky=6,
                                 latmin=tomo.grid.latmin,
                                 latmax=tomo.grid.latmax,
                                 lonmin=tomo.grid.lonmin,
                                 lonmax=tomo.grid.lonmax,
                                 cell_size=0.5,
                                 anom_amp=0.1,
                                 noise=0,
                                 rdamp=3e-3)

input_model = 1 / restest['synth_model']
input_mesh = restest['mesh']
retrieved_model = 1 / restest['retrieved_model']

tomo.plot_map(tomo.grid.mesh, retrieved_model)
tomo.plot_map(input_mesh, input_model)

Note

The resolution test (check also resolution_test() and spike_test()) should be performed after (and only after) the coefficients in the data-kernel (tomo.A matrix) have been complied, i.e., after having called compile_coefficients() (and eventually refine_parameterization(), which updates tomo.A).

Finally, we show how we can obtain a heat map of the raypaths intersecting each model parameter:

raypaths = tomo.raypaths_per_pixel()
img, cb = tomo.plot_map(mesh=tomo.grid.mesh, 
                        c=raypaths, 
                        cmap='cividis', 
                        style='contourf', 
                        levels=20,
                        show=False)
cb.set_label('Raycounts')

Hint

When the above argument show is False, the image and colorbar will be returned. That allow us to change the label on the colorbar.

Hint

For more control on the plots, consider creating your own instance of GeoAxesSubplot (see cartopy documentation) before calling the plot_map function (or, alternatively, colormesh(), contour(), or contourf() of SeismicTomography). For example:

from seislib.plotting import make_colorbar

fig = plt.figure(figsize=(12, 12))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Mercator())
ax.coastlines(resolution='50m', color='k', lw=1, zorder=100)
img = tomo.colormesh(mesh=tomo.grid.mesh, 
                    c=raypaths, 
                    ax=ax, 
                    cmap='cividis', 
                    shading='flat', 
                    edgecolors='face')
map_boundaries = (tomo.grid.lonmin, tomo.grid.lonmax, 
                tomo.grid.latmin, tomo.grid.latmax)
ax.set_extent(map_boundaries, ccrs.PlateCarree())  
cb = make_colorbar(ax, img, orientation='horizontal')
cb.set_label(label='Raycounts', labelpad=10, fontsize=22)

Attributes

verbosebool

If True (default), information about the grid and data will be displayed in console

gridseislib.EqualAreaGrid

Instance of seislib.tomography.grid.EqualAreaGrid, corresponding to an equal-area parameterization

lonmin_data, latmin_data, lonmax_data, latmax_datafloat

Minimum and maximum values of latitude and longitude in the data. Only available after function call add_data()

data_coordsndarray of shape (n, 4)

Lat1, lon1, lat2, lon2 of the pairs of receivers (or epicenter-receiver), in degrees (-180<lon<180, -90<lat<90). Only available after function call add_data()

velocityndarray of shape (n,)

Velocity (in m/s) measured between each pair of stations or epicenter-station. Only available after function call add_data()

refvelfloat

Reference velocity (in m/s), used in the inversion. Only available after function call add_data()

Andarray of shape (n, m)

Jacobian matrix. Where n is the number of data, m is the number of parameters (i.e., number of rows in self.grid.mesh). Only available after function call compile_coefficients()

Rndarray of shape (m, m)

Roughness operator, where m is the number of parameters (i.e., number of rows in self.grid.mesh). Only available after function call solve()

Lndarray of shape (m, m)

Norm damping operator, where m is the number of parameters (i.e., number of rows in self.grid.mesh). Only available after function call solve()

Methods

add_data([data, refvel, src])	Loads the data to be used for tomography.
azimuth_backazimuth(lat1, lon1, lat2, lon2)	Calculates the azimuth and backazimuth (in degrees) between coordinate points (in degrees).
checkerboard(ref_value[, kx, ky, latmin, ...])	Checkerboard-like pattern on Earth
checkerboard_test(kx, ky[, regular_grid, ...])	Resolution test, known as "checkerboard test".
colormesh(mesh, c, ax, **kwargs)	Adaptation of matplotlib.pyplot.pcolormesh to the (adaptive) equal-area grid.
compile_coefficients([refine, coeff_matrix, ...])	Compiles the matrix of coefficients (A) used to perform the inversion.
contour(mesh, c, ax[, smoothing])	Adaptation of matplotlib.pyplot.contour to the (adaptive) equal area grid
contourf(mesh, c, ax[, smoothing])	Adaptation of matplotlib.pyplot.contourf to the (adaptive) equal-area grid
delay_to_velocity(lat1, lon1, lat2, lon2, ...)	Converts a delay to velocity
gc_distance(lat1, lon1, lat2, lon2)	Calculates the great circle distance (in m) between coordinate points (in degrees).
lcurve([A, slowness, refvel, mesh, damping, ...])	L-curve analysis.
load(path)	Loads the pickle file at the specified path and returns the associated instance of SeismicTomography.
norm_operator(mesh[, return_sparse])	Computes the norm damping operator associated with the parameterization
plot_colored_rays([ax, show, cmap, vmin, ...])	Utility function to display the great-circle paths associated with pairs of data coordinates, colored according to their respective measurements
plot_map(mesh, c[, ax, projection, ...])	Utility function to display the lateral variations of some quantity on a equal-area grid
plot_rays([ax, show, stations_color, ...])	Utility function to display the great-circle paths associated with pairs of data coordinates
plot_stations([ax, show, oceans_color, ...])	Creates a maps of seismic receivers
raypaths_per_pixel()	Calculates the number of raypaths in each pixel of the mesh stored in the SeismicTomography instance (self.grid.mesh)
reduce_measurements([latmin, latmax, ...])	Any measurement not intersecting the specified region is removed from the data.
refine_parameterization([hitcounts, ...])	Halves the dimension of the pixels intersected by a number of raypaths >= hitcounts.
resolution_test(mesh, velocity_model[, ...])	Resolution test on an arbitrary input model
roughness_operator(mesh[, return_sparse])	Computes the roughness operator associated with the parameterization
save(path)	Saves the SeismicTomography instance to the specified path in the pickle compressed format.
solve([A, slowness, refvel, mesh, ndamp, rdamp])	Method for solving the (regularized) inverse problem Ax = b, where A is the matrix of coefficients, x the sought model of slowness, and b the array of measurements expressed in slowness.
spike(ref_value, x0, y0[, sigma_x, sigma_y, ...])	Spike-like pattern on Earth, created via a 3-D Gaussian
spike_test(x0, y0, sigma_x, sigma_y[, ...])	Resolution test, known as "spike test".
update_data_info([refvel])	Updates the information on minimum / maximum longitude and latitude of the data and the reference velocity
velocity_to_delay(lat1, lon1, lat2, lon2, ...)	Converts a velocity to delay

add_data(data=None, refvel=None, src=None, **kwargs)

Loads the data to be used for tomography.

If the argument data isn’t passed in, the data are loaded from src. If the instance of SeismicTomography already contains data, the new data will be added (concatenated) to the previously loaded ones.

Note

In the current implementation, the 6th column of the data (i.e., the standard deviation), even if present, is not employed in the inversion (see solve()). This behaviour will be changed in a future release.

Parameters

datandarray, optional

The array must be of shape (n, 5) or (n, 6). Cols 1 to 4 must contain the coordinates of the measurements (lat1, lon1, lat2, lon2), col 5 the velocity, and col 6 the standard deviation (optional). (-180<lon<180, -90<lat<90, velocity in m/s)

refvelfloat, optional

Reference velocity (in m/s). If None, the average velocity is used as reference

srcstr, optional

If data is None, src is used to load the data. The file extension must be .txt

**kwargsoptional

Arguments to be passed to numpy.loadtxt

Returns

None

The data are stored in self.__dict__ and can be accessed through self.data_coords, self.velocity, self.refvel, self.std (if present), self.latmin_data, self.lonmin_data, self.latmax_data, self.lonmax_data

classmethod azimuth_backazimuth(lat1, lon1, lat2, lon2)

Calculates the azimuth and backazimuth (in degrees) between coordinate points (in degrees). This function calls directly the obspy gps2dist_azimuth, it only extends its functionality through the numpy.vectorize decorator

Parameters

lat1, lon1, lat2, lon2float or array-like of shape (n,)

Coordinates of the points on the Earth’s surface, in degrees.

Returns

tuple of shape (2,) or ndarray of shape (n, 2)

The columns consist of azimuth azimuth and backazimuth

static checkerboard(ref_value, kx=10, ky=10, latmin=None, latmax=None, lonmin=None, lonmax=None, anom_amp=0.1)

Checkerboard-like pattern on Earth

Parameters

ref_valuefloat

Reference value with respect to which the anomalies are generated

kx, kyfloat

Frequency of anomalies along longitude (kx) and latitude (ky) within the boundaries defined by lonmin, lonmax, latmin, latmax. Defaults are 10.

latmin, latmax, lonmin, lonmaxfloat, optional

Boundaries of the checkerboard (in degrees).

anom_ampfloat

Amplitude of the anomalies. Default is 0.1 (i.e. 10% of the reference value)

Returns

Function

Generating function that should be called passing the values of longitude and latitude where the values of the checkerboard want to be retrieved.

Examples

>>> from seislib.tomography import SeismicTomography
>>> tomo = SeismicTomography(5)
>>> lat = np.mean(tomo.grid.mesh[:, :2], axis=1)
>>> lon = np.mean(tomo.grid.mesh[:, 2:], axis=1)
>>> pattern = tomo.checkerboard(ref_value=3, kx=6, ky=6)
>>> c = pattern(lon, lat)
>>> img, cb = tomo.plot_map(tomo.grid.mesh, c, projection='Robinson', 
...                         show=False)
>>> cb.set_label('Checkerboard Pattern')

checkerboard_test(kx, ky, regular_grid=False, latmin=None, latmax=None, lonmin=None, lonmax=None, cell_size=5, anom_amp=0.1, refvel=None, noise=0, ndamp=0, rdamp=0)

Resolution test, known as “checkerboard test”. The method first builds synthetic data (velocities) using a checkerboard() pattern as input (velocity) model.

Parameters

kx, kyint , float

Frequency of anomalies along longitude (kx) and latitude (ky) within the boundaries defined by lonmin, lonmax, latmin, latmax

regular_gridbool

If False (default), the study area is discretized using an equal-area parameterization. Otherwise, a regular grid is employed.

latmin, latmax, lonmin, lonmaxfloat, optional

Boundaries of the checkerboard, in degrees. (-180<lon<180, -90<lat<90)

cell_sizeint | float

Size of the checkerboard grid-cell sides (in degrees). Larger values are preferrable to avoid the resolution of an overly large inverse problem

anom_ampfloat

Intensity of the anomalies with respect to ref_value. A value of 0.1 means 10% of ref_value

refvelfloat, optional

Original values of the anomalies (black and white squares of a checkerboard). By default, their value is +-1. Therefore, if anom_extent is 0.1 (say) their final values will be +-0.1

noisefloat, optional

If > 0, random noise expressed as percentage of the synthetic data

ndampfloat

Norm damping coefficients (default is zero)

rdampfloat

Roughness damping coefficients (default is zero)

Returns

dict

The dictionary is structured as follows:

• ‘synth_data’: ndarray of shape (m,), where m corresponds to the number of rows in the A matrix (i.e., in the data kernel)

• ‘synth_model’: inverse of velocity_model, used to create the synthetic data

• ‘retrieved_model’: model retrieved from the inversion of the synthetic data

• ‘mesh’: same as the input mesh

classmethod colormesh(mesh, c, ax, **kwargs)

Adaptation of matplotlib.pyplot.pcolormesh to the (adaptive) equal-area grid.

Parameters

meshndarray (n, 4)

Equal area grid, consisting of n pixels described by lat1, lat2, lon1, lon2

clist of ndarray (n,)

Values to plot in each grid cell

axcartopy.mpl.geoaxes.GeoAxesSubplot

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

**kwargs

Additional inputs passed to seislib.plotting.colormesh

Returns

imgInstance of matplotlib.collections.QuadMesh

compile_coefficients(refine=False, coeff_matrix=None, keep_empty_cells=True)

Compiles the matrix of coefficients (A) used to perform the inversion. It calls the _compile_coefficients function, written in cython.

Parameters

refinebool

should be True if the user has refined the parameterization of the mesh (by calling refine_parameterization). This allows to reduce computation time. (Default is False)

coeff_matrixif None, a new A will be created and stored in the

SeismicTomography instance (self.A). A should be passed in after the parameterization refinement. (Default is None)

keep_empty_cellsbool

If False, cells intersected by no raypaths are removed

Returns

None

The SeismicTomography instance will be updated, so that A can be accessed typing self.A

classmethod contour(mesh, c, ax, smoothing=None, **kwargs)

Adaptation of matplotlib.pyplot.contour to the (adaptive) equal area grid

Parameters

meshndarray (n, 4)

Equal area grid, consisting of n pixels described by lat1, lat2, lon1, lon2

clist of ndarray (n,)

Values to plot in each grid cell

axcartopy.mpl.geoaxes.GeoAxesSubplot

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

smoothingfloat, optional

Value passed to scipy.ndimage.filters.gaussian_filter and used to obtain a smooth representation of c (default is None)

**kwargs

Additional inputs passed to seislib.plotting.contourf

Returns

imgInstance of matplotlib.contour.QuadContourSet

classmethod contourf(mesh, c, ax, smoothing=None, **kwargs)

Adaptation of matplotlib.pyplot.contourf to the (adaptive) equal-area grid

Parameters

meshndarray (n, 4)

Equal area grid, consisting of n pixels described by lat1, lat2, lon1, lon2

clist of ndarray (n,)

Values to plot in each grid cell

axcartopy.mpl.geoaxes.GeoAxesSubplot

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

smoothingfloat, optional

Value passed to scipy.ndimage.filters.gaussian_filter and used to obtain a smooth representation of c (default is None)

**kwargs

Additional inputs passed to seislib.plotting.contourf

Returns

imgInstance of matplotlib.contour.QuadContourSet

classmethod delay_to_velocity(lat1, lon1, lat2, lon2, delay, refvel)

Converts a delay to velocity

A negative delay dt corresponds to faster velocity with respect to the reference velocity v0:

dt = x/v - x/v0
v = x / (dt + x/v0)

where x denotes distance and v the observed velocity.

Parameters

lat1, lon1, lat2, lon2float or array-like of shape (n,)

Coordinates (in radians) associated with the dalay

delayfloat or array-like of shape (n,)

Measured delay (in seconds)

refvelfloat

Reference velocity used to calculate the delay (in m/s)

Returns

float or array-like of shape (n,)

classmethod gc_distance(lat1, lon1, lat2, lon2)

Calculates the great circle distance (in m) between coordinate points (in degrees). This function calls directly the obspy gps2dist_azimuth, it only extends its functionality through the numpy.vectorize decorator.

Parameters

lat1, lon1, lat2, lon2float or array-like of shape (n,)

Coordinates of the points on the Earth’s surface, in degrees.

Returns

Great-circle distance (in m)

If the input is an array (or list) of coordinates, an array of distances is returned

lcurve(A=None, slowness=None, refvel=None, mesh=None, damping='roughness', n=20, damp_min=1e-05, damp_max=0.1, logspace=True, show=True)

L-curve analysis. The function calls iteratively the solve method

Parameters

Andarray, shape (m, n), optional

Matrix of the coefficients. If None (default), the A matrix stored in the SeismicTomography instance is used

slownessndarray, shape (m,), optional

Inverse of the velocity measurements. Corresponds to b in the equation Ax = b. If None (default), the velocity stored in the SeismicTomography instance is used to compute slowness

refvelfloat, optional

Reference velocity (in m/s) used for the initial guess of the velocity model. If provided, should stabilize the inversion and make it converge faster. If None (default), the reference velocity stored in the SeismicTomography instance is used

meshndarray, shape (n, 4), optional

Grid cells in seislib format, consisting of n pixels described by lat1, lat2, lon1, lon2. If None (default), the mesh stored in the SeismicTomography instance is used

damping{‘norm’, ‘roughness’}

Damping criterion used in the inversion. If ‘norm’, a norm damping is applied, otherwise the inversion is regularized using roughness damping (default)

nint

Number of inversions performed in the analysis

damp_minfloat

Minimum damping term (default is 1e-5)

damp_maxfloat

Maximum damping term (default is 1e-1)

logspacebool

If True (default), the damping range is defined to be linearly spaced in logarithmic scale between damp_min and damp_max, i.e., numpy.logspace(np.log10(damp_min), np.log10(damp_ax), n). False corresponds to numpy.linspace(damp_min, damp_max, n)

showbool

If True (default), the result of the analysis is displayed

Returns

damp_rangendarray (n,)

Damping range used in the analysis

resulttuple of shape (2,)

(i) Residual norm |Ax - b| and (ii) norm of the damped solution |Gx|, where G = I (identity matrix) if the L-curve was performed using a norm damping (ndamp>0), or the roughness operator if the roughness damping was used instead (rdamp>0)

classmethod load(path)

Loads the pickle file at the specified path and returns the associated instance of SeismicTomography.

Parameters

pathstr

Absolute path of the file

static norm_operator(mesh, return_sparse=True)

Computes the norm damping operator associated with the parameterization

Parameters

meshndarray, shape (n, 4)

Grid cells in seislib format, consisting of n pixels described by lat1, lat2, lon1, lon2 (in degrees). If None (default), the mesh stored in the SeismicTomography instance is used

return_sparsebool

If True, the matrix of roughness is converted to scipy.sparse.csr_matrix

Returns

ndarray of shape (n,n)

Norm damping operator

plot_colored_rays(ax=None, show=True, cmap=<matplotlib.colors.ListedColormap object>, vmin=None, vmax=None, stations_color='k', oceans_color='lightgrey', lands_color='w', edgecolor='k', stations_alpha=None, paths_alpha=None, projection='Mercator', resolution='110m', map_boundaries=None, bound_map=True, paths_width=1, colorbar=True, cbar_dict={}, **kwargs)

Utility function to display the great-circle paths associated with pairs of data coordinates, colored according to their respective measurements

Parameters

axcartopy.mpl.geoaxes.GeoAxesSubplot

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

showbool

If True (default), the map will be showed once generated. Otherwise a GeoAxesSubplot instance is returned

cmapstr or Colormap

If str, it should be a valid matplotlib.cm.colormap name (see matplotlib documentation).

vmin, vmaxfloat

Boundaries of the colormap. If None, the minimum and maximum values of c will be taken

stations_colorstr

Color of the receivers and of the great-circle paths (see matplotlib documentation for valid color names). Defaults are ‘r’ (red) and ‘k’ (black)

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)

stations_alpha, paths_alphafloat, optional

Transparency of the stations and of the great-circle paths. Defaults are None and 0.3

paths_widthfloat

Linewidth of the great-circle paths

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.

map_boundarieslist or tuple of floats, shape (4,), optional

Lonmin, lonmax, latmin, latmax (in degrees) defining the extent of the map

bound_mapbool

If True, the map boundaries will be automatically determined. Ignored if map_boundaries is not None

colorbarbool

If True (default), a colorbar associated with c is displayed on the side of the map

cbar_dictdict

Keyword arguments passed to matplotlib.colorbar.ColorbarBase

**kwargs

Additional keyword arguments passed to ax.plot

Returns

If show is True

None

Otherwise

GeoAxesSubplot instance together with an instance of matplotlib.colorbar.Colorbar (if colorbar is True)

classmethod plot_map(mesh, c, ax=None, projection='Mercator', map_boundaries=None, bound_map=True, colorbar=True, show=True, style='colormesh', add_features=False, resolution='110m', cbar_dict={}, **kwargs)

Utility function to display the lateral variations of some quantity on a equal-area grid

Parameters

meshndarray, shape (n, 4)

Grid cells in seislib format, consisting of n pixels described by lat1, lat2, lon1, lon2 (in degrees)

cndarray, shape (n,)

Values associated with the grid cells (mesh)

axcartopy.mpl.geoaxes.GeoAxesSubplot

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

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’

map_boundarieslist or tuple of floats, shape (4,), optional

Lonmin, lonmax, latmin, latmax (in degrees) defining the extent of the map

bound_mapbool

If True, the map boundaries will be automatically determined. Ignored if map_boundaries is not None

colorbarbool

If True (default), a colorbar associated with c is displayed on the side of the map

showbool

If True (default), the map will be showed once generated

style{‘colormesh’, ‘contourf’, ‘contour’}

Possible options are ‘colormesh’, ‘contourf’, and ‘contour’, each corresponding to the homonymous method. Default is ‘colormesh’

add_featuresbool

If True, natural Earth features will be added to the GeoAxesSubplot. Default is False. If ax is None, it is automatically set to True

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’

cbar_dictdict

Keyword arguments passed to matplotlib.colorbar.ColorbarBase

**kwargs

Additional inputs passed to the ‘colormesh’, ‘contourf’, and ‘contour’ methods of SeismicTomography

Returns

If show is True

None

Otherwise

an instance of matplotlib.collections.QuadMesh, together with an instance of matplotlib.colorbar.Colorbar (if colorbar is True)

plot_rays(ax=None, show=True, stations_color='r', paths_color='k', oceans_color='water', lands_color='land', edgecolor='k', stations_alpha=None, paths_alpha=0.3, projection='Mercator', resolution='110m', map_boundaries=None, bound_map=True, paths_width=0.2, **kwargs)

Utility function to display the great-circle paths associated with pairs of data coordinates

Parameters

axcartopy.mpl.geoaxes.GeoAxesSubplot

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

showbool

If True (default), the map will be showed once generated. Otherwise a GeoAxesSubplot instance is returned

stations_color, paths_colorstr

Color of the receivers and of the great-circle paths (see matplotlib documentation for valid color names). Defaults are ‘r’ (red) and ‘k’ (black)

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)

stations_alpha, paths_alphafloat, optional

Transparency of the stations and of the great-circle paths. Defaults are None and 0.3

paths_widthfloat

Linewidth of the great-circle paths

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.

map_boundarieslist or tuple of floats, shape (4,), optional

Lonmin, lonmax, latmin, latmax (in degrees) defining the extent of the map

bound_mapbool

If True, the map boundaries will be automatically determined. Ignored if map_boundaries is not None

Returns

None if show is True. Otherwise a GeoAxesSubplot instance

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

Creates a maps of seismic receivers

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. Default is ‘110m’

kwargs :

Additional keyword arguments passed to matplotlib.pyplot.scatter

Returns

If show is True

None

Otherwise

ax, i.e. the GeoAxesSubplot

raypaths_per_pixel()

Calculates the number of raypaths in each pixel of the mesh stored in the SeismicTomography instance (self.grid.mesh)

Returns

ndarray of dimension self.grid.mesh.shape[0]

reduce_measurements(latmin=- 90, latmax=90, lonmin=- 180, lonmax=180)

Any measurement not intersecting the specified region is removed from the data. Consider using this function in presence of large datasets, for filtering out those rays that do not constrain the parameters (i.e., the grid cells) of interest in the inversion.

Parameters

lat1, lon1, lat2, lon2float (in degrees)

Boundaries of the region of interest. (-180<lon<180, -90<lat<90)

Returns

None

The instance of SeismicTomography is updated with the new data

refine_parameterization(hitcounts=100, keep_empty_cells=True, latmin=None, latmax=None, lonmin=None, lonmax=None)

Halves the dimension of the pixels intersected by a number of raypaths >= hitcounts.

Parameters

hitcountsint

Each parameter (grid cell) intersected by a number of raypaths equal or greater than this threshold is “refined”, i.e., splitted in four equal-area sub-parameters. Default is 100

keep_empty_cellsbool

If False, cells intersected by no raypaths are removed

lat1, lon1, lat2, lon2float

If specified (in degrees), only the region inside this boundary is considered for refining the parameterization. (-180<lon<180, -90<lat<90)

Returns

None

self.grid.mesh and self.A are updated

resolution_test(mesh, velocity_model, noise=0, ndamp=0, rdamp=0)

Resolution test on an arbitrary input model

This function is called under the hood to perform both the checkerboard and the spike test.

Parameters

meshndarray, shape (n, 4)

Grid cells in seislib format, consisting of n pixels described by lat1, lat2, lon1, lon2 (in degrees). (-180<lon<180, -90<lat<90)

velocity_modelndarray, shape (n,)

Synthetic velocity model associated with mesh

noisefloat, optional

If > 0, random noise expressed as percentage of the synthetic data

ndampfloat

Extent of norm damping (default is zero)

rdampfloat

Extent of roughness damping (default is zero)

Returns

dict

The dictionary is structured as follows:

• ‘synth_data’: ndarray of shape (m,), where m corresponds to the number of rows in the A matrix (i.e., in the data kernel)

• ‘synth_model’: inverse of velocity_model, used to create the synthetic data

• ‘retrieved_model’: model retrieved from the inversion of the synthetic data

• ‘mesh’: same as the input mesh

static roughness_operator(mesh, return_sparse=True)

Computes the roughness operator associated with the parameterization

Parameters

meshndarray, shape (n, 4)

Grid cells in seislib format, consisting of n pixels described by lat1, lat2, lon1, lon2 (in degrees). If None (default), the mesh stored in the SeismicTomography instance is used

return_sparsebool

If True, the matrix of roughness is converted to scipy.sparse.csr_matrix

Returns

ndarray of shape (n,n)

Roughness operator

save(path)

Saves the SeismicTomography instance to the specified path in the pickle compressed format.

Parameters

pathstr

Absolute path of the resulting file

solve(A=None, slowness=None, refvel=None, mesh=None, ndamp=0, rdamp=0)

Method for solving the (regularized) inverse problem Ax = b, where A is the matrix of coefficients, x the sought model of slowness, and b the array of measurements expressed in slowness. The method first computes the array of residuals r = b - Ax0, where x0 denotes the reference slowness passed to the function as refvel (and converted to slowness prior to the inversion). The residuals are then used to retrieve the least square solution

\[{\bf x}_{LS} = ({\bf A^{T}A} + \mu^2 {\bf I} + \rho^2 {\bf I G})^-1 \bf{A^{T} r},\]

where T denotes transposition, I the identity matrix, and \(\mu\) and \(\rho\) are norm and roughness damping, respectively. The square matrix G has the same dimensions of the model x. The above solution is then used to compute the final model \({\bf x = x0 + x}_{LS}\).

Parameters

Andarray, shape (m, n), optional

Matrix of the coefficients. If None (default), the A matrix stored in the SeismicTomography instance is used

slownessndarray, shape (m,), optional

Inverse of the velocity measurements. Corresponds to b in the equation Ax = b. If None (default), the velocity stored in the SeismicTomography instance is used to compute slowness

refvelfloat, optional

Reference velocity (in m/s) used for the initial guess of the velocity model. If provided, should stabilize the inversion and make it converge faster. If None (default), the reference velocity stored in the SeismicTomography instance is used

meshndarray, shape (n, 4), optional

Grid cells in seislib format, consisting of n pixels described by lat1, lat2, lon1, lon2 (in degrees). If None (default), the mesh stored in the SeismicTomography instance is used

ndampfloat

Norm damping coefficient (default is zero)

rdampfloat

Roughness damping coefficient (default is zero)

Returns

ndarray of shape (n,)

Least-square solution (slowness, in s/m)

static spike(ref_value, x0, y0, sigma_x=1, sigma_y=1, anom_amp=0.1)

Spike-like pattern on Earth, created via a 3-D Gaussian

Parameters

ref_valuefloat

Reference value with respect to which the anomalies are generated

x0, y0float

Central longitude and latitude of the anomaly

sigma_x, sigma_yfloat

Standard deviations in the longitude and latitude directions of the 3-D Gaussian. Default is 1

anom_ampfloat

Amplitude of the anomalies. Default is 0.1 (i.e. 10% of the reference value)

Returns

Function

Generating function that should be called passing the values of longitude and latitude where the values of the spike pattern want to be retrieved

Examples

>>> from seislib.tomography import SeismicTomography
>>> tomo = SeismicTomography(1, latmin=-20, latmax=20, lonmin=-20, lonmax=20)
>>> lat = np.mean(tomo.grid.mesh[:, :2], axis=1)
>>> lon = np.mean(tomo.grid.mesh[:, 2:], axis=1)
>>> pattern = tomo.spike(ref_value=3, 
...                      x0=0, 
...                      y0=0, 
...                      sigma_x=5, 
...                      sigma_y=5,
...                      anom_amp=0.2)
>>> c = pattern(lon, lat)
>>> img, cb = tomo.plot_map(tomo.grid.mesh, c, projection='Robinson', 
...                         show=False)
>>> cb.set_label('Spike Anomaly')

spike_test(x0, y0, sigma_x, sigma_y, regular_grid=False, latmin=None, latmax=None, lonmin=None, lonmax=None, cell_size=5, anom_amp=0.1, refvel=None, noise=0, ndamp=0, rdamp=0)

Resolution test, known as “spike test”. The method first builds synthetic data (velocities) using a spike pattern as an input (velocity) model. The spike is built using SeismicTomography.spike.

Parameters

x0, y0float

Central longitude and latitude of the anomaly. (-180<lon<180, -90<lat<90)

sigma_x, sigma_yfloat (default is 1)

Standard deviations in the longitude and latitude directions of the 3-D Gaussian

regular_gridbool

If False (default), the study area is discretized using an equal-area parameterization. Otherwise, a regular grid is employed.

latmin, latmax, lonmin, lonmaxfloat, optional

Boundaries of the checkerboard, in degrees

cell_sizeint | float

Size of the checkerboard grid-cell sides (in degrees). Larger values are preferrable to avoid the resolution of an overly large inverse problem

anom_ampfloat

Intensity of the anomalies with respect to ref_value. A value of 0.1 means 10% of ref_value

refvelfloat, optional

Original values of the anomalies (black and white squares of a checkerboard). By default, their value is +-1. Therefore, if anom_extent is 0.1 (say) their final values will be +-0.1

noisefloat, optional

If > 0, random noise expressed as percentage of the synthetic data

ndampfloat

Extent of norm damping (default is zero)

rdampfloat

Extent of roughness damping (default is zero)

Returns

Dictionary object structured as follows:

• ‘synth_data’: ndarray of shape (m,), where m corresponds to the number of rows in the A matrix (i.e., in the data kernel)

• ‘synth_model’: ndarray of shape (n,), synthetic model (slowness) used to create the synthetic data

• ‘retrieved_model’: model retrieved from the inversion of the synthetic data

• ‘mesh’: ndarray of shape (n, 4), parameterization associated with the synthetic model

update_data_info(refvel=None)

Updates the information on minimum / maximum longitude and latitude of the data and the reference velocity

Parameters

refvelfloat, optional

If not None, the reference velocity is updated along with minimum and maximum longitudes and latitudes

classmethod velocity_to_delay(lat1, lon1, lat2, lon2, velocity, refvel)

Converts a velocity to delay

Velocities (v) faster than the reference (v0) correspond to negative delays (dt):

dt = x/v - x/v0

where x denotes distance.

Parameters

lat1, lon1, lat2, lon2float or array-like of shape (n,)

Coordinates (in radians) associated with the velocity

velocityfloat or array-like of shape (n,)

Measured velocity (in m/s)

refvelfloat

Reference velocity to be used to calculate the delay (in m/s)

Returns

float or array-like of shape (n,)