import glob
from os.path import exists, join, dirname
import numpy as np
from matplotlib import pyplot as plt
from scipy.interpolate import interp1d, interpn
from seispy.utils import vs2vprho
import warnings
warnings.filterwarnings("ignore", "invalid value encountered in sqrt")
def _search_vel_file(mode_name):
"""
search vel file given by mode_name
and try to open it with np.loadtxt:
1. precise name of vel mod file
2. default vel mod (containing iasp01
default vel mod stored in seispy/data/modname.vel file
do basic value evaluation, cal rho if not given in files
Parameters
----------
mode_name
Returns
----------
model matrix [layers:4]:
depth vp vs rho
.... ... ... ....
"""
if not isinstance(mode_name, str):
raise TypeError('velmod should be in str type')
if exists(mode_name):
filename = mode_name
## if not found, search in default vel model fold
elif exists(join(dirname(dirname(__file__)), 'data', mode_name.lower() + '.vel')):
filename = join(dirname(dirname(__file__)), 'data', mode_name.lower() + '.vel')
else:
raise ValueError('No such file of velocity model')
try:
raw_model = np.loadtxt(filename)
except:
raise IOError
# check cols of file
if raw_model.shape[1] < 3 or raw_model.shape[1] > 4:
raise ValueError('The file should contain 3 or 4 columns')
# cal rho if rho is not given in vel files.
if raw_model.shape[1] == 3:
model = np.zeros((raw_model.shape[0], 4))
model[:, :3] = raw_model[:, :]
_p, rho = vs2vprho(raw_model[:, 2])
model[:, 3] = rho
return model
else:
return raw_model
def _layer2grid(dep_range, model):
"""
trans model from layer_model to layer_grid
leave stuffing and interp 2 discretize
dep_range : grids at depth axis, np.ndarray
h : thichness of each layer
vp, vs, rho
Returns
-------
"""
neo_model = np.zeros((len(dep_range), 4)).astype(float)
picks = np.searchsorted(model[:, 0], dep_range, side="left")
for _i, _j in enumerate(picks):
neo_model[_i, :] = model[_j, :]
return neo_model[:, 1], neo_model[:, 2], neo_model[:, 3]
def _intep_mod(model, depths_elev):
vp = interp1d(model[:,0], model[:,1], bounds_error=False,
fill_value=model[0,1])(depths_elev)
vs = interp1d(model[:,0], model[:,2], bounds_error=False,
fill_value=model[0,2])(depths_elev)
rho = interp1d(model[:,0], model[:,3], bounds_error=False,
fill_value=model[0,3])(depths_elev)
return vp, vs, rho
def _from_layer_model(dep_range, h, vs, vp=None, rho=None):
dep = 0
vp_dep = np.zeros_like(dep_range).astype(float)
vs_dep = np.zeros_like(dep_range).astype(float)
rho_dep = np.zeros_like(dep_range).astype(float)
for i, layer in enumerate(h):
if (layer == 0 and i == len(h)-1) or \
(dep+layer < dep_range[-1] and i == len(h)-1):
idx = np.where(dep_range>=dep)[0]
else:
idx = np.where((dep_range >= dep) & (dep_range < dep+layer))[0]
if idx.size == 0:
raise ValueError('The thickness of layer {} less than the depth interval'.format(i+1))
vs_dep[idx] = vs[i]
if vp is not None:
vp_dep[idx] = vp[i]
else:
vp_dep[idx], _ = vs2vprho(vs[i])
if rho is not None:
rho_dep[idx] = rho[i]
else:
_, rho_dep[idx] = vs2vprho(vs[i])
dep += layer
if dep > dep_range[-1]:
break
return np.vstack((dep_range, vp_dep, vs_dep, rho_dep)).T
[docs]
class DepModel(object):
"""
Class for constructing 1D velocity model and computing delay time of Ps, PpPs and PsPs Rays.
.. rubric:: Examples
>>> model = DepModel(np.array([0, 20.1, 35.1, 100]))
>>> print(model.dz)
[ 0. 20.1 15. 64.9]
>>> print(model.vp)
[5.8 6.5 8.04001059 8.04764706]
>>> print(model.vs)
[3.36 3.75 4.47003177 4.49294118]
"""
def __init__(self, dep_range, velmod='iasp91', elevation=0., layer_mod=False):
""" Initialize DepModel object
:type dep_range: numpy.ndarray
:param dep_range: Depth range for conversion
:type velmod: str, optional
:param velmod: Text file of 1D velocity model with first 3 columns of depth/thickness, Vp and Vs, by default 'iasp91'
:type elevation: float, optional
:param elevation: Elevation in km, by default 0.0
:type layer_mod: bool, optional
:param layer_mod: True for search, and False for interp1d, by default False
"""
self.isrho = False
self.elevation = elevation
self.layer_mod = layer_mod
# dep layer for CCP or other purpose, relative to sea level, dont contain any depth infomation
self.depths = dep_range.astype(float)
self.dep_val = np.average(np.diff(self.depths))
try:
self.model_array = _search_vel_file(velmod)
except (IOError, TypeError):
return
else:
self._elevation()
if layer_mod:
self.vp, self.vs, self.rho = \
_layer2grid(self.depths, self.model_array)
else:
self.vp, self.vs, self.rho = \
_intep_mod(self.model_array, self.depths_elev)
[docs]
@classmethod
def read_layer_model(cls, dep_range, h, vp, vs, rho=None, elevation=0):
""" Read layer model from given parameters
:type dep_range: numpy.ndarray
:param dep_range: Depth range for conversion
:type h: numpy.ndarray
:param h: Thickness of each layer
:type vp: numpy.ndarray
:param vp: P-wave velocity of each layer
:type vs: numpy.ndarray
:param vs: S-wave velocity of each layer
:type rho: numpy.ndarray
:param rho: Density of each layer, by default None
:type elevation: float
:param elevation: Elevation in km, by default 0
:rtype: DepModel
:return: DepModel object
.. rubric:: Examples
>>> dep_range = np.arange(100)
>>> h = np.array([20, 15., 0])
>>> vp = np.array([5.8, 6.5, 8.04])
>>> vs = np.array([3.36, 3.75, 4.47])
>>> model = DepModel.read_layer_model(dep_range, h, vp, vs)
"""
mod = cls(dep_range, velmod=None, layer_mod=True, elevation=elevation)
if rho is not None:
mod.isrho = True
mod._elevation()
mod.model_array = _from_layer_model(mod.depths, h, vp, vs, rho=rho)
mod.vp, mod.vs, mod.rho = _layer2grid(mod.depths, mod.model_array)
return mod
def _elevation(self):
"""
set all depth related values:
1. depths_elev: depth array contains layer above sea level(represent by value lower than 0
2. depths_extend: depth array contains layer above sea level( evaluate from 0 to original depth
3. dz: 0, thick_1, thick_2....
4. thickness: thick_1, thick_2, ... , 0
requires elevation, depths_range or other component
.. rubric:: Examples
>>> model = DepModel(np.array([0, 20.1, 35.1, 100]))
>>> model.depths_elev
array([ 0. , 20.1, 35.1, 100. ])
>>> model.depths_extend
array([ 0. , 20.1, 35.1, 100. ])
>>> model.dz
array([ 0. , 20.1, 15. , 64.9])
>>> model.thickness
array([20.1, 15. , 64.9, 0. ])
>>> model = DepModel(np.array([0, 20.1, 35.1, 100]),elevation=10.)
>>> print(model.depths_elev)
[-10. 10.1 25.1 90. 123.33333333]
>>> print(model.depths_extend)
[ 0. 20.1 35.1 100. 133.33333333]
>>> print(model.dz)
[ 0. 20.1 15. 64.9 33.33333333]
>>> print(model.thickness)
[20.1 15. 64.9 33.33333333 0. ]
"""
if self.elevation == 0:
self.depths_elev = self.depths
self.depths_extend = self.depths
else:
depths_append = np.arange(self.depths[-1]+self.dep_val,
self.depths[-1]+self.dep_val+np.floor(self.elevation/self.dep_val+1), self.dep_val)
self.depths_extend = np.append(self.depths, depths_append)
self.depths_elev = np.append(self.depths, depths_append) - self.elevation
self.dz = np.append(0, np.diff(self.depths_extend))
self.thickness = np.append(np.diff(self.depths_extend), 0.)
self.R = 6371.0 - self.depths_elev
[docs]
def plot_model(self, show=True):
"""
plot model with matplotlib
:type show: bool, optional
:param show: whether to show the plot, by default True
"""
plt.style.use('bmh')
if self.isrho:
self.model_fig = plt.figure(figsize=(6, 6))
fignum = 2
else:
self.model_fig = plt.figure(figsize=(4, 6))
fignum = 1
self.model_ax = self.model_fig.add_subplot(1, fignum, 1)
self.model_ax.step(self.vp, self.depths, where='pre', label='Vp')
self.model_ax.step(self.vs, self.depths, where='pre', label='Vs')
self.model_ax.legend()
self.model_ax.set_xlabel('Velocity (km/s)')
self.model_ax.set_ylabel('Depth (km)')
self.model_ax.set_ylim([self.depths[0], self.depths[-1]])
self.model_ax.invert_yaxis()
if self.isrho:
self.rho_ax = self.model_fig.add_subplot(1, fignum, 2)
self.rho_ax.step(self.rho, self.depths, where='pre', color='C2', label='Density')
self.rho_ax.legend()
self.rho_ax.set_xlabel('Density (km/s)')
self.rho_ax.set_ylim([self.depths[0], self.depths[-1]])
self.rho_ax.invert_yaxis()
if show:
plt.show()
[docs]
def tpds(self, rayps, raypp, sphere=True):
# generate docstring
"""
calculate travel time of Pds
:type rayps: float or numpy.ndarray
:param rayps: ray parameter of Ps wave
:type raypp: float or numpy.ndarray
:param raypp: ray parameter of P wave
:type sphere: bool, optional
:param sphere: whether to use sphere earth, by default True
:rtype: numpy.ndarray
:return: travel time of Pds
"""
if sphere:
radius = self.R
else:
radius = 6371.
tps = np.cumsum((np.sqrt((radius / self.vs) ** 2 - rayps ** 2) -
np.sqrt((radius / self.vp) ** 2 - raypp ** 2)) *
(self.dz / radius))
return tps
[docs]
def tpppds(self, rayps, raypp, sphere=True):
"""
calculate travel time of Ppds
:type rayps: float or numpy.ndarray
:param rayps: ray parameter of Ps wave
:type raypp: float or numpy.ndarray
:param raypp: ray parameter of P wave
:type sphere: bool, optional
:param sphere: whether to use sphere earth, by default True
:rtype: numpy.ndarray
:return: travel time of Ppds
"""
if sphere:
radius = self.R
else:
radius = 6371.
tps = np.cumsum((np.sqrt((radius / self.vs) ** 2 - rayps ** 2) +
np.sqrt((radius / self.vp) ** 2 - raypp ** 2)) *
(self.dz / radius))
return tps
[docs]
def tpspds(self, rayps, sphere=True):
"""
calculate travel time of Ppsds
:type rayps: float or numpy.ndarray
:param rayps: ray parameter of Ps wave
:type sphere: bool, optional
:param sphere: whether to use sphere earth, by default True
:rtype: numpy.ndarray
:return: travel time of Ppsds
"""
if sphere:
radius = self.R
else:
radius = 6371.
tps = np.cumsum(2 * np.sqrt((radius / self.vs) ** 2 - rayps ** 2) *
(self.dz / radius))
return tps
[docs]
def radius_s(self, rayp, phase='P', sphere=True):
"""
calculate horizontal radius from piercing point to station postion.
P for Sp and S for Ps.
:type rayp: float or numpy.ndarray
:param rayp: ray parameter
:type phase: str, optional
:param phase: phase name, by default 'P'
:type sphere: bool, optional
:param sphere: whether to use sphere earth, by default True
:rtype: numpy.ndarray
:return: horizontal radius
.. rubric:: Examples
>>> model = DepModel(np.array([0, 20.1, 35.1, 100]))
>>> model.dz
array([ 0. , 20.1, 15. , 64.9])
>>> model.R
array([6371. , 6350.9, 6335.9, 6271. ])
>>> model.radius_s(1.2,phase="S", sphere=False)*111.2
array([0. , 0.0002478 , 0.00046823, 0.00142685])
"""
if phase == 'P':
vel = self.vp
else:
vel = self.vs
if sphere:
radius = self.R
else:
radius = 6371.
hor_dis = np.cumsum((self.dz / radius) / np.sqrt((1. / (rayp ** 2. * (radius / vel) ** -2)) - 1))
#hor_dis = np.sqrt((1. / (rayp ** 2. * (radius / vel) ** -2)) - 1)
return hor_dis
[docs]
def save_tvel(self, filename):
"""
save vel mod in tvel format for taup
:type filename: str
:param filename: output file name
.. rubric:: Examples
>>> model = DepModel(np.array([0, 20.1, 35.1, 100]))
>>> model.save_tvel("test")
0.00 5.80 3.36 2.72
20.10 5.800 3.36 2.72
20.10 6.50 3.75 2.92
35.10 6.500 3.75 2.92
35.10 8.04 4.47 3.32
100.00 8.040 4.47 3.32
100.00 8.05 4.49 3.36
"""
outlines =[]
depth = self.depths
if depth[0] != 0:
depth[0] = 0
for _i in range(self.depths.shape[0]):
out = f'{depth[_i]:>8.2f} {self.vp[_i]:>8.2f} {self.vs[_i]:>8.2f} {self.rho[_i]:>8.2f}'
outlines.append(out)
if _i == self.depths.shape[0] - 1:
break
out = f'{depth[_i+1]:>8.2f} {self.vp[_i]:>8.3f} {self.vs[_i]:>8.2f} {self.rho[_i]:>8.2f}'
outlines.append(out)
# this switch for doctest
if filename == 'test':
for _i in range(len(outlines)):
print(outlines[_i])
else:
try:
f = open(filename, 'w')
f.write("temp_mod-P\ntemp_mod-S")
for i in range(len(outlines)):
f.write(outlines[i])
f.close()
except IOError:
raise IOError('cannot write to {}'.format(filename))
[docs]
def raylength(self, rayp, phase='P', sphere=True):
"""
calculate ray length, P for Sp and S for Ps
:type rayp: float or numpy.ndarray
:param rayp: ray parameter
:type phase: str, optional
:param phase: phase name, by default 'P'
:type sphere: bool, optional
:param sphere: whether to use sphere earth, by default True
:rtype: numpy.ndarray
:return: ray length
"""
if phase == 'P':
vel = self.vp
else:
vel = self.vs
if sphere:
radius = self.R
else:
radius = 6371.
raylen = (self.dz * radius) / (np.sqrt(((radius / self.vs) ** 2) - (rayp ** 2)) * vel)
return raylen
[docs]
@classmethod
def ccp_model(cls, dep_range = np.array([0, 20.1, 35.1, 100]),
elevation = 0, layerd = False, **kwargs):
"""
import ccp configure and init DepModel object for time2depth convertion
if any parameters given is wrong, return a default DepModel object
there's 3 types of input is allowed:
1. mod3d, stla, stlo: for 3d model( need modification
2. modfolder, staname: for dir
3. mod: for single file
"""
if kwargs.get("mod3d", None):
stla = kwargs.pop("stla",None)
stlo = kwargs.pop("stlo", None)
# if stla and stlo is not given return iasp91 model
if not stla or not stlo:
return cls(dep_range, "iasp91", elevation, layerd)
try:
mod = np.load(kwargs["mod3d"])
dep_tmp = DepModel(dep_range,"iasp91", elevation, layerd)
dep_tmp.vp, dep_tmp.vs = interp_depth_model(mod, stla, stlo, dep_range)
return dep_tmp
except:
return cls(dep_range, "iasp91", elevation, layerd)
elif kwargs.get("modfolder", None):
try:
mod = _load_mod(kwargs.pop("modfolder", "./"), kwargs.pop("staname", None))
except:
mod = "iasp91"
finally:
return cls(dep_range, mod, elevation, layerd)
else:
mod = kwargs.pop("mod", "iasp91")
return cls(dep_range,mod, elevation, layerd)
def _load_mod(datapath, staname):
"""Load 1D velocity model files with suffix of ".vel". The model file should be including 3 columns with depth, vp and vs.
:param datapath: Folder name with 1D velocity model files.
:type datapath: string
:param staname: The station name as a part of file name of 1D velocity model files.
:type staname: string
"""
expresion = join(datapath, "*"+staname+"*.vel")
modfiles = glob.glob(expresion)
if len(modfiles) == 0:
raise FileNotFoundError("The model file of {} were not found.".format(expresion))
elif len(modfiles) > 1:
raise ValueError('More then 1 file were found as the expresion: {}'.format(expresion))
else:
return modfiles[0]
[docs]
def interp_depth_model(model, lat, lon, new_dep):
""" Interpolate Vp and Vs from 3D velocity with a specified depth range.
:param mod3d: 3D velocity loaded from a ``.npz`` file
:type mod3d: :meth:`np.lib.npyio.NpzFile`
:param lat: Latitude of position in 3D velocity model
:type lat: float
:param lon: Longitude of position in 3D velocity model
:type lon: float
:param new_dep: 1D array of depths in km
:type new_dep: :meth:`np.ndarray`
:rtype: :meth:`np.ndarray`
:return: Vp and Vs in ``new_dep``
"""
# model = np.load(modpath)
points = [[depth, lat, lon] for depth in new_dep]
vp = interpn((model['dep'], model['lat'], model['lon']), model['vp'], points, bounds_error=False, fill_value=None)
vs = interpn((model['dep'], model['lat'], model['lon']), model['vs'], points, bounds_error=False, fill_value=None)
return vp, vs
if __name__ == "__main__":
import doctest
doctest.testmod()