Visualization

The recovered model in OUTPUT_FILES_ph_aniso/final_model.h5 contains the isotropic Vs field together with the azimuthal-anisotropy magnitude and fast-axis direction. This page shows how to plot each of them with PyGMT.

Model structure

An anisotropic SurfATT output contains:

Dataset	Shape	Meaning
`x`, `y`, `z`	1-D	Longitude, latitude, depth (km)
`vs`	`(nx, ny, nz)`	Isotropic shear-wave velocity (km/s)
`gc`	`(nx, ny, nz)`	Azimuthal-anisotropy coefficient $G_c$
`gs`	`(nx, ny, nz)`	Azimuthal-anisotropy coefficient $G_s$
`g0`	`(nx, ny, nz)`	Azimuthal-anisotropy magnitude (fraction of Vs)
`theta`	`(nx, ny, nz)`	Fast-axis direction (radians, measured from east, counter-clockwise positive)

g0 and theta are already the physical quantities you want to plot — they are related to the azimuthal-anisotropy coefficients $G_c$ and $G_s$ by

G_c = g_0 \cos(2\theta), \qquad G_s = g_0 \sin(2\theta)

or, inverted,

g_0 = \sqrt{G_c^2 + G_s^2}, \qquad \theta = \tfrac{1}{2}\arctan(G_s / G_c)

Load the model


import h5py
import numpy as np
import pygmt
from pygmt.clib import Session
import xarray as xr
import pandas as pd
import io
 
with h5py.File("OUTPUT_FILES_ph_aniso/initial_model.h5") as f:
    vs_init = f["vs"][:]
    x = f["x"][:]
    y = f["y"][:]
    z = f["z"][:]
 
with h5py.File("OUTPUT_FILES_ph_aniso/final_model.h5") as f:
    vs    = f["vs"][:]
    g0    = f["g0"][:]
    theta = f["theta"][:]
 
# depth-dependent isotropic velocity perturbation (%)
dv = np.zeros_like(vs)
for iz in range(z.size):
    dv[:, :, iz] = (vs[:, :, iz] - np.mean(vs[:, :, iz])) / np.mean(vs[:, :, iz]) * 100
 
depths = [6, 20, 40, 64]
region = [x[0], x[-1], y[0], y[-1]]
xx, yy = np.meshgrid(x, y, indexing="ij")

Two optional auxiliary files are referenced below:

usbound.dat: study-region polygon used with GMT clip to mask outside the domain
boundaries.gmt: geological province outlines drawn as an overlay

Isotropic Vs maps


cpt = "vel_norm.cpt"
 
fig = pygmt.Figure()
pygmt.config(MAP_FRAME_TYPE="plain")
 
with fig.subplot(
    ncols=2, nrows=2,
    figsize=("14c", "18.5c"),
    margins=("0.2c", "0.1c"),
    sharex="b", sharey="l",
):
    for i, dep in enumerate(depths):
        iz = int(np.argmin(np.abs(z - dep)))
 
        grid = pygmt.surface(
            x=xx.ravel(), y=yy.ravel(), z=vs[:, :, iz].ravel(),
            region=region, spacing="0.1",
        )
 
        fig.basemap(region=region, projection="M?", frame=["a5f1"], panel=True)
        fig.coast(area_thresh=1_000_000, resolution="h", land=True)
 
        vmin = np.min(vs[:, :, iz]) - 0.05
        vmax = np.max(vs[:, :, iz]) + 0.05
        pygmt.makecpt(cmap=cpt, series=[vmin, vmax, 0.02], continuous=True)
 
        # optional: clip outside the study region
        with Session() as lib:
            lib.call_module("clip", "bound.dat")
        fig.grdimage(grid=grid, cmap=True)
        with Session() as lib:
            lib.call_module("clip", "-C")
 
        fig.coast(
            area_thresh=1_000_000, resolution="h",
            shorelines="0.3p", borders="2/0.1p,150",
        )
        fig.plot("boundaries.gmt", pen="0.8p,255")
        fig.text(
            text=f"{dep} km", position="RB",
            font="12p", offset="-0.1c/0.1c", fill="255",
        )
 
        with pygmt.config(FONT_LABEL="20p", FONT_ANNOT_PRIMARY="16p"):
            fig.colorbar(
                position="JBC+w4c/0.25c+h",
                frame=["xafg", "y+l@:8p:Vs (km/s)"],
            )
 
        # annotate provinces on the first panel
        if dep == depths[0]:
            for text, lon, lat, angle in [
                ("Col",  -117.0, 46.0,  0),
                ("GB",   -115.5, 41.0,  0),
                ("YS",   -109.5, 41.5,  0),
                ("CP",   -110.0, 36.0,  0),
                ("SN",   -118.6, 36.7, -55),
                ("GV",   -121.8, 38.0, -55),
                ("SBR",  -112.8, 34.0, -36),
            ]:
                fig.text(text=text, x=lon, y=lat, angle=angle,
                         font="12p,1,255=0.5p,black")
 
fig.savefig("wus_vs_maps.png")
fig.show()

Interpreting the isotropic maps

Depth	Expected features
6 km	Slow sedimentary basins: Central Valley (GV), Salton Trough, Great Basin (GB) playas
20 km	Felsic upper-crustal roots of the Sierra Nevada (SN) and Colorado Plateau (CP)
40 km	Moho transition — fast cratonic crust east of ~106°W; slow, thinned crust in the Basin and Range (SBR/GB)
64 km	Upper mantle — slow Yellowstone (YS) / Snake River Plain anomaly, fast Juan de Fuca slab beneath Cascadia

Province labels above follow standard western-US nomenclature — Col: Columbia Plateau, GB: Great Basin, YS: Yellowstone, CP: Colorado Plateau, SN: Sierra Nevada, GV: Great Valley, SBR: Southern Basin and Range.

Azimuthal anisotropy

The fast-axis directions are plotted as short ellipse “bars” (GMT symbol e), with the azimuth taken from theta, the major axis proportional to the anisotropy magnitude g0, and the minor axis fixed to a small value so the ellipses render as thin bars.


def to_xarray(data, x, y, key="z"):
    """Convert a 2-D NumPy array to an xarray.DataArray indexed on (lat, lon)."""
    return xr.DataArray(data.T, coords=[y, x], dims=["lat", "lon"], name=key)
 
dn           = 6          # sub-sampling step for the bar field
ani_scaling  = 10         # stretch factor for visual bar length
min_ax       = 0.06       # minor-axis length (makes ellipse look like a bar)
barcolor     = "BROWN4"
barpen       = "0.1p,200"
 
fig = pygmt.Figure()
pygmt.config(MAP_FRAME_TYPE="plain")
 
with fig.subplot(
    ncols=2, nrows=2,
    figsize=("14c", "17c"),
    margins=("0.2c", "0.1c"),
    sharex="b", sharey="l",
):
    pygmt.makecpt(cmap="lapaz", series=[0, 0.08, 0.01],
                  continuous=True, reverse=True)
 
    for i, dep in enumerate(depths):
        iz = int(np.argmin(np.abs(z - dep)))
 
        fig.basemap(region=region, projection="M?",
                    frame=["a5f1"], panel=i)
 
        # background: anisotropy magnitude
        grid_g0 = pygmt.surface(
            x=xx.ravel(), y=yy.ravel(), z=g0[:, :, iz].ravel(),
            region=region, spacing="0.1",
        )
        fig.grdimage(grid=grid_g0, cmap=True)
 
        # foreground: fast-axis bars (one ellipse per sub-sampled cell)
        grid_g0_samp = to_xarray(g0[::dn, ::dn, iz],    x[::dn], y[::dn])
        grid_theta   = to_xarray(theta[::dn, ::dn, iz], x[::dn], y[::dn])
        xyz_theta = pygmt.grd2xyz(grid_theta)
        xyz_g0    = pygmt.grd2xyz(grid_g0_samp)
 
        ani = pd.concat(
            [xyz_theta["lon"], xyz_theta["lat"],
             xyz_theta["z"], xyz_g0["z"] * ani_scaling],
            axis=1,
        )
        ani["minax"] = np.ones_like(ani["lon"]) * min_ax
 
        fig.plot(ani, style="e", fill=barcolor, pen=barpen)
 
        fig.coast(
            area_thresh=1_000_000, resolution="h",
            shorelines="0.3p", borders="2/0.1p,150",
        )
        fig.plot("boundaries.gmt", pen="0.8p,255")
        fig.text(
            text=f"{dep} km", position="RB",
            font="12p", offset="-0.1c/0.1c", fill="255",
        )
 
        # scale-bar legend on the first two panels
        if dep == depths[0]:
            legend_spec = io.StringIO(f"""
S 0.2c e 45/{0.03 * ani_scaling:.1f}/{min_ax} {barcolor} {barpen} 0.6c 3%
S 0.2c e 45/{0.06 * ani_scaling:.1f}/{min_ax} {barcolor} {barpen} 0.6c 6%
""")
            fig.legend(spec=legend_spec,
                       position="jBL+w1.5c/1.1c+l1.17",
                       box="+p0.5p+gwhite")
        if dep == depths[1]:
            legend_spec = io.StringIO(f"""
S 0.2c e 45/{0.02 * ani_scaling:.1f}/{min_ax} {barcolor} {barpen} 0.6c 2%
S 0.2c e 45/{0.04 * ani_scaling:.1f}/{min_ax} {barcolor} {barpen} 0.6c 4%
""")
            fig.legend(spec=legend_spec,
                       position="jBL+w1.5c/1.1c+l1.17",
                       box="+p0.5p+gwhite")
 
    pygmt.config(FONT_LABEL="20p", FONT_ANNOT_PRIMARY="12p")
    fig.colorbar(position="JRM+w8c/0.4c",
                 frame=["xag+lAnisotropy strength"])
 
fig.savefig("wus_ani_maps.png")
fig.show()

style="e" draws ellipses of the form angle/major/minor. Because the minor axis (min_ax = 0.06) is much smaller than the major axis (g0 × ani_scaling), each ellipse collapses into a thin bar whose orientation matches the fast-axis direction. The bars have no arrowheads because fast-axis directions have 180° ambiguity — a bar at 30° and one at 210° are physically equivalent.

Interpreting the anisotropy maps

Crustal depths (6–20 km) — anisotropy often aligns with regional extensional fabric (E–W in the Basin and Range) or transform motion (NW–SE along the San Andreas system).
Uppermost mantle (40–64 km) — fast axes often track absolute plate motion or asthenospheric flow deflected by slabs; expect coherent E-to-NE fast directions south of the Mendocino Triple Junction and trench-parallel directions beneath Cascadia.

Misfit Reduction

objective_function.txt is a whitespace-separated table with an iter column and a misfit column:


import pandas as pd
import matplotlib.pyplot as plt
 
ms = pd.read_csv('./OUTPUT_FILES_ph_aniso/objective_function.txt', sep=r'\s+')
 
plt.style.use('seaborn-v0_8')
# print(plt.style.available)
fig, ax = plt.subplots(figsize=(3.5, 3.5))
 
# Plot line
ax.plot(ms_lbfgs['iter'], ms_lbfgs['misfit']/ms_lbfgs['misfit'].iloc[0], 
        label='L-BFGS', lw=2.5, color='#d62728', zorder=3)
 
# Labels
ax.set_xlabel('Iteration', fontweight='bold')
ax.set_ylabel('Normalized Misfit', fontweight='bold')
 
# Remove top and right spines for a cleaner look
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
 
# Add light grid
ax.grid(True, linestyle='--', alpha=0.6, zorder=0)
 
plt.tight_layout()
plt.savefig('misfit.png', dpi=600, bbox_inches='tight')

Plotting the normalized misfit (divided by its value at iteration 0) makes it easy to compare runs with different datasets or scalings. A healthy anisotropic run drops the normalized misfit to ~22% by 70 iterations.

Traveltime residual distribution

The traveltime residual (synthetic minus observed) summarises how well the model explains the data. Comparing the initial and final distributions shows the reduction in both bias and spread.


import pandas as pd
import matplotlib.pyplot as plt
 
data   = pd.read_csv('src_rec_data_wus.csv')
fwd0   = pd.read_csv('./OUTPUT_FILES_ph_aniso/src_rec_file_forward_PH_000.csv')
fwd39  = pd.read_csv('./OUTPUT_FILES_ph_aniso/src_rec_file_forward_PH_069.csv')
 
tr  = fwd0['tt']  - data['tt']
trf = fwd39['tt'] - data['tt']
 
print(f"Initial: mean = {tr.mean():.4f}, std = {tr.std():.4f}")
print(f"Final:   mean = {trf.mean():.4f}, std = {trf.std():.4f}")
 
fig, ax = plt.subplots(figsize=(4.5, 3.5))
 
ax.hist(tr,  bins=100, alpha=0.7, color='#1f77b4', edgecolor='none', label='Initial model', zorder=2)
ax.hist(trf, bins=100, alpha=0.7, color='#c7414c', edgecolor='none', label='Final model',   zorder=2)
 
ax.legend(frameon=False, loc='upper right')
ax.set_xlim(-15, 15)
ax.set_xlabel('Traveltime residual (s)', fontweight='bold')
ax.set_ylabel('Count', fontweight='bold')
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.grid(True, axis='y', linestyle='--', alpha=0.6, zorder=0)
 
plt.tight_layout()
plt.savefig('residual.png', dpi=600, bbox_inches='tight')
plt.show()