API

Modules

pytomofilt

pytomofilt.cli

pytomofilt.filter_models

Functions to filter geodynamic models

These are used to provide a command line interface to common filtering tasks.

class pytomofilt.filter_models.TomographicModelFiles(name: str, base_path: Path, coef_file: Path, weights_file: Path, evec_file: Path, damping: float | None = 0.002)

Input files for a tomographic model

pytomofilt.filter_models.ptf_reparam_filter_files(tomographic_model: Path, geodynamic_model: Path)

Reparameterize and filter a geodynamic model so it can be compared with tomography

pytomofilt.filter_models.tomographic_model_from_name(model_name: str) → TomographicModelFiles

tomographic_model_from_name builds a tomographic model with the name as input

Parameters:

model_name – Name of tomographic model

Returns:

Paths to each file needed for tomographic filtering

Return type:

TomographicModelFiles

pytomofilt.filter_models.tomographic_model_from_path(input_path: Path) → TomographicModelFiles

tomographic_model_from_path builds and validates a tomographic model exists in input_path

Parameters:

input_path (Path) – Location of tomographic files (directory)

Returns:

Paths to each file needed for tomographic filtering

Return type:

TomographicModelFiles

pytomofilt.filter

class pytomofilt.filter.Filter(evec_file, wght_file, damping, verbose=False)

Object to apply tomographic resolution operator to seismic model

The creation of tomographic models of seismic wave velocities such as S20RTS, S40RTS and SP12RTS can also provide a resolution operator that maps a true earth model into the tomography model, and fully describes the spatially heterogeneous resolution of the tomography. This object allows a resolution operator from the “RTS” family of models to be applied to a model of seismic velocities, such that the velocities can be fairly compared with the tomography. Software for doing this also exists in Fortran. The main difference between this implementation and the Fortran version is that this version holds the operator in memory allowing multiple models to be “filtered”. The Fortran implementation reads applies one “line” of the operator at a time reading the line from disk. This means that this implementation uses much more memory but is much faster if multiple models have to be filtered. On the other hand, the Fortran implementation can be used on low memory machines. Within LEMA where we make use of very large numbers of low resolution models holding the model in memory is clearly the correct choice (indeed, it is ~100 times faster to do it this way).

To use the filter first create an instance of this class while providing the file names of the files containing the resolution operator (this is quite slow, see the documentation of the __init__ method). Then convert your model of seismic velocities into an RTSJointShell object before calling the RTSJointShell.filter_sp_model method with the instance provided in the “filter” keyword. The RTSJointShell object is then updated to represent the “filtered” seismic velocities.

apply_filter(x)

Applies the tomography filter

Reimplementation of the guts of mk3d_res_ESEP.f but with array operations when this is easy.

Arguments: x. The input model as a 1D vector. If the model is represented by a set of values varying in radius, r, degree, l and order, m, with real and imaginary components the order of the elements (r, l, m, s) in the vector is given by: (0, 0, 0, r), (0, 1, 0, r), (0, 1, 1, r), (0, 1, 1, i), (0, 2, 0, r), (0, 2, 1, r), (0, 2, 1, i), (0, 2, 2, r), (0, 2, 2, i), (0, 3, 0, r), (0, 3, 1, r), (0, 3, 1, i) … (0, lmax, lamx, r), (0, lmax, lmax, i) … (rmax, lmax, lmax, r), (rmax, lmax, lmax, i). That is, imaginary coeffs are skipped when they do not exist and the data increments through l, m, then r, in that order.

Returns: xout. The filtered model (same structure as x)

read_eigenvec_file(filename)

Read evec file for SP12RTS filtering

This contains the eigenvectors and eigenvalues of the model resolution matrix. The file is an unformatted fortran file, so we read using scipy.io.FortranFile and note that this format does not form part of any standard so is (Fortran) processor dependent. We rely on the format being readable to us. There are three ways of terminating reading. (1) We have read all lenatd eigenvectors, (2) the eigenvalues (which are sorted) is smaller that a (damping dependent value) or (3) fewer than this number of vectors are in the file.

read_wgts_file(filename)

Read the weights file

Again unformatted. One interesting feature is that the type does not seem to be real*8, reads OK as integer… … so I assume that integer (default prec) is OK and then cast to a float in python (implicit type of twts should be real).

pytomofilt.model

class pytomofilt.model.RealLayer(radius: float, lats: list[float], lons: list[float], vals: list[float])

class pytomofilt.model.RealLayerModel(layers: list[pytomofilt.model.RealLayer])

plot_radial_slice(layer_no)

Calls the plotting.plot_grid function to plot velocity variations for one the layers.

Parameters

layer_noint

Layer number that is to be plotted.

pytomofilt.plotting

pytomofilt.plotting.plot_grid(lons: ~numpy.ndarray, lats: ~numpy.ndarray, vals: ~numpy.ndarray, r: float | None = None, title: str = '', quantity: str = '', cmap: str = 'RdBu', coast_color: str = 'k', projection: ~cartopy.crs.CRS = <Projected CRS: +proj=robin +a=6378137.0 +lon_0=0 +no_defs +type=c ...> Name: unknown Axis Info [cartesian]: - E[east]: Easting (metre) - N[north]: Northing (metre) Area of Use: - undefined Coordinate Operation: - name: unknown - method: Robinson Datum: unknown - Ellipsoid: unknown - Prime Meridian: Greenwich, scale_factor: float | None = 1.0, levels: ~numpy.ndarray | None = None, vmin: float | None = None, vmax: float | None = None, extend: str = 'neither', fig: ~matplotlib.figure.Figure | None = None, ax: ~matplotlib.axes._axes.Axes | None = None, figsize: ~typing.Tuple[int, int] = (15, 5), labelsize: int = 15) → Tuple[Figure, Axes, Any]

Plot a grid of values defined at (lons, lats). If fig and ax are supplied, then the plot is added to the provided matplotlib figure and axis handles, respectively.

Parameters

lons, latsarray-like

1D arrays of point coordinates (degrees).

valuesarray-like

Values corresponding to the coordinates.

rfloat, optional

Radius (km) to display in title.

titlestr

Title for the plot.

quantitystr

Label for colorbar.

cmapstr

Colormap for plotting grid filled contours. Default is RdBu

coast_colorstr

Color for coastlines. Default is black.

projectionccrs.CRS

Geographic project used for figure. Default is ccrs.Robinson()

scale_factor: float, optional

Scale factor multiplied to the grid before plotting. Default is 1

levels, vmin, vmax, extend :

Same as those in plt.contourf

figplt.Figure

Figure object of the plot.

axplt.Axes

Axes for the plot.

figsizetuple

Figure size. Default is (15,5)

labelsizeint

Labelsize for axes and colorbar ticklabels. Default is 15.

Returns

fig, ax, mappable

pytomofilt.plotting.plot_heatmap(image: ndarray, yticks: list, aspect_ratio: float = 0.005, title: str = '', figsize: Tuple[int, int] = (12, 12), fontsize: int = 12)

Plot a 2D raster image, used for plotting the correlation or spectra of models in spherical harmonics. The x-axis represents the spherical harmonic degree (l) and the y-axis represents radius (km). The image and yticks should be ordered in layers from shallowest to deepest.

Parameters

imagenp.ndarray

2D array representing the raster image to be displayed.

yticks: list (n_layers,)

A list of labels for the y-axis ticks corresponding to each layer.

aspect_ratiofloat, optional

Aspect ratio for the image display. Default is 0.005.

titlestr, optional

Title for the heatmap plot. Default is an empty string.

figsizeTuple[int, int], optional

Figure size as (width, height) in inches. Default is (12, 12).

fontsizeint, optional

Font size for axis labels and title. Default is 12.

pytomofilt.plotting.plot_shcoefs(coefs: ~numpy.ndarray, r: float | None = None, title: str = '', quantity: str = '', cmap: str = 'RdBu', coast_color: str = 'k', projection: ~cartopy.crs.CRS = <Projected CRS: +proj=robin +a=6378137.0 +lon_0=0 +no_defs +type=c ...> Name: unknown Axis Info [cartesian]: - E[east]: Easting (metre) - N[north]: Northing (metre) Area of Use: - undefined Coordinate Operation: - name: unknown - method: Robinson Datum: unknown - Ellipsoid: unknown - Prime Meridian: Greenwich, scale_factor: float | None = 1.0, lmax: int | None = None, levels: ~numpy.ndarray | None = None, vmin: float | None = None, vmax: float | None = None, extend: str = 'neither', fig: ~matplotlib.figure.Figure | None = None, ax: ~matplotlib.axes._axes.Axes | None = None, figsize: ~typing.Tuple[int, int] = (15, 5), labelsize: int = 15) → Tuple[Figure, Axes, Any]

Plot a set of spherical harmonic coefficients coefs. Label the radius with r (km) and the value with quantity. Harmonics are truncated at degree lmax. If fig and ax are supplied, then the plot is added to the provided matplotlib figure and axis handles, respectively.

Parameters

coefsnp.ndarray

RTS-format coefficients (as in your code). Shape expected to match sh_tools.rts_to_sh.

rfloat, optional

Radius (km) to display in title.

titlestr

Title for the plot.

quantitystr

Label for colorbar.

cmapstr

Colormap for plotting grid filled contours. Default is RdBu

coast_colorstr

Color for coastlines. Default is black.

projectionccrs.CRS

Geographic project used for figure. Default is ccrs.Robinson()

scale_factor: float, optional

Scale factor multiplied to the grid before plotting. Default is 1

levels, vmin, vmax, extend :

Same as those in plt.contourf

figplt.Figure

Figure object of the plot.

axplt.Axes

Axes for the plot.

figsizetuple

Figure size. Default is (15,5)

labelsizeint

Labelsize for axes and colorbar ticklabels. Default is 15.

Returns

fig, ax, mappable