# -*- coding: utf-8 -*-
"""
Module to produce traveltime lookup tables defined on a Cartesian grid.
:copyright:
2020, QuakeMigrate developers.
:license:
GNU General Public License, Version 3
(https://www.gnu.org/licenses/gpl-3.0.html)
"""
import logging
import os
import pathlib
import struct
from shutil import rmtree
import numpy as np
import pyproj
from scipy.interpolate import interp1d
import quakemigrate.util as util
from .lut import LUT
[docs]def read_nlloc(path, stations, phases=["P", "S"], fraction_tt=0.1, log=False):
"""
Read in a traveltime lookup table that is saved in the NonLinLoc format.
Parameters
----------
path : str
Path to directory containing .buf and .hdr files.
stations : `pandas.DataFrame`
DataFrame containing station information (lat/lon/elev).
phases : list of str, optional
List of seismic phases for which to calculate traveltimes.
fraction_tt : float, optional
An estimate of the uncertainty in the velocity model as a function of
a fraction of the traveltime. (Default 0.1 == 10%)
log : bool, optional
Toggle for logging - default is to only print information to stdout.
If True, will also create a log file.
Returns
-------
lut : :class:`~quakemigrate.lut.LUT` object
Lookup table populated with traveltimes from the NonLinLoc files.
"""
path = pathlib.Path(path)
util.logger(pathlib.Path.cwd() / "logs" / "lut", log)
logging.info("Loading NonLinLoc traveltime lookup tables for...")
for i, phase in enumerate(phases):
logging.info(f"\t...phase: {phase}...")
for j, station in enumerate(stations["Name"].values):
logging.info(f"\t\t...station: {station}")
file = path / f"layer.{phase}.{station}.time"
if i == 0 and j == 0:
gridspec, transform, traveltimes = _read_nlloc(file)
node_count = np.array(gridspec[0])
grid_origin = np.array(gridspec[1])
node_spacing = np.array(gridspec[2])
gproj, cproj, gproj_string = transform
if gproj is None:
raise NotImplementedError(f"Projection type {gproj_string}"
" not supported.")
# Transform from grid projection origin to a coord origin
ll_corner = pyproj.transform(gproj, cproj, *grid_origin)
# Calculate the ur corner
ur_corner = (np.array(grid_origin)
+ (node_count - 1) * node_spacing)
ur_corner = pyproj.transform(gproj, cproj, *ur_corner)
# Need to initialise the grid
lut = LUT(ll_corner=ll_corner, ur_corner=ur_corner,
node_spacing=node_spacing, grid_proj=gproj,
coord_proj=cproj, fraction_tt=fraction_tt)
else:
_, _, traveltimes = _read_nlloc(file)
lut.traveltimes.setdefault(station, {}).update(
{phase: traveltimes})
lut.station_data = stations
lut.phases = phases
return lut
[docs]def compute_traveltimes(grid_spec, stations, method, phases=["P", "S"],
fraction_tt=0.1, save_file=None, log=False, **kwargs):
"""
Top-level method for computing traveltime lookup tables.
This function takes a grid specification and is capable of computing
traveltimes for an arbitrary number of phases using a variety of
techniques.
Parameters
----------
grid_spec : dict
Dictionary containing all of the defining parameters for the underlying
3-D grid on which the traveltimes are to be calculated. For expected
keys, see :class:`~quakemigrate.lut.lut.Grid3D`.
stations : `pandas.DataFrame`
DataFrame containing station information (lat/lon/elev).
method : str
Method to be used when computing the traveltime lookup tables.
"homogeneous" - straight line velocities
"1dfmm" - 1-D fast-marching method using scikit-fmm
"1dsweep" - a 2-D traveltime grid for a 1-D velocity model is swept
over the 3-D grid using a bilinear interpolation scheme
phases : list of str, optional
List of seismic phases for which to calculate traveltimes.
fraction_tt : float, optional
An estimate of the uncertainty in the velocity model as a function of
a fraction of the traveltime. (Default 0.1 == 10%)
filename : str, optional
Path to location to save pickled lookup table.
log : bool, optional
Toggle for logging - default is to only print information to stdout.
If True, will also create a log file.
kwargs : dict
Dictionary of all keyword arguments passed to compute when called.
For lists of valid arguments, please refer to the relevant method.
Returns
-------
lut : :class:`~quakemigrate.lut.LUT` object
Lookup table populated with traveltimes from the NonLinLoc files.
"""
util.logger(pathlib.Path.cwd() / "logs" / "lut", log)
lut = LUT(**grid_spec, fraction_tt=fraction_tt)
lut.station_data = stations
lut.phases = phases
if method == "homogeneous":
logging.info("Computing homogeneous traveltimes for...")
lut.velocity_model = "Homogeneous velocity model:"
for phase in phases:
velocity = kwargs.get(f"v{phase.lower()}")
if velocity is None:
raise TypeError(f"Missing argument: 'v{phase.lower()}'")
lut.velocity_model += f"\n\tV{phase.lower()} = {velocity:5.2f} m/s"
logging.info(f"\t...phase: {phase}...")
_compute_homogeneous(lut, phase, velocity)
if method == "1dfmm":
logging.info("Computing 1-D fast-marching traveltimes for...")
lut.velocity_model = vmodel = kwargs.get("vmod")
if vmodel is None:
raise TypeError("Missing argument: 'vmod'")
for phase in phases:
logging.info(f"\t...phase: {phase}...")
_compute_1d_fmm(lut, phase, vmodel)
if method == "3dfmm":
raise NotImplementedError("Feature coming soon - please contact the "
"QuakeMigrate developers")
if method == "1dsweep":
logging.info("Computing 1-D sweep traveltimes for...")
lut.velocity_model = vmodel = kwargs.get("vmod")
if vmodel is None:
raise TypeError("Missing argument: 'vmod'")
for phase in phases:
logging.info(f"\t...phase: {phase}...")
_compute_1d_sweep(lut, phase, vmodel, **kwargs)
if save_file is not None:
lut.save(save_file)
return lut
def _compute_homogeneous(lut, phase, velocity):
"""
Calculate the traveltime lookup table for a station in a homogeneous
velocity model.
Parameters
----------
lut : :class:`~quakemigrate.lut.LUT` object
Defines the grid on which the traveltimes are to be calculated.
phase : str
The seismic phase for which to calculate traveltimes.
velocity : float
Seismic phase velocity.
"""
grid_xyz = lut.grid_xyz
stations_xyz = lut.stations_xyz
for i, station in enumerate(lut.station_data["Name"].values):
logging.info(f"\t\t...station: {station} - {i+1} of "
f"{stations_xyz.shape[0]}")
dx, dy, dz = [grid_xyz[j] - stations_xyz[i, j] for j in range(3)]
dist = np.sqrt(dx**2 + dy**2 + dz**2)
lut.traveltimes.setdefault(station, {}).update(
{phase: dist / velocity})
def _compute_1d_fmm(lut, phase, vmodel):
"""
Calculate traveltime lookup tables for each station in a 1-D velocity
model using the fast-marching method.
Parameters
----------
lut : :class:`~quakemigrate.lut.LUT` object
Defines the grid on which the traveltimes are to be calculated.
phase : str
The seismic phase for which to calculate traveltimes.
vmodel : `pandas.DataFrame` object
DataFrame containing the velocity model to be used to generate the LUT.
Columns:
"Depth" of each layer in model (positive down)
"V<phase>" velocity for each layer in model (e.g. "Vp")
"""
try:
depths, vmodel = vmodel[["Depth", f"V{phase.lower()}"]].values.T
except KeyError:
raise util.InvalidVelocityModelHeader(f"V{phase.lower()}")
finfo = np.finfo(float)
depths = np.insert(np.append(depths, finfo.max), 0, finfo.min)
vmodel = np.insert(np.append(vmodel, vmodel[-1]), 0, vmodel[0])
grid_xyz = lut.grid_xyz
stations_xyz = lut.stations_xyz
# Interpolate the velocity model in the Z-dimension
f = interp1d(depths, vmodel)
int_vmodel = f(grid_xyz[2])
for i, station in enumerate(lut.station_data["Name"].values):
logging.info(f"\t\t...station: {station} - {i+1} of "
f"{stations_xyz.shape[0]}")
lut.traveltimes.setdefault(station, {}).update(
{phase: _eikonal_fmm(grid_xyz, lut.node_spacing,
int_vmodel, stations_xyz[i])})
def _eikonal_fmm(grid_xyz, node_spacing, velocity_grid, station_xyz):
"""
Calculates the traveltime lookup tables by solving the eikonal equation
using an implementation of the fast-marching algorithm.
Traveltime calculation can only be performed between grid nodes: the
station location is therefore taken as the closest grid node. Note that
for large node spacings this may cause a modest error in the calculated
traveltimes.
.. warning:: Requires the scikit-fmm python package.
Parameters
----------
grid_xyz : array-like
[X, Y, Z] coordinates of each node.
node_spacing : array-like
[X, Y, Z] distances between each node.
velocity_grid : array-like
Contains the speed of interface propagation at each point in the
domain.
station_xyz : array-like
Station location (in grid xyz).
Returns
-------
traveltimes : array-like, same shape as grid_xyz
Contains the traveltime from the zero contour (zero level set) of phi
to each point in the array given the scalar velocity field speed. If
the input array speed has values less than or equal to zero the return
value will be a masked array.
"""
try:
import skfmm
except ImportError:
raise ImportError("Unable to import skfmm - you need to install "
"scikit-fmm to use this method.\nSee the "
"installation instructions in the documentation"
"for more details.")
phi = -np.ones(grid_xyz[0].shape)
# Find closest grid node to true station location
indx = np.argmin(abs(grid_xyz[0] - station_xyz[0])
+ abs(grid_xyz[1] - station_xyz[1])
+ abs(grid_xyz[2] - station_xyz[2]))
phi[np.unravel_index(indx, grid_xyz[0].shape)] = 1.0
return skfmm.travel_time(phi, velocity_grid, dx=node_spacing)
def _compute_1d_sweep(lut, phase, vmodel, **kwargs):
"""
Calculate 3-D traveltime lookup tables from a 1-D velocity model.
NonLinLoc Grid2Time is used to generate a 2-D lookup table which is then
swept around the full range of azimuths, centred on the station location,
to populate the 3-D traveltime grid.
.. warning:: Requires NonLinLoc to be installed, and `Vel2Grid` and
`Grid2Time` to be in the user's path. Alternatively, a custom path to
these executables can be specified with the kwarg `nlloc_path`.
Parameters
----------
lut : :class:`~quakemigrate.lut.LUT` object
Defines the grid on which the traveltimes are to be calculated.
phase : str
The seismic phase for which to calculate traveltimes.
vmodel : `pandas.DataFrame` object
DataFrame containing the velocity model to be used to generate the LUT.
Columns:
"Depth" of each layer in model (positive down)
"V<phase>" velocity for each layer in model (e.g. "Vp")
kwargs : dict
Can contain:
nlloc_dx : float, optional
NLLoc 2D grid spacing (default: 0.1 km). Note: units must be km.
nlloc_path : str, optional
Path to NonLinLoc executables Vel2Grid and Grid2Time (default: "").
block_model : bool, optional
Toggle to choose whether to interpret velocity model with constant
velocity blocks or a linear gradient (default: False).
retain_nll_grids : bool, optional
Toggle to choose whether to keep the 2-D traveltime grids created
by NonLinLoc (default: False).
"""
from subprocess import check_output, STDOUT
# Unpack kwargs
nlloc_dx = kwargs.get("nlloc_dx", 0.1)
nlloc_path = pathlib.Path(kwargs.get("nlloc_path", ""))
block_model = kwargs.get("block_model", False)
retain_nll_grids = kwargs.get("retain_nll_grids", False)
# Check nlloc_path is valid and contains Vel2Grid and Grid2Time
if kwargs.get("nlloc_path", "") != "":
if not (nlloc_path / "Vel2Grid").exists() \
or not (nlloc_path / "Grid2Time").exists():
raise FileNotFoundError(f"Incorrect nlloc_path? - Grid2Time and "
f"Vel2Grid not found in {nlloc_path}")
# For NonLinLoc, distances/velocities must be in km - use conversion factor
km_cf = 1000 / lut.unit_conversion_factor
grid_xyz = [g / km_cf for g in lut.grid_xyz]
stations_xyz = lut.stations_xyz / km_cf
ll, *_, ur = lut.grid_corners / km_cf
vmodel = vmodel / km_cf
# Make folders in which to run NonLinLoc
cwd = pathlib.Path.cwd()
(cwd / "time").mkdir(exist_ok=True)
(cwd / "model").mkdir(exist_ok=True)
for i, station in enumerate(lut.station_data["Name"].values):
logging.info(f"\t\t...running NonLinLoc - station: {station:5s} - "
f"{i+1} of {stations_xyz.shape[0]}")
dx, dy = [grid_xyz[j] - stations_xyz[i, j] for j in range(2)]
distances = np.sqrt(dx**2 + dy**2).flatten()
depths = grid_xyz[2].flatten()
max_dist = np.max(distances)
# NLLoc needs the station to lie within the 2-D section -> we pick the
# depth extent of the 2-D grid from the max. possible extent of the
# station and the grid - [min_z, max_z]
depth_span = [np.min([ll[2], stations_xyz[i, 2]]),
np.max([ur[2], stations_xyz[i, 2]])]
# Allow 2 nodes on depth extent as a computational buffer
_write_control_file(stations_xyz[i], station, max_dist, vmodel,
depth_span, phase, nlloc_dx, block_model)
for mode in ["Vel2Grid", "Grid2Time"]:
out = check_output([str(nlloc_path / mode),
"control.in"], stderr=STDOUT)
if b"ERROR" in out:
raise Exception(f"{mode} Error", out)
to_read = cwd / "time" / f"layer.{phase}.{station}.time"
gridspec, _, traveltimes = _read_nlloc(to_read, ignore_proj=True)
lut.traveltimes.setdefault(station, {}).update(
{phase: _bilinear_interpolate(np.c_[distances, depths],
gridspec[1, 1:],
gridspec[2, 1:],
traveltimes[0, :, :]
).reshape(lut.node_count)})
# Tidy up: remove control file and nll model and time files
os.remove(cwd / "control.in")
if not retain_nll_grids:
for file in (cwd / "time").glob(f"layer.{phase}.{station}.time*"):
file.unlink()
for file in (cwd / "time").glob(f"layer.{phase}.mod.*"):
file.unlink()
if not retain_nll_grids:
rmtree(cwd / "model")
# Check time directory is empty before removing (might have
# saved grids from previous runs)
if not os.listdir(cwd / "time"):
rmtree(cwd / "time")
else:
logging.info("Warning: time directory not empty; does it contain "
"traveltime grids from a previous run?\nNot removed.")
def _write_control_file(station_xyz, station, max_dist, vmodel, depth_span,
phase, dx, block_model):
"""
Write out a control file for NonLinLoc.
Parameters
----------
station_xyz : array-like
Station location expressed in the coordinate space of the grid, in km.
station : str
Station name.
max_dist : float
Maximum distance between the station and any point in the grid, in km.
vmodel : `pandas.DataFrame` object
DataFrame containing the velocity model to be used to generate the LUT.
Columns:
"Depth" of each layer in model (positive down), in km.
"V<phase>" velocity for each layer in model (e.g. "Vp"), in km / s.
depth_span : array-like
Minimum/maximum extent of the grid in the z-dimension, in km.
phase : str
The seismic phase for which to calculate traveltimes.
dx : float
NLLoc 2D grid spacing, in km.
block_model : bool
Toggle to choose whether to interpret velocity model with constant
velocity blocks or a linear gradient.
"""
control_string = ("CONTROL 0 54321\n"
"TRANS NONE\n\n"
"VGOUT {model_path:s}\n"
"VGTYPE {phase:s}\n\n"
"VGGRID {grid:s} SLOW_LEN\n\n"
"{vmodel:s}\n\n"
"GTFILES {model_path:s} {time_path:s} {phase:s}\n"
"GTMODE GRID2D ANGLES_NO\n\n"
"GTSRCE {station:s} XYZ {x:f} {y:f} {z:f} 0.0\n\n"
"GT_PLFD 1.0E-3 0")
cwd = pathlib.Path.cwd()
out = control_string.format(phase=phase,
grid=_grid_string(max_dist, depth_span, dx),
vmodel=_vmodel_string(vmodel, block_model,
phase),
model_path=str(cwd / "model" / "layer"),
time_path=str(cwd / "time" / "layer"),
station=station, x=station_xyz[0],
y=station_xyz[1], z=station_xyz[2])
with open(cwd / "control.in", "w") as f:
f.write(out)
def _read_nlloc(fname, ignore_proj=False):
"""
Read traveltime lookup tables saved in the NonLinLoc format.
Parses in the header of a NonLinLoc file for the grid specification and the
coordinate transforms. Then unpacks the binary buffer file containing the
traveltimes.
Parameters
----------
fname : str
Path to file containing NonLinLoc traveltime lookup tables, without
the extension.
ignore_proj : bool (optional)
Flag to suppress the "No projection specified message".
Returns
-------
traveltimes : array-like
Traveltimes for the station.
transform : array of `pyproj.Proj` objects
Array containing the grid and coordinate projections, respectively.
gridspec : array-like
Details on the NonLinLoc grid specification. Contains the number of
nodes, the grid origin and the node spacings.
"""
with open(f"{fname}.hdr", "r") as f:
# Read the grid definition line
line = f.readline().split()
nx, ny, nz = int(line[0]), int(line[1]), int(line[2])
x0, y0, z0 = float(line[3]), float(line[4]), float(line[5])
dx, dy, dz = float(line[6]), float(line[7]), float(line[8])
# Read the station information line
_ = f.readline().split()
# Read the transform definition line
line = f.readline().split()
cproj = pyproj.Proj(proj="longlat", ellps="WGS84", datum="WGS84",
no_defs=True)
gproj = None
if line[1] == "NONE" and not ignore_proj:
logging.info("\tNo projection selected.")
elif line[1] == "SIMPLE":
orig_lat = float(line[3])
orig_lon = float(line[5])
# The simple projection is the Plate Carree/Equidistant Cylindrical
gproj = pyproj.Proj(proj="eqc", lat_0=orig_lat, lon_0=orig_lon,
units="km")
elif line[1] == "LAMBERT":
proj_ellipsoid = line[3]
orig_lat = float(line[5])
orig_lon = float(line[7])
parallel_1 = float(line[9])
parallel_2 = float(line[11])
if proj_ellipsoid == "WGS-84":
proj_ellipsoid = "WGS84"
elif proj_ellipsoid == "GRS-80":
proj_ellipsoid = "GRS80"
elif proj_ellipsoid == "WGS-72":
proj_ellipsoid = "WGS72"
elif proj_ellipsoid == "Australian":
proj_ellipsoid = "aust_SA"
elif proj_ellipsoid == "Krasovsky":
proj_ellipsoid = "krass"
elif proj_ellipsoid == "International":
proj_ellipsoid = "intl"
elif proj_ellipsoid == "Hayford-1909":
proj_ellipsoid = "intl"
elif proj_ellipsoid == "Clarke-1880":
proj_ellipsoid = "clrk80"
elif proj_ellipsoid == "Clarke-1866":
proj_ellipsoid = "clrk66"
elif proj_ellipsoid == "Airy":
proj_ellipsoid = "airy"
elif proj_ellipsoid == "Bessel":
proj_ellipsoid = "bessel"
elif proj_ellipsoid == "Hayford-1830":
proj_ellipsoid = "evrst30"
elif proj_ellipsoid == "Sphere":
proj_ellipsoid = "sphere"
else:
logging.info(f"Projection Ellipsoid {proj_ellipsoid} not "
"supported!\n\tPlease notify the QuakeMigrate "
"developers.\n\n\tWGS-84 used instead...\n")
proj_ellipsoid = "WGS84"
gproj = pyproj.Proj(proj="lcc", lon_0=orig_lon, lat_0=orig_lat,
lat_1=parallel_1, lat_2=parallel_2, units="km",
ellps=proj_ellipsoid, datum="WGS84")
elif line[1] == "TRANS_MERC":
orig_lat = float(line[5])
orig_lon = float(line[7])
gproj = pyproj.Proj(proj="tmerc", lon=orig_lon, lat=orig_lat,
units="km")
transform = [gproj, cproj, line[1]]
with open(f"{fname}.buf", "rb") as f:
npts = nx * ny * nz
buf = f.read(npts * 4)
traveltimes = struct.unpack("f" * npts, buf)
traveltimes = np.array(traveltimes).reshape((nx, ny, nz))
gridspec = np.array([[nx, ny, nz], [x0, y0, z0], [dx, dy, dz]])
return gridspec, transform, traveltimes
def _bilinear_interpolate(xz, xz_origin, xz_dimensions, traveltimes):
"""
Perform a bi-linear interpolation between 4 data points on the input
2-D lookup table to calculate the traveltime to nodes on the 3-D grid.
Note: x is used to denote the horizontal dimension over which the
interpolation is performed. Due to the NonLinLoc definition, this
corresponds to the y component of the grid.
Parameters
----------
xz : array-like
Column-stacked array of distances from the station and depths for all
points in grid.
xz_origin : array-like
The x (actually y) and z values of the grid origin.
xz_dimensions : array-like
The x (actually y) and z values of the node spacing.
traveltimes : array-like
A slice through the traveltime grid at x = 0, on which to perform the
interpolation.
Returns
-------
int_traveltimes : array-like
Interpolated 3-D traveltime lookup table.
"""
# Get indices of the nearest node
i, k = np.floor((xz - xz_origin) / xz_dimensions).astype(int).T
# Get fractional distance along each axis
x_d, z_d = (np.remainder(xz, xz_dimensions) / xz_dimensions).T
# Get 4 data points of surrounding square
c00 = traveltimes[i, k]
c10 = traveltimes[i+1, k]
c11 = traveltimes[i+1, k+1]
c01 = traveltimes[i, k+1]
# Interpolate along x-axis
c0 = c00 * (1 - x_d) + c10 * x_d
c1 = c01 * (1 - x_d) + c11 * x_d
# Interpolate along z-axis (c = c0 if np.all(z_d == 0.))
return c0 * (1 - z_d) + c1 * z_d
def _vmodel_string(vmodel, block_model, phase):
"""
Creates a string representation of the velocity model for the control file.
Parameters
----------
vmodel : `pandas.DataFrame` object
DataFrame containing the velocity model to be used to generate the LUT.
Columns:
"Depth" of each layer in model (positive down), in km.
"V<phase>" velocity for each layer in model (e.g. "Vp"), in km / s.
block_model : bool
Toggle to choose whether to interpret velocity model with constant
velocity blocks or a linear gradient.
phase : str
The seismic phase for which to calculate traveltimes.
Returns
-------
str_vmodel : str
NonLinLoc formatted string describing the velocity model.
"""
string_template = "LAYER {0:f} {1:f} {2:f} {1:f} {2:f} 0.0 0.0"
out = []
i = 0
while i < len(vmodel):
if not block_model:
try:
dvdx = _velocity_gradient(i, vmodel, phase)
except KeyError:
dvdx = 0.0
else:
dvdx = 0.0
out.append(string_template.format(vmodel["Depth"][i],
vmodel[f"V{phase.lower()}"][i],
dvdx))
i += 1
return "\n".join(out)
def _velocity_gradient(i, vmodel, phase):
"""
Calculate the gradient of the velocity model between two layers.
Parameters
----------
i : int
Index of upper layer.
vmodel : `pandas.DataFrame` object
DataFrame containing the velocity model to be used to generate the LUT.
Columns:
"Depth" of each layer in model (positive down), in km.
"V<phase>" velocity for each layer in model (e.g. "Vp"), in km / s.
phase : str
The seismic phase for which to calculate traveltimes.
Returns
-------
dvdx : float
Velocity gradients for the block model.
"""
d_depth = vmodel["Depth"][i+1] - vmodel["Depth"][i]
d_vel = vmodel[f"V{phase.lower()}"][i+1] - vmodel[f"V{phase.lower()}"][i]
return d_vel / d_depth
def _grid_string(max_dist, depth_limits, dx):
"""
Creates a string representation of the grid for the control file.
Parameters
----------
max_dist : float
Maximum distance between the station and any point in the grid, in km.
depth_limits : array-like
Minimum/maximum extent of the grid in the z-dimension, in km.
dx : float
NLLoc 2D grid spacing (default: 0.1 km).
Returns
-------
str_grid : str
NonLinLoc formatted string describing the grid.
"""
string_template = "2 {0:d} {1:d} 0.0 0.0 {2:f} {3:f} {3:f} {3:f}"
max_x = int(np.ceil(max_dist / dx)) + 5
max_z = int(np.ceil((depth_limits[1] - depth_limits[0]) / dx)) + 5
return string_template.format(max_x, max_z, depth_limits[0], dx)