2.6.1. Building the traveltime lookup table

For this synthetic example, we use a hypothetical earthquake region situated at (0.0°, 0.0°), spanning ±0.15° in longitude and latitude, and 31 km in depth, spanning from -1 km to 30 km. A randomly distributed network of seismometers are generated, using a fixed random seed, to ensure the same random distribution is generated each time the example is run.

The station file is created by the following section of code:

station_file = "./inputs/synthetic_stations.txt"

# --- Build station file ---
rng = np.random.default_rng(13)  # Fix seed for reproducible results
stations = pd.DataFrame()
stations["Network"] = ["SC"] * 10
stations["Name"] = [f"STA{i}" for i in range(10)]
stations["Longitude"] = rng.uniform(low=-0.15, high=0.15, size=10)
stations["Latitude"] = rng.uniform(low=-0.15, high=0.15, size=10)
stations["Elevation"] = rng.uniform(low=-0.0, high=1.0, size=10)
stations.to_csv(station_file, index=False)

Note

You can change the seed (or remove it) to explore how the synthetic example performs for different distributions of seismometers in the network.

A simple 1-D velocity model is used, with P- and S-phase velocities for depths spanning a range greater than the targeted search depths. The input velocity model is stored in a .csv file and read in using a utility function provided by QuakeMigrate that performs some initial sanity checks, e.g., the .csv must contain at least a Depth and a Vp column.

vmodel_file = "./inputs/velocity_model.csv"

# --- Read in the velocity model file ---
vmodel = read_vmodel(vmodel_file)

We employ the Lambert Conformal Conic projection as the basis of our transformation between the geographical coordinate space and a Cartesian coordinate space, on which the grid search is performed. These projections are defined using the pyproj package, which provides an interface to the PROJ library.

# --- Define the input and grid projections ---
gproj = Proj(
    proj="lcc",
    units="km",
    lon_0=0.0,
    lat_0=0.0,
    lat_1=-0.10,
    lat_2=0.101,
    datum="WGS84",
    ellps="WGS84",
    no_defs=True,
)
cproj = Proj(proj="longlat", datum="WGS84", ellps="WGS84", no_defs=True)

These projections are combined with a specification of the search grid (per the geographical specifications at the top of this tutorial) using the AttribDict class, which is an obspy utility class that can be used to construct an object that behaves like a dictionary (i.e., can be accessed using dict[key]) as well as providing .-style access to the (key, value) pairs.

# --- Define the grid specifications ---
# AttribDict behaves like a Python dict, but also has '.'-style access.
grid_spec = AttribDict()
grid_spec.ll_corner = [-0.15, -0.15, -1.0]
grid_spec.ur_corner = [0.15, 0.15, 30.0]
grid_spec.node_spacing = [0.5, 0.5, 0.5]
grid_spec.grid_proj = gproj
grid_spec.coord_proj = cproj

Finally, we bring all of these parts together to compute the traveltime grids—that is, for every station, we compute the traveltime from the position of the station to each node in the grid, for each phase (here P and S). In this instance, we make use of the NonLinLoc package to compute these traveltimes. For more information, please refer to the dedicated lookup table documentation.

lut_out = "./outputs/lut/example.LUT"

# --- Homogeneous LUT generation ---
lut = compute_traveltimes(
    grid_spec,
    stations,
    method="1dnlloc",
    vmod=vmodel,
    phases=["P", "S"],
    log=True,
    save_file=lut_out,
)
print()
print(lut)

The final line will print out an overview of the travetime lookup table.

The full script looks like this:

"""
This script generates the traveltime lookup table for the synthetic example described in
the tutorial in the online documentation.

:copyright:
    2020–2024, QuakeMigrate developers.
:license:
    GNU General Public License, Version 3
    (https://www.gnu.org/licenses/gpl-3.0.html)

"""

# Stop numpy using all available threads (these environment variables must be
# set before numpy is imported for the first time).
import os

os.environ.update(
    OMP_NUM_THREADS="1",
    OPENBLAS_NUM_THREADS="1",
    NUMEXPR_NUM_THREADS="1",
    MKL_NUM_THREADS="1",
)

import numpy as np
from obspy.core import AttribDict
import pandas as pd
from pyproj import Proj
from quakemigrate.io import read_vmodel
from quakemigrate.lut import compute_traveltimes


# Build synthetic lookup table
station_file = "./inputs/synthetic_stations.txt"
vmodel_file = "./inputs/velocity_model.csv"
lut_out = "./outputs/lut/example.LUT"

# --- Build station file ---
rng = np.random.default_rng(13)  # Fix seed for reproducible results
stations = pd.DataFrame()
stations["Network"] = ["SC"] * 10
stations["Name"] = [f"STA{i}" for i in range(10)]
stations["Longitude"] = rng.uniform(low=-0.15, high=0.15, size=10)
stations["Latitude"] = rng.uniform(low=-0.15, high=0.15, size=10)
stations["Elevation"] = rng.uniform(low=-0.0, high=1.0, size=10)
stations.to_csv(station_file, index=False)

# --- Read in the velocity model file ---
vmodel = read_vmodel(vmodel_file)

# --- Define the input and grid projections ---
gproj = Proj(
    proj="lcc",
    units="km",
    lon_0=0.0,
    lat_0=0.0,
    lat_1=-0.10,
    lat_2=0.101,
    datum="WGS84",
    ellps="WGS84",
    no_defs=True,
)
cproj = Proj(proj="longlat", datum="WGS84", ellps="WGS84", no_defs=True)

# --- Define the grid specifications ---
# AttribDict behaves like a Python dict, but also has '.'-style access.
grid_spec = AttribDict()
grid_spec.ll_corner = [-0.15, -0.15, -1.0]
grid_spec.ur_corner = [0.15, 0.15, 30.0]
grid_spec.node_spacing = [0.5, 0.5, 0.5]
grid_spec.grid_proj = gproj
grid_spec.coord_proj = cproj

# --- Homogeneous LUT generation ---
lut = compute_traveltimes(
    grid_spec,
    stations,
    method="1dnlloc",
    vmod=vmodel,
    phases=["P", "S"],
    log=True,
    save_file=lut_out,
)
print()
print(lut)