Plot models

Here, we demonstrate how to visualize the model file in HDF5 format. First, we load necessary Python modules


import pygmt
pygmt.config(FONT="16p", IO_SEGMENT_MARKER="<<<")
 
from pytomoatt.model import ATTModel
from pytomoatt.data import ATTData
import numpy as np
 
import os
try: 
    os.mkdir('img')
except:
    pass

1. Read model files

Using the PyTomoATT module, we read the initial model and output model as ATTModel objects


# ---------------- read model files ----------------
# file names
init_model_file  = 'input_files/model_init_N61_61_61.h5'        # initial model file
inv_model_file  = 'OUTPUT_FILES/final_model.h5'               # final model file
par_file    = 'input_files/input_params.yaml'                   # parameter file
 
# read initial and final model file
model = ATTModel.read(init_model_file, par_file)
init_model = model.to_xarray()
 
model = ATTModel.read(inv_model_file, par_file)
inv_model = model.to_xarray()

Users can easily check the 3D numpy arrays of

velocity vel;
anti-clock angle with respect of the east direction phi;
magnitude of azimuthal anisotropy epsilon;
original azimuthally anisotropic parameters xi and eta.


# ---------------- access 3D model parameters ----------------
 
# we can access 3D dataset with keys:
# 1. "vel" for velocity
# 2. "phi" fast velocity direction, anti-clock angle w.r.t the east direction. (only available for anisotropic model) 
# 3. "epsilon" for anisotropic magnitude (only available for anisotropic model)
# 4. "xi" and "eta" for anisotropic parameters: xi = epsilon * cos(phi), eta = epsilon * sin(phi)
vel_3d_array = inv_model["vel"]
phi_3d_array = inv_model["phi"]
epsilon_3d_array = inv_model["epsilon"]
xi_3d_array = inv_model["xi"]
eta_3d_array = inv_model["eta"]
 
print("3D array of model parameters. \n vel: ", vel_3d_array.shape, " \n phi: ", phi_3d_array.shape, 
      " \n epsilon: ", epsilon_3d_array.shape, " \n xi: ", xi_3d_array.shape, " \n eta: ", eta_3d_array.shape)

2. Plot horizontal profile of velocity

Here, we plot the perturbation of the output velocity relative to the initial velocity at 20 km depth.

First, we extract the 2D profiles of the initial and output velocity at 20 km depth as vel_init and vel_inv.


# interp vel at depth = 20 km
depth = 20.0
vel_init    = init_model.interp_dep(depth, field='vel')    # vel_init[i,:] are (lon, lat, vel)
vel_inv   = inv_model.interp_dep(depth, field='vel')    
 
print("vel_depth at depth = ", depth, " km. vel_depth:", vel_init.shape, ", (lon, lat, vel)")

Second, we generate the grid object for PyGMT plotting


x = vel_init[:,0];  # longitude
y = vel_init[:,1];  # latitude
value = (vel_inv[:,2] - vel_init[:,2])/vel_init[:,2] * 100    # velocity perturbation relative to the initial model
grid = pygmt.surface(x=x, y=y, z=value, spacing=0.04,region=[0,2,0,2])

Finally, we plot the horizontal profile, and output the figure as img/1_dep_vel.png.


# plot
fig = pygmt.Figure()
fig.basemap(region=[0,2,0,2], frame=["xa1","ya1","+tVelocity perturbation"], projection="M10c")   # base map
pygmt.makecpt(cmap="../utils/svel13_chen.cpt", series=[-20, 20], background=True, reverse=False)    # colorbar
 
fig.grdimage(grid = grid)   # plot figure
fig.text(text="%d km"%(depth), x = 0.2 , y = 0.1, font = "14p,Helvetica-Bold,black", fill = "white") 
 
# colorbar
fig.shift_origin(xshift=0, yshift=-1.5)  
fig.colorbar(frame = ["a20","y+ldlnVp (%)"], position="+e+w4c/0.3c+h") 
 
fig.savefig("img/1_dep_vel.png")

The figure is shown

3. Plot horizontal profile of azimuthal anisotropy

Here, we plot the azimuthal anisotropy of the output model at 20 km depth.

If the maximum magnitude of anisotropy is 0 across this profile (isotropic structure), an error will be reported, as the direction of anisotropy is not defined.

First, First, we extract the 2D profiles of anisotrop magnitude as epsilon_inv


# ---------------- 2D depth profile of azimuthal anisotropy ----------------
 
# interp magnitude of anisotropy at depth = 20 km
depth = 20.0
epsilon_inv   = inv_model.interp_dep(depth, field='epsilon')    # epsilon_inv[i,:] are (lon, lat, epsilon)
 
print("epsilon_inv at depth = ", depth, " km. epsilon_inv:", epsilon_inv.shape, ", (lon, lat, epsilon)")

Second, we build an array to plot fast velocity direction in PyGMT


# generate fast velocity direction (anisotropic arrow)
samp_interval = 3
length = 10
width = 0.1
ani_thd = 0.02
ani_phi = inv_model.interp_dep(depth, field='phi', samp_interval=samp_interval)
ani_epsilon = inv_model.interp_dep(depth, field='epsilon', samp_interval=samp_interval)
ani_arrow = np.hstack([ani_phi, ani_epsilon[:,2].reshape(-1, 1)*length, np.ones((ani_epsilon.shape[0],1))*width])   # lon, lat, angle, length, width
idx = np.where(ani_epsilon[:,2] > ani_thd)
ani_arrow = ani_arrow[idx[0],:]
 
print("ani_arrow at depth = ", depth, " km. ani_arrow:", ani_arrow.shape, ", (lon, lat, angle, length, width)")

samp_interval controls the density to plot fast velocity direction (FVD) bars. Based on the forward mesh grid, we select one node for plotting every samp_interval nodes in each direction;
length: length scaling of the FVD bars;
width: width of the FVD bars
ani_thd: FVD bars whose anisotropic magnitude is less than ani_thd are ignore.

Third, we generate the grid object of anisotropic magnitude for PyGMT plotting


x = epsilon_inv[:,0];  # longitude
y = epsilon_inv[:,1];  # latitude
value = epsilon_inv[:,2] # magnitude of anisotropy
grid = pygmt.surface(x=x, y=y, z=value, spacing=0.04,region=[0,2,0,2])

Finally, we plot the horizontal profile, and output the figure as img/1_dep_ani.png.


# plot
fig = pygmt.Figure()
fig.basemap(region=[0,2,0,2], frame=["xa1","ya1","+tAzimuthal Anisotropy"], projection="M10c")   # base map
pygmt.makecpt(cmap="cool", series=[0, 0.1], background=True, reverse=False)    # colorbar
 
fig.grdimage(grid = grid)   # plot magnitude of anisotropy
fig.plot(ani_arrow, style='j', fill='yellow1', pen='0.5p,black')    # plot fast velocity direction
 
fig.text(text="%d km"%(depth), x = 0.2 , y = 0.1, font = "14p,Helvetica-Bold,black", fill = "white") 
 
# colorbar
fig.shift_origin(xshift=0, yshift=-1.5)  
fig.colorbar(frame = ["a0.1","y+lAnisotropy"], position="+e+w4c/0.3c+h") 
 
fig.savefig("img/1_dep_ani.png")

The figure is shown

3. Plot vertical profile of velocity

First, we extract the 2D vertical profiles of the initial and output velocity from (0, 0.75) to (2, 0.75) as vel_init_sec and vel_inv_sec.


# ---------------- 2D vertical profile of velocity perturbation ----------------
 
# interp vel from [0,0.75] in lon-lat to [2,0.75] in lon-lat, gap = 1 km
start = [0,0.75]; end = [2,0.75]; gap = 1
vel_init_sec    = init_model.interp_sec(start, end, field='vel', val = gap)      # vel_init_sec[i,:] are (lon, lat, dis, dep, vel)
vel_inv_sec   = inv_model.interp_sec(start, end, field='vel', val = gap)       
 
print("vel_init_sec:", vel_init_sec.shape, ", (lon, lat, distance, depth, vel)")

Second, we generate the grid object of anisotropic magnitude for PyGMT plotting


x = vel_init_sec[:,0];  # longitude
y = vel_init_sec[:,3];  # depth
value = (vel_inv_sec[:,4] - vel_init_sec[:,4])/vel_init_sec[:,4] * 100    # velocity perturbation relative to the initial model
grid = pygmt.surface(x=x, y=y, z=value, spacing="0.04/1",region=[0,2,0,40])

Finally, we plot the vertical profile, and output the figure as img/1_sec_vel.png.


# plot
fig = pygmt.Figure()
fig.basemap(region=[0,2,0,40], frame=["xa1+lLongitude","ya20+lDepth (km)","+tVelocity perturbation"], projection="X10c/-4c")   # base map
pygmt.makecpt(cmap="../utils/svel13_chen.cpt", series=[-20, 20], background=True, reverse=False)    # colorbar
 
fig.grdimage(grid = grid)   # plot figure
fig.text(text="A", x = 0.1 , y = 5, font = "14p,Helvetica-Bold,black", fill = "white") 
fig.text(text="A@+'@+", x = 1.9 , y = 5, font = "14p,Helvetica-Bold,black", fill = "white") 
 
# colorbar
fig.shift_origin(xshift=0, yshift=-2)  
fig.colorbar(frame = ["a20","y+ldlnVp (%)"], position="+e+w4c/0.3c+h") 
 
fig.savefig("img/1_sec_vel.png")

The figure is shown