Plot sensitivity kernels

Here, we demonstrate how to visualize the sensitivity kernel. First, we load necessary Python modules


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

We choose to plot the kernel at 7-th iteration. Related file names are


# ---------------- read files ----------------
Nstep = "0007"
kernel_list = {}    # dictionary to store all the kernels
 
# file names
data_file   = 'OUTPUT_FILES/out_data_sim_group_0.h5'        # data file
par_file    = 'input_files/input_params.yaml'               # parameter file
grid_file   = 'OUTPUT_FILES/out_data_grid.h5'               # grid file
group       = 'model'

We will add all types of kernels to the dictionary kernel_list in the following subsections.

1. Original misfit kernel without modifications

We extract the misfit kernel without any modification, that is, the Frechet derivatives of the objective function with respect to unknown parameters.

The misfit kernel with respect to slowness: kernel_list['Ks_inv_0007'];
The misfit kernel with respect to xi: kernel_list['Kxi_inv_0007'];
The misfit kernel with respect to eta: kernel_list['Keta_inv_0007'].


# (Option 1) read original sensitivity kernel 
# Ks: kernel w.r.t. slowness at the 7-th iteration
dataset     = 'Ks_inv_%s'%(Nstep)                      
data = ATTData.read(data_file, par_file, grid_file, group, dataset)
kernel_list[dataset] = data.to_xarray()
 
# Kxi: kernel w.r.t. xi (anisotropic parameter) at the 7-th iteration
dataset     = 'Kxi_inv_%s'%(Nstep)                      
data = ATTData.read(data_file, par_file, grid_file, group, dataset)
kernel_list[dataset] = data.to_xarray()
 
# Keta: kernel w.r.t. eta (anisotropic parameter) at the 7-th iteration
dataset     = 'Keta_inv_%s'%(Nstep)                      
data = ATTData.read(data_file, par_file, grid_file, group, dataset)
kernel_list[dataset] = data.to_xarray()

2. Kernel density

We extract the misfit kernel density, which is obtained by setting all adjoint source to 1.0. The kernel density reflects the distribution of sensitivity. It can be used to normalize the original kernel for higher convergent rates (see Section 3.3 in Chen et al., 2025 for details).

The kernel density with respect to slowness: kernel_list['Ks_density_inv_0007'];
The kernel density with respect to xi: kernel_list['Kxi_density_inv_0007'];
The kernel density with respect to eta: kernel_list['Keta_density_inv_0007'].


# (Option 2) read kernel_density
# Ks_den: kernel density w.r.t. slowness at the 7-th iteration
dataset     = 'Ks_density_inv_%s'%(Nstep)                      
data = ATTData.read(data_file, par_file, grid_file, group, dataset)
kernel_list[dataset] = data.to_xarray()
 
# Kxi_den: kernel density w.r.t. xi (anisotropic parameter) at the 7-th iteration
dataset     = 'Kxi_density_inv_%s'%(Nstep)                      
data = ATTData.read(data_file, par_file, grid_file, group, dataset)
kernel_list[dataset] = data.to_xarray()
 
# Keta_den: kernel density w.r.t. eta (anisotropic parameter) at the 7-th iteration
dataset     = 'Keta_density_inv_%s'%(Nstep)                      
data = ATTData.read(data_file, par_file, grid_file, group, dataset)
kernel_list[dataset] = data.to_xarray()

3. Smoothed kernel

We extract the smoothed kernels, which are obtained by applying multiple-grid parameterization and kernel density normalization to the original misfit kernel (see Section 3.3 in Chen et al., 2025 for details).

When using the gradient descent method, the model update at this step equals step_length * smoothed_kernel / max(abs(smoothed_kernel))


# (Option 3) read normalized kernel smoothed by multigrid parameterization
# Ks_update: smoothed normalized kernel w.r.t. slowness at the 7-th iteration
dataset     = 'Ks_update_inv_%s'%(Nstep)                      
data = ATTData.read(data_file, par_file, grid_file, group, dataset)
kernel_list[dataset] = data.to_xarray()
 
# Kxi_update: smoothed normalized kernel w.r.t. xi (anisotropic parameter) at the 7-th iteration
dataset     = 'Kxi_update_inv_%s'%(Nstep)                      
data = ATTData.read(data_file, par_file, grid_file, group, dataset)
kernel_list[dataset] = data.to_xarray()
 
# Keta_update: smoothed normalized kernel w.r.t. eta (anisotropic parameter) at the 7-th iteration
dataset     = 'Keta_update_inv_%s'%(Nstep)                      
data = ATTData.read(data_file, par_file, grid_file, group, dataset)
kernel_list[dataset] = data.to_xarray()

4. Plot horizontal profiles of kernels

Users can easily check the 3D kernel array and the coordinates in three directions


# ---------------- access 3D array ----------------
# we can access 3D dataset:
 
dep_1d_array = kernel_list['Ks_inv_0007']["dep"]
lat_1d_array = kernel_list['Ks_inv_0007']["lat"]
lon_1d_array = kernel_list['Ks_inv_0007']["lon"]
 
print("3D array of coordinates. \n dep: ", dep_1d_array.shape, " \n lat: ", lat_1d_array.shape, " \n lon: ", lon_1d_array.shape)
 
array = kernel_list['Ks_inv_0007']["Ks_inv_0007"]
 
print("3D array of kernel. \n Ks: ", array.shape)

The horizontal and vertical profiles of all kernels are plot via


# %%
# ---------------- 2D depth profile of kernels ----------------
 
for dataset in kernel_list:
 
    # interp vel at depth = 20 km
    depth = 20.0
    kernel = kernel_list[dataset].interp_dep(depth, field=dataset)    # kernel[i,:] are (lon, lat, kernel)
 
    print("kernel at depth = ", depth, " km. kernel:", kernel.shape, ", (lon, lat, kernel)")
 
    # plot
    fig = pygmt.Figure()
    fig.basemap(region=[0,2,0,2], frame=["xa1","ya1","+t%s"%(dataset)], projection="M10c")   # base map
    pygmt.makecpt(cmap="jet", series=[-0.5, 0.5], background=True, reverse=True)    # colorbar
 
    x = kernel[:,0];      # longitude
    y = kernel[:,1];      # latitude
    value = kernel[:,2]/np.nanmax(np.abs(kernel[:,2]))   # traveltime
    grid = pygmt.surface(x=x, y=y, z=value, spacing=0.04,region=[0,2,0,2])
 
    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 = ["a0.5","y+l%s"%(dataset)], position="+e+w4c/0.3c+h") 
 
    fig.savefig("img/3_dep_%s.png"%(dataset))

4. Plot vertical profiles of kernels


# %%
# ---------------- 2D vertical profile of kernels ----------------
 
for dataset in kernel_list:
 
    # interp from [0,0.6] in lon-lat to [2,0.6] in lon-lat, gap = 1 km
    start = [0,0.6]; end = [2,0.6]; gap = 1
    kernel_sec    = kernel_list[dataset].interp_sec(start, end, field=dataset, val = gap)      # kernel_sec[i,:] are (lon, lat, dis, dep, kernel)
 
    print("kernel_sec:", kernel_sec.shape, ", (lon, lat, distance, depth, kernel)")
 
    # plot
    fig = pygmt.Figure()
    fig.basemap(region=[0,2,0,40], frame=["xa1+lLongitude","ya20+lDepth (km)","+t%s"%(dataset)], projection="X10c/-2c")   # base map
    pygmt.makecpt(cmap="jet", series=[-0.5, 0.5], background=True, reverse=True)    # colorbar
 
    x = kernel_sec[:,0];  # longitude
    y = kernel_sec[:,3];  # depth
    value = kernel_sec[:,4]/np.nanmax(np.abs(kernel_sec[:,4])) # traveltime
    grid = pygmt.surface(x=x, y=y, z=value, spacing="0.04/1",region=[0,2,0,40])
 
    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 = ["a0.5","y+l%s"%(dataset)], position="+e+w4c/0.3c+h") 
 
    fig.savefig("img/3_sec_%s.png"%(dataset))