salvus.toolbox.toolbox

def create_events_from_segy_file(
    segy_file: pathlib.Path | str,
    component: str,
    fields: typing.Iterable[str] | str,
    source_callback_function: typing.Callable,
    receiver_callback_function: typing.Callable,
    dim: int,
    event_suffix: str = "",
    discrimination_points: set[int] | int | None = None,
    scale_coordinates_by: float = 1.0,
    event_discriminator: str = "EnergySourcePoint",
    source_burial_depth_in_meters: float = 0.0,
    receiver_burial_depth_in_meters: float = 0.0,
) -> EventDataCollection: ...

Create a set of events from a SEG-Y file (experimental).

Currently only works for single-component receivers. Currently only burial depths of 0m are supported.

This experimental action adds events from a SEG-Y file to a project. An “event” in this context is the collection of sources and receivers which share a common trace header (event_discriminator); multiple events can be extracted from a single SEG-Y file. Source locations are derived from the trace headers “SourceX” and “SourceY”. Receiver locations are derived from “GroupX” and “GroupY”. Events will be named by following the pattern “event_{EnergySourcePoint}_{event_suffix}”. Both sources and receivers will be of the “SideSet” variety, with side-set tags “z1” (or “y1” in 2-D). “SourceY” and “GroupY” are ignored if the source is 2-D.

The ‘segyio’ package is required to use this functionality. In practice, SEG-Y files can be extremely large, and you can expect this function to take a bit of time in those cases.

def detect_fluid(
    um: unstructured_mesh.UnstructuredMesh,
) -> unstructured_mesh.UnstructuredMesh: ...

Automatically detect fluid regions.

Automatically detects where fluid regions are in a coupled mesh, or marks the entire mesh as fluid if no elastic-only parameters are found.

def gen_taper_spherical(
    ds: xr.Dataset,
    deg_lat: float | None,
    deg_lon: float | None,
    use_nan_mask: str | None = None,
) -> xr.Dataset: ...

Generate an approximate lat/lon taper in spherical coordinates.

Often times models do not taper smoothly into the background model being used. This function can help you generate a taper function which will be used during the model interpolation to smoothly transition from the 3-D model to the values already in the mesh. Two modes are accessible here: either one which generates a taper from the spherical bounding box of the model, or one which uses ‘not-a-number’ (NaN) values in the dataset to restrict the taper to regions where the model is not defined. This can be useful when considering models from the IRIS ERC, as these models often contain NaNs as a means of masking out certain regions.

In the case of a CF model with NaNs (such as those supplied by the IRIS EMC), one could consider just passing 0.0 for both deg_lat and deg_lon. This simply masks out model regions where the model is not defined, and has been found to be quite robust even in models with complicated and unstructured distributions of NaNs. Nevertheless, the taper widths are always model dependent, and it helps to try a few values before settling on something final.

def generate_mesh_from_splines_2d(
    x_min: float,
    x_max: float,
    splines: list[tuple[interp1d, interp1d]],
    maximum_frequency: float,
    slowest_velocities: list[float],
    model_order: int = 2,
    absorbing_boundaries: tuple[list[str], float] | None = None,
    elements_per_wavelength: float = 1.5,
    use_refinements: bool = True,
    extrude_splines: bool = False,
) -> tuple[list[unstructured_mesh.UnstructuredMesh], float]: ...

Generate a 2-D mesh with interfaces defined by 1-D splines.

A common use case is to run simulations on a 2-D model with deformed internal discontinuities. This function helps to mesh such a model, and takes care to add space for absorbing boundaries, and adjust the vertical elemenet size within each layer. Given a list of n splines, a mesh with n

• 1 layers (chunks) will be returned. These chunks can then be summed together to create a continuous mesh.

The splines must be provided in a top down order. That is, the first set of points defines the top of the mesh (sideset y1), while the last set defines the bottom (sideset y0).

def geographic_to_geocentric_latitude(
    lat: float, flattening: float = 0.0033528106647474805
) -> float: ...

Return a geocentric latitude when given geographic latitude.

def get_interpolating_splines(
    layers_x: list[np.ndarray],
    layers_y: list[np.ndarray],
    kind: list[str] | str,
) -> list[tuple[interp1d, interp1d]]: ...

Get a list of tuples of interpolating splines from x and y values.

The list contains pairs of interpolating functions defining the top and bottom of a given layer. layers_x should be a 1-D numpy array which contains the x values of spline points, while layers_y should be 2-D numpy array which contains the y values definining each layer. If layers_y is an array with shape (n, m) then we have n layers which are each defined by m tie points.

For this function, horizons should not cross each other.

def get_linear_feature_mask(
    mesh: unstructured_mesh.UnstructuredMesh,
    normal: list[float] | npt.NDArray,
    point: list[float] | npt.NDArray,
    length: float,
    buffer: float,
) -> npt.NDArray: ...

Return an array which can be used to mask a linear feature.

It is often desirable to refine the mesh around one or more linear features in the model. For instance, one may be interested in a region of anomalously slow velocites, or an unconformity (i.e. a fault). Given a parameterized line, this function returns a mask which can later be used to refine, or otherwise modify, a mesh in the vicinity of the linear feature.

def get_shotgather(
    filename: str | pathlib.Path, field: str, cmp: int = 0, axis: int = 0
) -> tuple[np.ndarray, float, tuple[float, float, float, float]]: ...

Get a shotgather from a Salvus HDF5 receiver file.

Populates a 2-D numpy array with one component from a Salvus HDF5 receiver file.

def interpolate_cartesian(
    data: xr.Dataset,
    mesh: unstructured_mesh.UnstructuredMesh,
    extrapolate: bool = True,
    mode: str | None = None,
    mask: np.ndarray | tuple[str, int] | list[tuple[str, int]] | None = None,
) -> unstructured_mesh.UnstructuredMesh: ...

Interpolate a volumetric model onto a Cartesian mesh.

Either scipy.interpolate.RectBivariateSpline, or scipy.interpolate.RegularGridInterpolater, will be used, depending on the value passed for ‘mode’. "spline" is only accepted for 2-D meshes and models.

def interpolate_spherical(
    data: xr.Dataset,
    mesh: unstructured_mesh.UnstructuredMesh,
    mode: str = "perturbation_to_background",
    spherical_max_radius: float = 6371000.0,
    layers: list[int] | str | None = None,
    globe: bool = False,
    ignore: tuple[str, int] | None = None,
    interpolate_parameters: list[str] | None = None,
    verbosity: int = 0,
) -> unstructured_mesh.UnstructuredMesh: ...

Do a model interpolation in spherical coordinates.

Use scipy.interpolate.RectSphereBivariateSpline to create a spherical spline interpolation over the surface of each radial layer. Then use linear interpolation between these radial layers to get values at a general unstructured depth. Function is vectorized, so it should be pretty quick.

Note that the interpolation is done assuming the radius from a 1-D background model. This means factors such as surface topography are ignored when performing the model interpolation.

def iris_emc_to_xarray(
    filename: str | pathlib.Path | xr.Dataset,
    planet_radius_in_km: float = 6371.0,
    fill_na: bool = True,
) -> tuple[xr.Dataset, str]: ...

Read a model from the IRIS EMC.

Most 3-D Earth models from the IRIS EMC are released as a NetCDF file in the CF-1.0 convention. This function will read in such a model, and return it as an Xarray Dataset to be used for interpolation onto a Salvus mesh. The exact formatting of models within the EMC is not entirely consistent, so it is possible that this function will fail for some models. It has been tested and found to work for many of the EMC models, and should be good at least as a starting point for contributing your own read function. The objects returned from this function can be used directly in the interpolate_spherical function, perhaps with some tapering required first (see gen_taper_spherical for information).

def load_cartesian_model(
    model: pathlib.Path | str | xr.Dataset,
) -> xr.Dataset: ...

Load a cartesian model.

def mesh_from_xarray(
    data: xr.Dataset,
    slowest_velocity: str,
    maximum_frequency: float,
    elements_per_wavelength: float,
    model_order: int = 4,
    absorbing_boundaries: tuple[list[str], float] | None = None,
    interpolation_mode: str | None = None,
) -> unstructured_mesh.UnstructuredMesh: ...

Generate a Salvus mesh from an xarray.Dataset object.

This function is meant to be used as an entry point and way to encapsulate the variety of mesh generation algorithms contained in SalvusMesh, and is optimized for regularly-gridded data. The data can be either 2- or 3-D, and must be contained within a xarray.Dataset object. This object needs to contain coordinate arrays labelled ‘x’ and ‘y’, or ‘x’, ‘y’, and ‘z’ if working in a Cartesian coordinate system. If working in a spherical coordinate system as is common in regional to global scale seismology, one must provide the parameters ‘radius’, ‘latitude’, and ‘longitude’. Material paramters stored as DataArrays within the Dataset should be named according to the Salvus parameter naming conventions.

This function will check the sanity of inputs, and will raise if something incorrect is detected. The sophistication and varitey of meshing algorithms contained within this function will evolve over time, so ensure to check back for updates.

def read_elastic_marmousi(
    directory: pathlib.Path, parameters: list[str], decimate: int = 10
) -> xr.Dataset: ...

Read the ‘AGL Elastic Marmousi model’ into an Xarray Dataset.

The resulting dataset suitable for use by SalvusMesh. To use another SEG-Y file, RSF file, or even your own regularly-gridded format, you need to modify the function below to accept your data, while ensuring that the structure of the returned xarray.Dataset remains the same.

def skeleton_regular_2d(
    data: xr.Dataset,
    slowest_velocity: str,
    maximum_frequency: float,
    elements_per_wavelength: float,
    absorbing_boundaries: tuple[list[str], float] | None = None,
    model_order: int = 1,
) -> unstructured_mesh.UnstructuredMesh: ...

Return a simple rectilinear skeleton mesh from Salvus’ mesher.

Given a string specifying the parameter with the slowest velocity, this function will return a mesh which ensures that at least the specified number of elements per wavelength are specified given the specified maximum frequency. This returns the simplest mesh possible — no local optimizations are done here to improve solver efficiency.

The mesh will be returned with elemental nodal parameters initialized to the minimum values of each of the parameters in data. As such, the resultant mesh is not ready for immediate consumption, and is instead meant to serve as a skeleton on which to interpolate a variable model.

def visualize_wavefield_2d(
    filename: str | pathlib.Path, field: str, cmp: int = 0
) -> tuple[tri.Triangulation, np.ndarray]: ...

Get an object that can be used to plot simple 2-D wavefields.

This method works by using the tricontourf method in matplotlib, and so comes with some limitations. For instance, if the mesh is not convex, the triangulation will fill in any holes and thus lead to unexpected results. This function is meant to quickly plot wavefields on simple domains for use in jupyter-notebooks and the like — for more detailed volumetric plotting (including in 3-D) please use Paraview.

def xy_to_xarray(model: str | pathlib.Path | xr.Dataset) -> xr.Dataset: ...

def xyz_to_xarray(model: str | pathlib.Path | xr.Dataset) -> xr.Dataset: ...

Read a file in NetCDF format into an Xarray Dataset.

When storing 3-D models parameterized on a regular grid, a simple NetCDF file in the xyz format is all that is required. A typical file should have something like the following structure (output of ncdump -h):

    dimensions:
            x = 100 ;
            y = 100 ;
            z = 10 ;
    variables:
            float vp(z, y, x) ;
                vp:units = "mps" ;
            float vs(z, y, x) ;
                vs:units = "mps" ;
            float rho(z, y, x) ;
                rho:units = "kgpcm" ;
            float x(x) ;
                x:units = "m" ;
                x:standard_name = "projection_x_coordinate" ;
                x:long_name = "x coordinate of projection" ;
            float y(y) ;
                y:units = "m" ;
                y:standard_name = "projection_y_coordinate" ;
                y:long_name = "y coordinate of projection" ;
            float z(z) ;
                z:units = "m" ;
                z:standard_name = "projection_z_coordinate" ;
                z:long_name = "z coordinate of projection" ;

Such a file can also be opened directly in Paraview using the built-in NetCDF reader. You can see an example of how to build such a file in Salvus’ online documentation, or the test suite for this package.