Generate HDF5 model files

Given parameter setting file in YAML format input_params/input_params.yaml

users can generate a series of HDF5 model files by


python 1_generate_models.py

The model files are output to the folder models.

Constant velocity models models/constant_velocity_N61_51_101_loop.h5 and models/constant_velocity_N61_51_101_PyTomoATT.h5. These two models are the same, but are generated by using simple loop and PyTomoATT module, respectively.
Linear velocity models models/linear_velocity_N61_51_101_loop.h5 and models/linear_velocity_N61_51_101_PyTomoATT.h5.
Checkerboard models models/linear_velocity_ckb_N61_51_101_loop.h5 and models/linear_velocity_ckb_N61_51_101_PyTomoATT.h5.
Checkerboard models with varying checker sizes models/linear_velocity_ckb_flex_N61_51_101.h5. This model is built by using simple loop.

Here we will demonstrate how to build these models step-by-step. First, please load necessary modules in Python Script:


from pytomoatt.model import ATTModel
from pytomoatt.checkerboard import Checker
import os
import numpy as np
import h5py

1. Read YAML file for mesh grid infomation

We read the parameter setting file to get the information of forward mesh grid:

the grid number in depth, latitude, and longtitude n_rtp.
1D depth array am_depths with the size of n_rtp[0];
1D latitude array am_latitudes with the size of n_rtp[1];
1D longitude array am_longitudes with the size of n_rtp[2];

In addition, we can initialize a model object att_model via PyTomoATT by att_model = ATTModel(par_file)


output_path = "models"
os.makedirs(output_path, exist_ok=True)
 
par_file = "input_params/input_params.yaml"
 
att_model = ATTModel(par_file)
 
n_rtp          = att_model.n_rtp       # grid node numbers in r (depth), t (longtitude), p (latitude) directions
am_depths      = att_model.depths      # depths in km
am_latitudes   = att_model.latitudes   # latitudes in degrees
am_longitudes  = att_model.longitudes  # longitudes in degrees
 
print("grid node numbers (N_r, N_t, N_p):", n_rtp)
print("depths (km):", am_depths)
print("latitudes (degree):", am_latitudes)
print("longitudes (degree):", am_longitudes)

The model file in HDF5 format for TomoATT consists of three dataset. We can use


h5ls model_name

to show the structure


eta                      Dataset {61, 51, 101}
vel                      Dataset {61, 51, 101}
xi                       Dataset {61, 51, 101}

Only three 3D arrays with the size of n_rtp are required. eta and xi represent azimuthally anisotropic parameters. vel represents 3D seismic velocity values.

Thus, our aim is to generate such a HDF5 file. Some examples are provided as follwos.

1. Build constant velocity model

Using PyTomoATT

We can assign a constant velocity and zero anisotropy to the model object att_model


# case 1. ---------- generate a constant velocity model using PyTomoATT module -------------
 
# set the velocity model to a constant value
constant_v = 6.0 # constant velocity (km/s)
att_model.vel[:,:,:] = constant_v
att_model.xi[:,:,:] = 0.0
att_model.eta[:,:,:] = 0.0

and write the model to a HDF5 file models/constant_velocity_N61_51_101_PyTomoATT.h5


# write the model to a file
fname = "%s/constant_velocity_N%d_%d_%d_PyTomoATT.h5"%(output_path, n_rtp[0], n_rtp[1], n_rtp[2])
att_model.write(fname)
print("generate model using PyTomoATT:", fname)

Using loop

The model file in HDF5 format has a simple structure. We can simply use h5py module to build the model file using loop.

First, we initialize the 3D numpy array of velocity vel, and anisotropy xi and eta.


vel = np.zeros(n_rtp)
xi  = np.zeros(n_rtp)
eta = np.zeros(n_rtp)

Second, we assign a constant velocity to the velocity array


for ir in range(n_rtp[0]):
    for it in range(n_rtp[1]):
        for ip in range(n_rtp[2]):
            vel[ir, it, ip] = constant_v

Finally, we can use h5py module to create the HDF5 model file as models/constant_velocity_N61_51_101_loop.h5


fname = "%s/constant_velocity_N%d_%d_%d_loop.h5"%(output_path, n_rtp[0], n_rtp[1], n_rtp[2])
 
with h5py.File(fname, 'w') as f:
    f.create_dataset('vel', data=vel)
    f.create_dataset('xi',  data=xi)
    f.create_dataset('eta', data=eta)
print("generate model using plain loop:", fname)

2. Build linear velocity model

Here, we also provide two ways to generate a model with linear velocity. Specifically, the model satisfies:

a constant velocity of 5.0 km/s above 0 km depth.
a constant velocity gradient from 5.0 km/s at 0 km to 8.0 km/s at 40 km.
a constant velocity of 8.0 km/s below 40 km depth.

Using PyTomoATT

The model file is built as models/linear_velocity_N61_51_101_PyTomoATT.h5


# case 1. ---------- generate a constant velocity model using PyTomoATT module -------------
 
# set the velocity model to a constant value
idx = np.where((am_depths >= 0.0) & (am_depths <= 30.0))
depth = am_depths[idx]
att_model.vel[idx,:,:] = 5.0 + 0.1 * depth[:, np.newaxis, np.newaxis]   # velocity increases linearly from 5.0 to 8.0 km/s
att_model.vel[np.where(am_depths > 30.0),:,:] = 8.0                     # velocity is constant at 8.0 km/s below 30.0 km depth
att_model.vel[np.where(am_depths <  0.0),:,:] = 5.0                     # velocity is constant at 5.0 km/s above 0.0 km depth
 
att_model.xi[:,:,:] = 0.0
att_model.eta[:,:,:] = 0.0
 
# write the model to a file
fname = "%s/linear_velocity_N%d_%d_%d_PyTomoATT.h5"%(output_path, n_rtp[0], n_rtp[1], n_rtp[2])
att_model.write(fname)
print("generate model using PyTomoATT:", fname)

Using loop

The model file is built as models/linear_velocity_N61_51_101_loop.h5


# case 2. ---------- generate a linear velocity model using plain loop (h5py module is required) -------------

# set the velocity model to a linear value
vel = np.zeros(n_rtp)
xi  = np.zeros(n_rtp)
eta = np.zeros(n_rtp)

for ir in range(n_rtp[0]):
    for it in range(n_rtp[1]):
        for ip in range(n_rtp[2]):
            if am_depths[ir] < 0.0:
                vel[ir, it, ip] = 5.0
            elif am_depths[ir] <= 30.0:
                vel[ir, it, ip] = 5.0 + 0.1 * am_depths[ir]
            else:
                vel[ir, it, ip] = 8.0
fname = "%s/linear_velocity_N%d_%d_%d_loop.h5"%(output_path, n_rtp[0], n_rtp[1], n_rtp[2])

with h5py.File(fname, 'w') as f:
    f.create_dataset('vel', data=vel)
    f.create_dataset('xi',  data=xi)
    f.create_dataset('eta', data=eta)
print("generate model using plain loop:", fname)

3. Build checkerboard model

Here, we build the checkerboard model by assigning positive and negative velocity perturbations to the isotropic linear velocity model in a staggered way. In addition, azimuthal anisotropy with orthogonal fast velocity directions are also added to the isotropic linear velocity model.


lim_x       = [0.5, 1.5]    # longitude limits of the checkerboard
lim_y       = [0.25, 0.75]  # latitude limits of the checkerboard
lim_z       = [0, 30]       # depth limits of the checkerboard
pert_vel    = 0.1           # amplitude of velocity perturbation (%)
pert_ani    = 0.05          # amplitude of anisotropy perturbation (fraction)
ani_dir     = 60.0          # fast velocity direction (anti-clockwise from x-axis, in degrees)
n_pert_x    = 4             # number of checkers in x (lon) direction
n_pert_y    = 2             # number of checkers in y (lat) direction
n_pert_z    = 3             # number of checkers in z (dep) direction
 
size_x      = (lim_x[1] - lim_x[0]) / n_pert_x  # size of each checker in x direction
size_y      = (lim_y[1] - lim_y[0]) / n_pert_y  # size of each checker in y direction
size_z      = (lim_z[1] - lim_z[0]) / n_pert_z  # size of each checker in z direction

Related parameter settings denoted here

The range of assigned perturbations is 0.5-1.5 degree in longitude (lim_x), 0.25-0.75 degree in latitude (lim_y), and 0-30 km in depth (lim_z);
The maximum perturbation of the velocity is 10% (pert_vel = 0.1);
The maximum amplitude of anisotropy is 0.05 (pert_ani);
The directions of fast velocity are 60 degree and 150 degree (anti-clockwise from x-axis) (ani_dir = 60.0)
The number of checkers are 4 in longitude (n_pert_x), 2 in latitude (n_pert_y), and 3 in depth (n_pert_z).
Correspondingly, the sizes of checkers are 0.25 degree in longitude (size_x) , 0.25 degree in latitude (size_y) , and 10 km in depth (size_z)

Using PyTomoATT

The model file is built as models/linear_velocity_ckb_N61_51_101_PyTomoATT.h5


# case 1. ---------- generate a constant velocity model using PyTomoATT module -------------
 
# file name of the background model
bg_model_fname = "%s/linear_velocity_N%d_%d_%d_PyTomoATT.h5" % (output_path, n_rtp[0], n_rtp[1], n_rtp[2])
 
ckb = Checker(bg_model_fname, para_fname=par_file)
# n_pert_x, n_pert_y, n_pert_z: number of checkers in x (lon), y (lat), z (dep) directions
# pert_vel: amplitude of velocity perturbation (km/s)
# pert_ani: amplitude of anisotropy perturbation (fraction)
# ani_dir: fast velicty direction (anti-cloclkwise from x-axis, in degrees)
# lim_x, lim_y, lim_z: limits of the checkerboard in x (lon), y (lat), z (dep) directions
ckb.checkerboard(
    n_pert_x=n_pert_x, n_pert_y=n_pert_y, n_pert_z=n_pert_z,         
    pert_vel=pert_vel, pert_ani=pert_ani, ani_dir=ani_dir,   
    lim_x=lim_x, lim_y=lim_y, lim_z=lim_z
)
 
fname = "%s/linear_velocity_ckb_N%d_%d_%d_PyTomoATT.h5" % (output_path, n_rtp[0], n_rtp[1], n_rtp[2])
ckb.write(fname)
print("generate checkerboard model based on the linear velocity model using PyTomoATT:", fname)

Using loop

The model file is built as models/linear_velocity_ckb_N61_51_101_loop.h5. The following script is a little complex, which does the same thing as PyTomoATT above.


# case 2. ---------- generate a checkerboard model using plain loop (h5py module is required) -------------
 
# read the background model
bg_model = np.zeros(n_rtp)
with h5py.File(bg_model_fname, 'r') as f:
    bg_model = f['vel'][:]
 
# set the checkerboard model
vel = np.zeros(n_rtp)
xi  = np.zeros(n_rtp)
eta = np.zeros(n_rtp)
 
for ir in range(n_rtp[0]):
    for it in range(n_rtp[1]):
        for ip in range(n_rtp[2]):
            depth = am_depths[ir]
            lat = am_latitudes[it]
            lon = am_longitudes[ip]
 
            # check if the current grid node is within the checkerboard limits
            if (lim_x[0] <= lon <= lim_x[1]) and (lim_y[0] <= lat <= lim_y[1]) and (lim_z[0] <= depth <= lim_z[1]):
                
                sigma_vel = np.sin(np.pi * (lon - lim_x[0])/size_x) * np.sin(np.pi * (lat - lim_y[0])/size_y) * np.sin(np.pi * (depth - lim_z[0])/size_z)
                sigma_ani = np.sin(np.pi * (lon - lim_x[0])/size_x) * np.sin(np.pi * (lat - lim_y[0])/size_y) * np.sin(np.pi * (depth - lim_z[0])/size_z)
 
                if (sigma_ani > 0):
                    psi = ani_dir / 180.0 * np.pi  # convert degrees to radians
                elif (sigma_ani < 0):
                    psi = (ani_dir + 90.0) / 180.0 * np.pi
                else:
                    psi = 0.0
 
            else:
                sigma_vel = 0.0
                sigma_ani = 0.0
                psi = 0.0
 
            # set the velocity and anisotropy
            vel[ir, it, ip] = bg_model[ir, it, ip] * (1.0 + pert_vel * sigma_vel)
            xi[ir, it, ip]  = pert_ani * abs(sigma_ani) * np.cos(2*psi)
            eta[ir, it, ip] = pert_ani * abs(sigma_ani) * np.sin(2*psi)
 
# write the model to a file
fname = "%s/linear_velocity_ckb_N%d_%d_%d_loop.h5" % (output_path, n_rtp[0], n_rtp[1], n_rtp[2])
with h5py.File(fname, 'w') as f:
    f.create_dataset('vel', data=vel)
    f.create_dataset('xi',  data=xi)
    f.create_dataset('eta', data=eta)
print("generate checkerboard model based on the linear velocity model using plain loop:", fname)

4. Build flexible checkerboard model

Here, we build a checkboard model with varying checker sizes:

0 - 8 km, checker sizes are 0.2 by 0.2 for velocity and 0.3 by 0.3 for anisotropy;
8 - 20 km, checker sizes are 0.3 by 0.3 for velocity and 0.4 by 0.4 for anisotropy;
20 - 36 km, checker sizes are 0.4 by 0.4 for velocity and 0.5 by 0.5 for anisotropy;

The model file is built as models/linear_velocity_ckb_flex_N61_51_101.h5.


# file name of the background model
bg_model_fname = "%s/linear_velocity_N%d_%d_%d_PyTomoATT.h5" % (output_path, n_rtp[0], n_rtp[1], n_rtp[2])
 
# read the background model
bg_model = np.zeros(n_rtp)
with h5py.File(bg_model_fname, 'r') as f:
    bg_model = f['vel'][:]
 
# set the checkerboard model
vel = np.zeros(n_rtp)
xi  = np.zeros(n_rtp)
eta = np.zeros(n_rtp)
 
for ir in range(n_rtp[0]):
    for it in range(n_rtp[1]):
        for ip in range(n_rtp[2]):
            depth = am_depths[ir]
            lat = am_latitudes[it]
            lon = am_longitudes[ip]
 
            if ((depth >= 0.0) and (depth <= 8.0)):
                size_vel = 0.2
                size_ani = 0.3
 
                sigma_vel = np.sin(np.pi * lon/size_vel) * np.sin(np.pi * lat/size_vel) * np.sin(np.pi * depth/8.0)
                sigma_ani = np.sin(np.pi * lon/size_ani) * np.sin(np.pi * lat/size_ani) * np.sin(np.pi * depth/8.0)
 
            elif ((depth > 8.0) and (depth <= 20.0)):
 
                size_vel = 0.3
                size_ani = 0.4
 
                sigma_vel = np.sin(np.pi * lon/size_vel) * np.sin(np.pi * lat/size_vel) * np.sin(np.pi * (depth - 8.0)/12.0 + np.pi)
                sigma_ani = np.sin(np.pi * lon/size_ani) * np.sin(np.pi * lat/size_ani) * np.sin(np.pi * (depth - 8.0)/12.0 + np.pi)
 
            elif ((depth > 20.0) and (depth <= 36.0)):
 
                size_vel = 0.4
                size_ani = 0.5
 
                sigma_vel = np.sin(np.pi * lon/size_vel) * np.sin(np.pi * lat/size_vel) * np.sin(np.pi * (depth - 20.0)/16.0 + 2*np.pi)
                sigma_ani = np.sin(np.pi * lon/size_ani) * np.sin(np.pi * lat/size_ani) * np.sin(np.pi * (depth - 20.0)/16.0 + 2*np.pi)
 
            else:
                sigma_vel = 0.0
                sigma_ani = 0.0
 
            if (sigma_ani > 0):
                psi = ani_dir / 180.0 * np.pi  # convert degrees to radians
            elif (sigma_ani < 0):
                psi = (ani_dir + 90.0) / 180.0 * np.pi
            else:
                psi = 0.0
 
            # set the velocity and anisotropy
            vel[ir, it, ip] = bg_model[ir, it, ip] * (1.0 + pert_vel * sigma_vel)
            xi[ir, it, ip]  = pert_ani * abs(sigma_ani) * np.cos(2*psi)
            eta[ir, it, ip] = pert_ani * abs(sigma_ani) * np.sin(2*psi)
 
# write the model to a file
fname = "%s/linear_velocity_ckb_flex_N%d_%d_%d.h5" % (output_path, n_rtp[0], n_rtp[1], n_rtp[2])
with h5py.File(fname, 'w') as f:
    f.create_dataset('vel', data=vel)
    f.create_dataset('xi',  data=xi)
    f.create_dataset('eta', data=eta)
 
print("generate flexible checkerboard model based on the linear velocity model using plain loop:", fname)