Mondaic Documentation
www.mondaic.com
Mondaic
Version:
No results found.
Copy
%matplotlib inline

import numpy as np
import xarray as xr
import salvus.namespace as sn
from salvus.mesh.unstructured_mesh_utils import extract_model_to_regular_grid

Working with variable models and extracting slices

Micro tutorial

In this tutorial we'll make a few simple meshes, and see how we can extract the model parameters on those meshes for visualization and analysis. This will allow us to visualize the meshed models along different depth slices, as well as convert them to a regularly-sampled representation that may be easier to share with others.
Salvus uses xarray to handle models defined on regular grid. This allows for a simple way to interpolate models onto unstructured grids. We recommend reading xarray's documentation for a deeper understanding how its indexing, plotting, and interpolation functions work.
On the Salvus side, the extraction algorithm consists of two critical steps:
  1. The elements enclosing each point are found, and
  2. The material parameters are interpolated to each point using the mesh's GLL basis.
Since the mesh may be fully unstructured, the enclosing element cannot be found via a simple coordinate lookup, and as such the runtime of the extraction routine scales with the size of the inputs. For many cases (>10k elements, >100k points), however, a runtime of several seconds can be expected. A progress bar is available by passing verbose=True to the extract_model_to_regular_grid function.
Points passed which are outside of the mesh are labeled with np.nan, following the standard conventions for missing data. Note that having many points outside of the mesh may negatively affect the enclosing-element algorithm, so it is good to keep these more-or-less in line with the dimensions of your domain. This is not always the case when topography is present. In these cases you may see a message reporting "N points were not claimed by enclosing elements". You may want to adjust the extraction points to keep these at a a minimum, but with deformed meshes or those not aligned with coordinate axes some points outside the mesh are inevitable, and the message need not be a concern.
Each section here -- 2-D, 3-D, and 3-D Spherical -- can be understood independently of the others.
# Gaussian
x, y = np.linspace(0.0, 1.0, 101), np.linspace(0.0, 1.0, 101)
xx, yy = np.meshgrid(x, y, indexing="xy")
g = np.exp(-(((xx - 0.5) ** 2 + (yy - 0.5) ** 2) / (2 * 0.2**2)))

# Pars
vp = 2 * g + 1
rho = vp / 2

# Xarray dataset
ds_model_2d = xr.Dataset(
    coords={"x": x, "y": y},
    data_vars={"VP": (["x", "y"], vp), "RHO": (["x", "y"], rho)},
)

# Salvus wrapper.
m = sn.model.volume.cartesian.GenericModel(name="blob", data=ds_model_2d)

# Plot
m.ds.VP.plot()
<matplotlib.collections.QuadMesh at 0x7e1f015db5d0>
p_2d = sn.Project.from_volume_model("proj_2d", m, True)
The simulation configuration here will also be used for the 3-D example.
sc = sn.SimulationConfiguration(
    name="blob_fmax",
    event_configuration=None,
    elements_per_wavelength=2.0,
    max_frequency_in_hertz=10.0,
    model_configuration=sn.ModelConfiguration(
        background_model=None, volume_models="blob"
    ),
)

p_2d.add_to_project(sc, overwrite=True)
mesh_0 = p_2d.simulations.get_mesh("blob_fmax")
mesh_0
[2025-04-15 19:37:04,782] INFO: Creating mesh. Hang on.

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7e1f00318d50>
You can also use xarray to extract a model from a mesh onto a regular grid you define. In this example, we pull out onto a regular grid the same model we put in.
ds_extract_2d = xr.Dataset(
    coords={"x": np.linspace(0, 1.0, 101), "y": np.linspace(0, 1.0, 101)}
)
ds_extract_2d = extract_model_to_regular_grid(
    mesh_0, ds_extract_2d, ["VP", "RHO"]
)
ds_extract_2d.VP.plot()
<matplotlib.collections.QuadMesh at 0x7e1f1b6ec5d0>
The xarray interface can also be used to extract parameters along a general profile through the model.
ds_extract_2d_line = xr.Dataset(
    coords={"x": np.linspace(0.0, 1.0, 101), "y": [0.5]}
)
ds_extract_2d_line = extract_model_to_regular_grid(
    mesh_0, ds_extract_2d_line, ["VP", "RHO"]
)

ds_extract_2d_line.VP.plot()
[<matplotlib.lines.Line2D at 0x7e1f1b05ced0>]
In the 3-D the process is much the same as 2-D.
# Gaussian
x, y, z = (
    np.linspace(0.0, 1.0, 101),
    np.linspace(0.0, 1.0, 101),
    np.linspace(0.0, 1.0, 101),
)

xx, yy, zz = np.meshgrid(x, y, z, indexing="xy")
g = np.exp(
    -(((xx - 0.5) ** 2 + (yy - 0.5) ** 2 + (zz - 0.5) ** 2) / (2 * 0.2**2))
)

# Pars
vp_3d = 2 * g + 1
rho_3d = vp_3d / 2

# Xarray dataset
ds_model_3d = xr.Dataset(
    coords={"x": x, "y": y, "z": z},
    data_vars={
        "VP": (["x", "y", "z"], vp_3d),
        "RHO": (["x", "y", "z"], rho_3d),
    },
)

# Salvus wrapper.
m_3d = sn.model.volume.cartesian.GenericModel(name="blob", data=ds_model_3d)
p_3d = sn.Project.from_volume_model("proj_3d", m_3d, True)
Note we are using the same simulation configuration as in 2-D.
p_3d.add_to_project(sc, overwrite=True)
mesh_3d = p_3d.simulations.get_mesh("blob_fmax")
mesh_3d
[2025-04-15 19:37:05,700] INFO: Creating mesh. Hang on.

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7e1f1b5069d0>
Similar to the 2-D example, here we define a 3-D regular grid. We can interpret the below as defining 11 slices along the x-y plane at varying z-values.
ds_extract_3d = xr.Dataset(
    coords={
        "x": np.linspace(0, 1.0, 101),
        "y": np.linspace(0, 1.0, 101),
        "z": np.linspace(0, 1.0, 11),
    }
)
ds_extract_3d = extract_model_to_regular_grid(
    mesh_3d, ds_extract_3d, ["VP", "RHO"]
)

# Get an xy slice at z == 0.5.
ds_extract_3d.VP.sel(z=0.5).plot()
<matplotlib.collections.QuadMesh at 0x7e1f1b5b4510>
The process is much the same again in spherical coordinates. First we define our model, and wrap it in the Salvus model handling objects. Note here that we define separate mantle and crustal models. These can be used together on the same mesh.
# Spherical coordinates
lat = np.linspace(-5.0, 5.0, 101)
lon = np.linspace(-5.0, 5.0, 101)
dep = np.linspace(0.0, 660e3, 101)

# Gaussian
xx, yy, zz = np.meshgrid(lat, lon, dep, indexing="xy")
g = np.zeros_like(xx)
for lats in lat[25:-25:25]:
    for lons in lon[::25]:
        g += np.exp(-(((xx - lats) ** 2 + (yy - lons) ** 2) / (2 * 0.5**2)))

# Xarray dataset
ds = xr.Dataset(
    coords={"latitude": lat, "longitude": lon, "depth": dep},
    data_vars={
        "vp": (["latitude", "longitude", "depth"], g * 20, {"units": "%"}),
    },
    attrs={
        "geospatial_lon_units": "degrees",
        "geospatial_lat_units": "degrees_north",
        "geospatial_vertical_units": "m",
    },
)

ds_mantle = ds.copy()
ds_mantle["vp"] *= -1

# Salvus wrapper.
mc = sn.model.volume.seismology.CrustalModel(name="crust", data=ds)
mm = sn.model.volume.seismology.MantleModel(name="mantle", data=ds)

# Plot depth slice.
mm.ds.vp.sel(depth=500e3, method="nearest").plot()
<matplotlib.collections.QuadMesh at 0x7e1f00143b50>
d_sph = sn.domain.dim3.SphericalChunkDomain(
    lat_center=0.0,
    lon_center=0.0,
    lat_extent=5.0,
    lon_extent=5.0,
    radius_in_meter=6371e3,
)

# Add mantle and crustal model.
p_sph = sn.Project.from_domain("proj_sph", d_sph, True)
for m in [mc, mm]:
    p_sph.add_to_project(m, overwrite=True)
Note here that we are specifying a new simulation configuration. This is because we are now meshing for 50 seconds, rather than 10 hertz as above.
p_sph.add_to_project(
    sn.SimulationConfiguration(
        name="blobs",
        tensor_order=2,
        event_configuration=None,
        max_depth_in_meters=660e3,
        min_period_in_seconds=50.0,
        elements_per_wavelength=2.0,
        model_configuration=sn.ModelConfiguration(
            background_model=sn.model.background.one_dimensional.BuiltIn(
                name="prem_iso_one_crust"
            ),
            volume_models=["mantle", "crust"],
        ),
    )
)
mesh_sph = p_sph.simulations.get_mesh("blobs")
mesh_sph
[2025-04-15 19:37:07,992] INFO: Creating mesh. Hang on.
Interpolating model: mantle.

Interpolating model: crust.

<salvus.mesh.data_structures.unstructured_mesh.unstructured_mesh.UnstructuredMesh object at 0x7e1f1b029050>
In this case, we define our xarray using spherical coordinates. Here we are extracting 11 depth slices, regularly spaced from 0 to 660 km in depth.
ds_extract_3d_sph = xr.Dataset(
    coords={
        "longitude": np.linspace(-5.0, +5.0, 101),
        "latitude": np.linspace(-5.0, 5.0, 101),
        "depth": np.linspace(0, 660e3, 11),
    },
    attrs={"radius_in_meters": 6371e3},
)
ds_extract_3d_sph = extract_model_to_regular_grid(
    mesh_sph, ds_extract_3d_sph, ["VP", "RHO"], verbose=True
)

# Get a lat/lon slice at depth == 500e3.
ds_extract_3d_sph.VP.sel(depth=500e3, method="nearest").plot()
[2025-04-15 19:37:11,168] WARNING - salvus.mesh.algorithms.unstructured_mesh.utils: 1936 points were not claimed by enclosing elements. Depending on your use case, this may not be an issue.

<matplotlib.collections.QuadMesh at 0x7e1eff7c8190>
Previous
Layered models
Next
Marmousi
We use cookies to analyze our traffic, but won't track you until explicit consent is given. Using our contact form will create relevant cookies and transmit data to third party services to keep track of your request. Read more about how we use personal information in our Privacy Policy.
PAGE CONTENTS