Plot inversion grids

TomoATT use the multiple-grid parameterization to descritize model parameter udpate on multiple coarse mesh grids. This post-processing method decreases the number of unknowns, which enhances the robustness of the inversion. (See Tong et al., 2019 for details)

Setting the inversion grids is essential for inversion. Theoretically, given sufficient data coverage, anomalies whose sizes are smaller than a half of the inversion grid spacing can not be resolved at all, while those whose sizes are greater than two times of the inversion grid spacing can be fully resolved. The anomalies whose sizes are between these two thresholds can be partially resolved.

Thus, choosing appropriate inversion grid spacing, which can sufficiently small to resolve important anomalies and should be sufficiently large to fit the limited data coverage, signficantly influences the imaging result and its reliability.

The inversion grid nodes are output in OUTPUT_path/inversion_grid.txt. We first read the inversion grid nodes


import pygmt
pygmt.config(FONT="16p", IO_SEGMENT_MARKER="<<<")
 
import sys
sys.path.append('../utils')
import functions_for_data as ffd
 
import os
try: 
    os.mkdir('img')
except:
    pass
 
inv_grid_vel, inv_grid_ani = ffd.read_inversion_grid_file("OUTPUT_FILES")
 
print("inversion grid for velocity: ", inv_grid_vel.shape)
print("inversion grid for anisotropy: ", inv_grid_vel.shape)

Then, we plot the grid nodes


Nset    = inv_grid_vel.shape[0]
Ngrid   = inv_grid_vel.shape[1]
 
colorlist = ["green","blue","red","purple","orange","yellow","black","gray","pink","cyan"]
 
# %%
# plot earthquakes and locations
 
fig = pygmt.Figure()
 
pygmt.makecpt(cmap="jet", series=[0, 40], background=True, reverse=True)    # colorbar
 
# -------- horizontal view (x-y) --------
fig.basemap(region=[0,2,0,2], frame=["xa1","ya1","NsWe"], projection="M10c")   # base map
# plot inversion grid
for igrid in range(Nset):
    x = inv_grid_vel[igrid,:,0]
    y = inv_grid_vel[igrid,:,1]
    fig.plot(x=x, y=y, style="c0.1c", fill=colorlist[igrid])
 
# -------- vertical view (x-z) --------
fig.shift_origin(xshift=0, yshift=-3)  
fig.basemap(region=[0,2,0,40], frame=["xa1","ya20+lDepth (km)","NsWe"], projection="X10c/-2c")   # base map
# plot inversion grid
for igrid in range(Nset):
    x = inv_grid_vel[igrid,:,0]
    y = inv_grid_vel[igrid,:,2]
    fig.plot(x=x, y=y, style="c0.1c", fill=colorlist[igrid])
 
 
# -------- vertical view (z-y) --------
fig.shift_origin(xshift=11, yshift=3)  
fig.basemap(region=[0,40,0,2], frame=["xa20+lDepth (km)","ya1","NsWe"], projection="X2c/10c")   # base map
# plot inversion grid
for igrid in range(Nset):
    x = inv_grid_vel[igrid,:,2]
    y = inv_grid_vel[igrid,:,1]
    fig.plot(x=x, y=y, style="c0.1c", fill=colorlist[igrid])
 
 
fig.savefig("img/7_inversion_grid.png")

Here, 5 sets of inversion grids are used in the inversion. Dots in the same color belong to one inversion grid.