Result visualization

Make sure the Python module Pygmt is loaded.

This part plot the output files for result visualization. Users can directly run the following codes to plot the results automatically.


python plot_output.py

Figures are output as

img/model_setting.png for earthquake location (dots), station location (triangles), velocity perturbations of the initial and true models, and azimuthal anisotropy (yellow bar).

img/model_inv.png for the velocity perturbation and azimuthal anisotropy of the inversion results.

Here we demonstrate the plotting script step-by-step

1. Load necessary Python modules


import pygmt
import os
from pytomoatt.model import ATTModel
from pytomoatt.data import ATTData
from pytomoatt.src_rec import SrcRec
import numpy as np

2. Read the initial model and the checkerboard model


Ngrid = [61,61,61]
data_file = '2_models/model_init_N%d_%d_%d.h5'%(Ngrid[0],Ngrid[1],Ngrid[2])
par_file = '3_input_params/input_params_signal.yaml'  
model = ATTModel.read(data_file, par_file)
initial_model = model.to_xarray()
 
data_file = '2_models/model_ckb_N%d_%d_%d.h5'%(Ngrid[0],Ngrid[1],Ngrid[2])
model = ATTModel.read(data_file, par_file)
ckb_model = model.to_xarray()

3. Slice profiles

Choose a horizontal profile at 10 km depth and a vertical profile from (1.25, 0) to (1.25, 2) in longitude and latitude.


depth = 10.0
vel_init = initial_model.interp_dep(depth, field='vel')
vel_ckb = ckb_model.interp_dep(depth, field='vel')
 
start = [1.25,0]; end = [1.25,2]
vel_init_sec = initial_model.interp_sec(start, end, field='vel', val = 1)
vel_ckb_sec = ckb_model.interp_sec(start, end, field='vel', val = 1)

4. Generate fast velocity direction arrows for the checkboard model


samp_interval = 3   # control the density of fast velocity direction arrows
length = 7
width = 0.1
ani_thd = 0.02  # arrows with magnitude smaller ani_thd is ignored
 
ani_ckb_phi = ckb_model.interp_dep(depth, field='phi', samp_interval=samp_interval)
ani_ckb_epsilon = ckb_model.interp_dep(depth, field='epsilon', samp_interval=samp_interval)
ani_ckb = np.hstack([ani_ckb_phi, ani_ckb_epsilon[:,2].reshape(-1, 1)*length, np.ones((ani_ckb_epsilon.shape[0],1))*width])   # lon, lat, angle, length, width
idx = np.where(ani_ckb_epsilon[:,2] > ani_thd)
ani_ckb = ani_ckb[idx[0],:]

The magnitude smaller 0.02 are ignored by


ani_thd = 0.02  # arrows with magnitude smaller ani_thd is ignored

5. Read earthquakes and stations


sr = SrcRec.read('1_src_rec_files/src_rec_config.dat')
station = sr.receivers[['stlo','stla','stel']].values.T
true_loc = sr.sources[['evlo','evla','evdp']].values.T
earthquake = true_loc
 
# categorize earthquakes
ev_idx1 = []
ev_idx2 = []
ev_idx3 = []
for i in range(earthquake.shape[1]):
    dep = earthquake[2,i]
    if dep < 15:
        ev_idx1.append(i)
    elif dep < 25:
        ev_idx2.append(i)
    elif dep < 35:
        ev_idx3.append(i)

6. Plot perturbation of the initial and true models relative to the initial model


fig = pygmt.Figure()
 
region = [0,2,0,2]
frame  = ["xa1","ya1","nSWe"]
projection = "M10c"
spacing = 0.04
 
vel_range = 20
 
# -------------- initial model and earthquake location --------------
fig.basemap(region=region, frame=["xa1","ya1","+tInitial model and locations"], projection=projection)
# velocity perturbation
pygmt.makecpt(cmap="../utils/svel13_chen.cpt", series=[-vel_range, vel_range], background=True, reverse=False)
x = vel_init[:,0]; y = vel_init[:,1]; value = (vel_init[:,2] - vel_init[:,2])/vel_init[:,2] * 100
grid = pygmt.surface(x=x, y=y, z=value, spacing=spacing,region=region)
fig.grdimage(grid = grid)
# earthquakes
fig.plot(x = true_loc[0,ev_idx1], y = true_loc[1,ev_idx1], style = "c0.1c", fill = "red")
fig.plot(x = true_loc[0,ev_idx2], y = true_loc[1,ev_idx2], style = "c0.1c", fill = "green")
fig.plot(x = true_loc[0,ev_idx3], y = true_loc[1,ev_idx3], style = "c0.1c", fill = "black")
 
# stations
fig.plot(x = station[0,:], y = station[1,:], style = "t0.4c", fill = "blue", pen = "black", label = "Station")
 
# # anisotropic arrow
# fig.plot(ani_ckb, style='j', fill='yellow1', pen='0.5p,black')
 
fig.shift_origin(xshift=11)
 
fig.basemap(region=[0,40,0,2], frame=["xa20+lDepth (km)","ya1","Nswe"], projection="X2c/10c")
x = vel_init_sec[:,3]; y = vel_init_sec[:,1]; value = (vel_init_sec[:,4] - vel_init_sec[:,4])/vel_init_sec[:,4] * 100
grid = pygmt.surface(x=x, y=y, z=value, spacing="1/0.04",region=[0,40,0,2])
fig.grdimage(grid = grid)
 
# earthquakes
fig.plot(x = true_loc[2,ev_idx1], y = true_loc[1,ev_idx1], style = "c0.1c", fill = "red")
fig.plot(x = true_loc[2,ev_idx2], y = true_loc[1,ev_idx2], style = "c0.1c", fill = "green")
fig.plot(x = true_loc[2,ev_idx3], y = true_loc[1,ev_idx3], style = "c0.1c", fill = "black")
 
 
fig.shift_origin(xshift=4)
 
 
 
 
# -------------- checkerboard model --------------
fig.basemap(region=region, frame=["xa1","ya1","+tTrue model and locations"], projection=projection)
# velocity perturbation
pygmt.makecpt(cmap="../utils/svel13_chen.cpt", series=[-vel_range, vel_range], background=True, reverse=False)
x = vel_ckb[:,0]; y = vel_ckb[:,1]; value = (vel_ckb[:,2] - vel_init[:,2])/vel_init[:,2] * 100
grid = pygmt.surface(x=x, y=y, z=value, spacing=spacing,region=region)
fig.grdimage(grid = grid)
# earthquakes
fig.plot(x = earthquake[0,ev_idx1], y = earthquake[1,ev_idx1], style = "c0.1c", fill = "red")
fig.plot(x = earthquake[0,ev_idx2], y = earthquake[1,ev_idx2], style = "c0.1c", fill = "green")
fig.plot(x = earthquake[0,ev_idx3], y = earthquake[1,ev_idx3], style = "c0.1c", fill = "black")
 
# stations
# fig.plot(x = loc_st[0,:], y = loc_st[1,:], style = "t0.4c", fill = "blue", pen = "black", label = "Station")
 
# anisotropic arrow
fig.plot(ani_ckb, style='j', fill='yellow1', pen='0.5p,black')
 
fig.shift_origin(xshift=11)
 
fig.basemap(region=[0,40,0,2], frame=["xa20+lDepth (km)","ya1","Nswe"], projection="X2c/10c")
x = vel_ckb_sec[:,3]; y = vel_ckb_sec[:,1]; value = (vel_ckb_sec[:,4] - vel_init_sec[:,4])/vel_init_sec[:,4] * 100
grid = pygmt.surface(x=x, y=y, z=value, spacing="1/0.04",region=[0,40,0,2])
fig.grdimage(grid = grid)
 
# earthquakes
fig.plot(x = earthquake[2,ev_idx1], y = earthquake[1,ev_idx1], style = "c0.1c", fill = "red")
fig.plot(x = earthquake[2,ev_idx2], y = earthquake[1,ev_idx2], style = "c0.1c", fill = "green")
fig.plot(x = earthquake[2,ev_idx3], y = earthquake[1,ev_idx3], style = "c0.1c", fill = "black")
 
 
# ------------------- colorbar -------------------
fig.shift_origin(xshift=-11, yshift=-1.5)
fig.colorbar(frame = ["a%f"%(vel_range),"x+ldlnVp (%)"], position="+e+w4c/0.3c+h") 
 
fig.shift_origin(xshift=6, yshift=-1)
fig.basemap(region=[0,1,0,1], frame=["wesn"], projection="X6c/1.5c")
ani = [
    [0.2, 0.6, 45, 0.02*length, width],    # lon, lat, phi, epsilon, size
    [0.5, 0.6, 45, 0.05*length, width],
    [0.8, 0.6, 45, 0.10*length, width],
    ]
fig.plot(ani, style='j', fill='yellow1', pen='0.5p,black')
fig.text(text=["0.02", "0.05", "0.10"], x=[0.2,0.5,0.8], y=[0.2]*3, font="16p,Helvetica", justify="CM")
fig.shift_origin(xshift= 11, yshift=2.5)
 
fig.show()
fig.savefig('img/model_setting.png', dpi=300)

7. Plot output model perturbations

Similarly, we can plot the model perturbation for the imaging results


# plot the inversion result
 
# read models
tag = "inv"
data_file = "OUTPUT_FILES/OUTPUT_FILES_%s/final_model.h5"%(tag)
model = ATTModel.read(data_file, par_file)
inv_model = model.to_xarray()
vel_inv = inv_model.interp_dep(depth, field='vel') # lon = [:,0], lat = [:,1], vel = [:,2]  
x = vel_inv[:,0]; y = vel_inv[:,1]; value = (vel_inv[:,2] - vel_init[:,2])/vel_init[:,2] * 100
vel_inv_sec = inv_model.interp_sec(start, end, field='vel', val = 1)
x_sec = vel_inv_sec[:,3]; y_sec = vel_inv_sec[:,1]; value_sec = (vel_inv_sec[:,4] - vel_init_sec[:,4])/vel_init_sec[:,4] * 100
 
ani_inv_phi = inv_model.interp_dep(depth, field='phi', samp_interval=samp_interval)
ani_inv_epsilon = inv_model.interp_dep(depth, field='epsilon', samp_interval=samp_interval)
ani_inv = np.hstack([ani_inv_phi, ani_inv_epsilon[:,2].reshape(-1, 1)*length, np.ones((ani_inv_epsilon.shape[0],1))*width])   # lon, lat, angle, length, width
idx = np.where(ani_inv_epsilon[:,2] > ani_thd)
ani_inv = ani_inv[idx[0],:]
 
# plot the inversion result
 
fig = pygmt.Figure()
 
region = [0,2,0,2]
frame  = ["xa1","ya1","+tInversion results"]
projection = "M10c"
spacing = 0.04
 
vel_range = 20
 
# -------------- checkerboard model --------------
fig.basemap(region=region, frame=frame, projection=projection)
# velocity perturbation
pygmt.makecpt(cmap="../utils/svel13_chen.cpt", series=[-vel_range, vel_range], background=True, reverse=False)
x = vel_inv[:,0]; y = vel_inv[:,1]; value = (vel_inv[:,2] - vel_init[:,2])/vel_init[:,2] * 100
grid = pygmt.surface(x=x, y=y, z=value, spacing=spacing,region=region)
fig.grdimage(grid = grid)
 
# anisotropic arrow
fig.plot(ani_inv, style='j', fill='yellow1', pen='0.5p,black')
 
fig.shift_origin(xshift=11)
 
fig.basemap(region=[0,40,0,2], frame=["xa20+lDepth (km)","ya1","Nswe"], projection="X2c/10c")
x = vel_inv_sec[:,3]; y = vel_inv_sec[:,1]; value = (vel_inv_sec[:,4] - vel_init_sec[:,4])/vel_init_sec[:,4] * 100
grid = pygmt.surface(x=x, y=y, z=value, spacing="1/0.04",region=[0,40,0,2])
fig.grdimage(grid = grid)
 
# ------------------- colorbar -------------------
fig.shift_origin(xshift=-11, yshift=-1.5)
fig.colorbar(frame = ["a%f"%(vel_range),"x+ldlnVp (%)"], position="+e+w4c/0.3c+h") 
 
fig.shift_origin(xshift=6, yshift=-1)
fig.basemap(region=[0,1,0,1], frame=["wesn"], projection="X6c/1.5c")
ani = [
    [0.2, 0.6, 45, 0.02*length, width],    # lon, lat, phi, epsilon, size
    [0.5, 0.6, 45, 0.05*length, width],
    [0.8, 0.6, 45, 0.10*length, width],
    ]
fig.plot(ani, style='j', fill='yellow1', pen='0.5p,black')
fig.text(text=["0.02", "0.05", "0.10"], x=[0.2,0.5,0.8], y=[0.2]*3, font="16p,Helvetica", justify="CM")
fig.shift_origin(xshift= 11, yshift=2.5)
 
 
fig.show()
fig.savefig('img/model_%s.png'%(tag), dpi=300)