Plot models

Here, we demonstrate how to visualize the traveltime field and adjoint field. 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

1. Read model files

Using the PyTomoATT module, we read the time field and adjoint_field of the source JC05 at the 7-th iteration as time_field and adjoint_field, repsectively.

/src_${src_name}/time_field_inv_${Nstep} and /src_${src_name}/adjoint_field_inv_${Nstep} are the paths of traveltime and adjoint fields in the output HDF5 file out_data_sim_group_XX.h5.


# ---------------- read files ----------------
src_name = 'JC05'
Nstep = "0007"
 
# file names
data_file   = 'OUTPUT_FILES/out_data_sim_group_5.h5'        # data file
par_file    = 'input_files/input_params.yaml'               # parameter file
grid_file   = 'OUTPUT_FILES/out_data_grid.h5'               # grid file
group       = 'src_%s'%(src_name)                           # src_${src_name}     
 
# read traveltime field
dataset_time     = 'time_field_inv_%s'%(Nstep)                   # time_field_inv_${Nstep}
data = ATTData.read(data_file, par_file, grid_file, group, dataset_time)
time_field = data.to_xarray()
 
# read adjoint field
dataset_adjoint  = 'adjoint_field_inv_%s'%(Nstep)                # adjoint_field_inv_${Nstep}
data = ATTData.read(data_file, par_file, grid_file, group, dataset_adjoint)
adjoint_field = data.to_xarray()

Since source parallelization is applied, fileds of each source is store in a specific file, named as out_data_sim_group_XX.h5. Users can check which sources are contained in the specific out_data_sim_group_XX.h5 file, e.g.,


h5ls OUTPUT_FILES/out_data_sim_group_5.h5

The output log is


src_JC05                 Group
src_JC13                 Group
src_JC21                 Group

src_${source_name} means that the source is contained in the file out_data_sim_group_0.h5.

Furthermore, using


h5ls OUTPUT_FILES/out_data_sim_group_5.h5/src_JC05

yields


adjoint_field_inv_0000   Dataset {226981}
adjoint_field_inv_0001   Dataset {226981}
adjoint_field_inv_0002   Dataset {226981}
adjoint_field_inv_0003   Dataset {226981}
adjoint_field_inv_0004   Dataset {226981}
adjoint_field_inv_0005   Dataset {226981}
adjoint_field_inv_0006   Dataset {226981}
adjoint_field_inv_0007   Dataset {226981}
adjoint_field_inv_0008   Dataset {226981}
adjoint_field_inv_0009   Dataset {226981}
time_field_inv_0000      Dataset {226981}
time_field_inv_0001      Dataset {226981}
time_field_inv_0002      Dataset {226981}
time_field_inv_0003      Dataset {226981}
time_field_inv_0004      Dataset {226981}
time_field_inv_0005      Dataset {226981}
time_field_inv_0006      Dataset {226981}
time_field_inv_0007      Dataset {226981}
time_field_inv_0008      Dataset {226981}
time_field_inv_0009      Dataset {226981}

It contains traveltime fields and adjoint fields at all iterations.

Users can easily check the 3D numpy arrays of

traveltime field time_3d_array;
adjoint time field adjoint_3d_array; and 1D numpy arrays of
depth coordinates dep_1d_array;
latitude cooridinates lat_1d_array;
longitude coordinates lon_1d_array.


# ---------------- access 3D time field and adjoint field ----------------
# we can access 3D dataset:
 
dep_1d_array = time_field["dep"]
lat_1d_array = time_field["lat"]
lon_1d_array = time_field["lon"]
 
print("3D array of coordinates. \n dep: ", dep_1d_array.shape, " \n lat: ", lat_1d_array.shape, " \n lon: ", lon_1d_array.shape)
 
time_3d_array = time_field[dataset_time]
adjoint_3d_array = adjoint_field[dataset_adjoint]
 
print("3D array of fields. \n time: ", time_3d_array.shape, " \n adjoint: ", adjoint_3d_array.shape)

2. Plot horizontal profile of the traveltime field

Here, we plot the horizontal profile of traveltime field at 20 km depth.

First, we extract the 2D profiles of the traveltime field at 20 km depth as time.


# interp vel at depth = 20 km
depth = 20.0
time    = time_field.interp_dep(depth, field=dataset_time)    # time[i,:] are (lon, lat, vel)
 
print("time at depth = ", depth, " km. time:", time.shape, ", (lon, lat, time)")

Second, we generate the grid object for PyGMT plotting


x = time[:,0];      # longitude
y = time[:,1];      # latitude
value = time[:,2]   # traveltime
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/2_dep_time.png.


# plot
fig = pygmt.Figure()
fig.basemap(region=[0,2,0,2], frame=["xa1","ya1","+tTraveltime"], projection="M10c")   # base map
pygmt.makecpt(cmap="jet", series=[0, 30], background=True, reverse=True)    # colorbar
 
fig.grdimage(grid = grid)   # plot figure
fig.contour(x=x, y=y, z=value, levels=5, pen="1.5p,white")   # contour
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+lTraveltime (s)"], position="+e+w4c/0.3c+h") 
 
fig.savefig("img/2_dep_time.png")

The figure is shown as below. The time contour is not a perfect circle due to velocity perturbation and azimuthal anisotrop (see models).

2. Plot vertical profile of the traveltime field

Here, we plot the vertical profile of traveltime field

First, we extract the 2D profiles of the traveltime field from (0,0.6) to (2,0.6). as time_sec.


# 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
time_sec    = time_field.interp_sec(start, end, field=dataset_time, val = gap)      # time_sec[i,:] are (lon, lat, dis, dep, time)
 
print("time_sec:", time_sec.shape, ", (lon, lat, distance, depth, time)")

Second, we generate the grid object for PyGMT plotting


x = time_sec[:,0];  # longitude
y = time_sec[:,3];  # depth
value = time_sec[:,4] # traveltime
grid = pygmt.surface(x=x, y=y, z=value, spacing="0.04/1",region=[0,2,0,40])

Finally, we plot the horizontal profile, and output the figure as img/2_sec_time.png.


# plot
fig = pygmt.Figure()
fig.basemap(region=[0,2,0,40], frame=["xa1+lLongitude","ya20+lDepth (km)","+tTraveltime"], projection="X10c/-2c")   # base map
pygmt.makecpt(cmap="jet", series=[0, 30], background=True, reverse=True)    # colorbar
 
fig.grdimage(grid = grid)   # plot figure
fig.contour(x=x, y=y, z=value, levels=5, pen="1.5p,white")   # contour
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+lTraveltime (s)"], position="+e+w4c/0.3c+h") 
 
fig.savefig("img/2_sec_time.png")

The figure is shown as below.

3. Plot horizontal profile of the adjoint field

Here, we plot the horizontal profile of adjoint field at 20 km depth.

First, we extract the 2D profiles of the adjoint field at 20 km depth as adjoint.


# interp vel at depth = 20 km
depth = 20.0
adjoint = adjoint_field.interp_dep(depth, field=dataset_adjoint)    
 
print("time at depth = ", depth, " km. time:", time.shape, ", (lon, lat, time)")

Second, we generate the grid object for PyGMT plotting


x = time[:,0];      # longitude
y = time[:,1];      # latitude
value = adjoint[:,2]   # traveltime
value = value/np.nanmax(np.abs(value))
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/2_dep_adjoint.png.


# plot
fig = pygmt.Figure()
fig.basemap(region=[0,2,0,2], frame=["xa1","ya1","+tAdjoint field"], projection="M10c")   # base map
pygmt.makecpt(cmap="jet", series=[-0.5, 0.5], background=True, reverse=False)    # colorbar
 
fig.grdimage(grid = grid)   # plot figure
fig.contour(x=x, y=y, z=time[:,2], levels=5, pen="1.5p,white")   # contour
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+lAdjoint field"], position="+e+w4c/0.3c+h") 
 
fig.savefig("img/2_dep_adjoint.png")

The figure is shown as below. The kernel looks a little bit thick due to the smoothing of coarse discretized mesh grid.

4. Plot vertical profile of the adjoint field

Here, we plot the vertical profile of adjoint field.

First, we extract the 2D profiles of the adjoint field from (0,0.6) to (2,0.6). as adjoint_sec.


# 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
adjoint_sec   = adjoint_field.interp_sec(start, end, field=dataset_adjoint, val = gap)       
 
print("adjoint_sec:", time_sec.shape, ", (lon, lat, distance, depth, adjoint)")

Second, we generate the grid object for PyGMT plotting


x = adjoint_sec[:,0];  # longitude
y = adjoint_sec[:,3];  # depth
value = adjoint_sec[:,4] # traveltime
value = value/np.nanmax(np.abs(value))
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/2_sec_adjoint.png.


# plot
fig = pygmt.Figure()
fig.basemap(region=[0,2,0,40], frame=["xa1+lLongitude","ya20+lDepth (km)","+tAdjoint field"], projection="X10c/-2c")   # base map
pygmt.makecpt(cmap="jet", series=[-0.5, 0.5], background=True, reverse=False)    # colorbar
 
fig.grdimage(grid = grid)   # plot figure
fig.contour(x=x, y=y, z=time_sec[:,4], levels=5, pen="1.5p,white")   # contour
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+lAdjoint"], position="+e+w4c/0.3c+h") 
 
fig.savefig("img/2_sec_adjoint.png")

The figure is shown as below.