Mondaic
This tutorial is presented as Python code running inside a Jupyter Notebook, the recommended way to use Salvus. To run it yourself you can copy/type each individual cell or directly download the full notebook, including all required files.
%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_gridIn 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:
- The elements enclosing each point are found, and
- 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 0x7fa8e790e410>
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[2024-03-15 09:08:22,206] INFO: Creating mesh. Hang on.
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fa8d99e3450>
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 0x7fa8d1b31b10>
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 0x7fa8d1bbf210>]
# 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[2024-03-15 09:08:24,051] INFO: Creating mesh. Hang on.
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fa8d1a35090>
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 0x7fa8d1aada90>
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 0x7fa8d0244810>
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[2024-03-15 09:08:26,610] INFO: Creating mesh. Hang on.
Interpolating model: mantle.
Interpolating model: crust.
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7fa8d19cfe90>
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()[2024-03-15 09:08:29,740] WARNING - salvus.mesh.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 0x7fa8ce4ed110>
PAGE CONTENTS