Visualization

After 10 iterations of inversion, we can visualize model model_M09.h5 in the optimize directory. The final model is stored as a 3-D regular grid in HDF5 format, identical to the initial model. Here we use PyGMT to visualize the recovered velocity perturbation.

Read HDF5 model

Read the final and initial models and compute Vs perturbation relative to the starting model.

Python


import h5py
import numpy as np
import xarray as xr
import pygmt
 
with h5py.File('./optimize/model_M09.h5', 'r') as f:
    x = f['x'][:]
    y = f['y'][:]
    z = f['z'][:]
    vsf = f['vs'][:]
    vpf = f['vp'][:]
 
with h5py.File('./optimize/model_M00.h5', 'r') as f:
    vsi = f['vs'][:]
    vpi = f['vp'][:]
 
dvs = 100 * (vsf - vsi) / vsi
dvp = 100 * (vpf - vpi) / vpi

Plot depth slices of velocity perturbation

Select depth slices at 20 km and 50 km to visualize the recovered Vs anomaly.

Python


depth = [-20000, -50000]
 
fig = pygmt.Figure()
with fig.subplot(nrows=1, ncols=2, figsize=("14c", "7c"),
                 margins=["0.5c", "0.5c"],
                 frame=['xaf+lX (m)', 'yaf+lY (m)', 'WSne'],
                 sharey=True):
    for i, dep in enumerate(depth):
        iz = np.argmin(np.abs(z - dep))
        fig.basemap(
            region=[x.min(), x.max(), y.min(), y.max()],
            projection="x?c",
            panel=True
        )
        grid = xr.DataArray(
            dvs[:, :, iz].T,
            coords=[y, x],
            dims=["y", "x"],
        )
        pygmt.makecpt(
            cmap="polar",
            series=[-10, 10, 0.1],
            continuous=True,
            reverse=True,
            background=True
        )
        fig.grdimage(
            grid=grid,
            cmap=True,
        )
        fig.text(
            text="%.1f km" % (dep / -1000),
            position="TR",
            font="12p,Helvetica-Bold",
            offset="-0.2c/-0.2c",
        )
fig.colorbar(frame=["a5f1", "x+ldVs (%)"])
fig.show()

Plot cross-sections of velocity perturbation

Select a cross-section at the centre of the model (Y = middle of the array) to examine the depth structure of the recovered anomaly.

Python


ysec = y[len(y) // 2]
 
fig = pygmt.Figure()
fig.basemap(
    region=[x.min(), x.max(), z.min(), z.max()],
    projection="x0.00005c/0.0001c",
    frame=["WSne", "xaf+lX (m)", "yaf+lZ (m)"],
)
iy = np.argmin(np.abs(y - ysec))
grid = xr.DataArray(
    dvs[:, iy, :].T,
    coords=[z, x],
    dims=["y", "x"],
)
pygmt.makecpt(
    cmap="polar",
    series=[-10, 10, 0.1],
    continuous=True,
    reverse=True,
    background=True
)
fig.grdimage(
    grid=grid,
    cmap=True,
)
fig.colorbar(frame=["a5f1", "x+ldVs (%)"])
fig.text(
    text=f"Y = {ysec:.0f} m",
    position="BR",
    font="12p,Helvetica-Bold",
    offset="-0.2c/0.2c",
)
fig.show()