Skip to content

Model Interpolation

Jupyter Notebook

This page is automatically generated from a Jupyter Notebook.

Note that everything here is fully automatically tested against - and thus guaranteed to work - only for the latest versions of all the Salvus packages. So please update if anything works differently on your machine.

Click here for a version that is compatbile with the 0.9.X releases.

Model Interpolation

This notebook shows how to interpolate a given model defined on a regular grid to a Salvus mesh using the Marmousi model as an example.

# initialize notebook
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import obspy
from salvus_mesh.structured_grid_2D import StructuredGrid2D

plt.style.use('ggplot')
plt.rcParams['figure.figsize'] = (10, 8)

The Marmousi model is provided in a segy file, which we open with obspy and map to numpy arrays.

def get_marmousi():
    nx, ny = 2301, 751
    dx, dy = 4., 4.

    rho = np.empty((ny, nx))
    st = obspy.read('data/marmousi_density.segy.gz')
    for _i, tr in enumerate(st):
        rho[_i, :] = tr.data

    vp = np.empty((ny, nx))
    st = obspy.read('data/marmousi_velocity.segy.gz')
    for _i, tr in enumerate(st):
        vp[_i, :] = tr.data

    x = np.arange(nx) * dx
    y = np.arange(ny) * dy

    return x, y, vp, rho

Plot the model to verify we got it in correctly:

x, y, vp, rho = get_marmousi()
fig, axes = plt.subplots(2, 1)
axes[0].imshow(vp)
axes[1].imshow(rho)
<matplotlib.image.AxesImage at 0x7f6e669046d8>

png

Build a rectilinear mesh for a subregion of the model. While this might be simplistic, note that the model interpolation afterwards would be the same for any unstructured mesh.

fmax = 10.                    # maximum frequency in Hz
elements_per_wavelength = 2.  # resolution criterion

# region in the marmousi model to mesh
mesh_min_x = x.max() * 0.2
mesh_max_x = x.max() * 0.7
mesh_min_y = y.min()
mesh_max_y = y.max()

hmax = vp.min() / fmax / elements_per_wavelength
nelem_x = int((mesh_max_x - mesh_min_x) / hmax) + 1
nelem_y = int((mesh_max_y - mesh_min_y) / hmax) + 1

sg = StructuredGrid2D.rectangle(nelem_x=nelem_x, nelem_y=nelem_y, 
                                min_x=mesh_min_x, max_x=mesh_max_x, 
                                min_y=mesh_min_y, max_y=mesh_max_y)
m = sg.get_unstructured_mesh()
m.change_tensor_order(4)

Interpolate the model to the mesh's element nodes using a scipy function and attach it to the mesh for plotting:

from scipy.interpolate import RectBivariateSpline

f_vp = RectBivariateSpline(x, y, vp.T[:, ::-1])
f_rho = RectBivariateSpline(x, y, rho.T[:, ::-1])

nodes = m.get_element_nodes()

nodes_x = nodes[:, :, 0]
nodes_y = nodes[:, :, 1]

vp_nodes = f_vp(nodes_x, nodes_y, grid=False)
rho_nodes = f_rho(nodes_x, nodes_y, grid=False)

m.attach_field('VP', vp_nodes)
m.attach_field('RHO', rho_nodes)
m.attach_field('fluid', np.ones(m.nelem))
m

<salvus_mesh.unstructured_mesh.UnstructuredMesh at 0x7f6e670e3c88>