提交 275299af 编辑于 作者: Mijian Xu's avatar Mijian Xu 😷
浏览文件

patch for 1.2.9

上级 20b64047
此差异已折叠。
......@@ -9,6 +9,7 @@ from seispy.psrayp import get_psrayp
from seispy.rfani import RFAni
from seispy.slantstack import SlantStack
import matplotlib.pyplot as plt
from seispy.utils import DepModel, tpds, Mod3DPerturbation
import warnings
import glob
......@@ -172,19 +173,6 @@ class RFStation(SACStation):
super().__init__(data_path, only_r=only_r)
class DepModel(object):
def __init__(self, YAxisRange, velmod='iasp91'):
VelocityModel = np.loadtxt(from_file(velmod))
self.depthsraw = VelocityModel[:, 0]
self.vpraw = VelocityModel[:, 1]
self.vsraw = VelocityModel[:, 2]
self.vp = interp1d(self.depthsraw, self.vpraw, bounds_error=False, fill_value='extrapolate')(YAxisRange)
self.vs = interp1d(self.depthsraw, self.vsraw, bounds_error=False, fill_value='extrapolate')(YAxisRange)
self.depths = YAxisRange
self.dz = np.append(0, np.diff(self.depths))
self.R = 6371 - self.depths
def _imag2nan(arr):
StopIndex = np.where(np.imag(arr) == 1)[0]
if StopIndex.size != 0:
......@@ -192,16 +180,6 @@ def _imag2nan(arr):
return arr
def from_file(mode_name):
if exists(mode_name):
filename = mode_name
elif exists(join(dirname(__file__), 'data', mode_name.lower()+'.vel')):
filename = join(dirname(__file__), 'data', mode_name.lower()+'.vel')
else:
raise ValueError('No such velocity mode')
return filename
def moveoutcorrect_ref(stadatar, raypref, YAxisRange, sampling=None, shift=None, velmod='iasp91'):
"""
:param stadatar: data class of SACStation
......@@ -226,10 +204,8 @@ def moveoutcorrect_ref(stadatar, raypref, YAxisRange, sampling=None, shift=None,
for i in range(stadatar.ev_num):
# x_s[i] = np.cumsum((dep_mod.dz / dep_mod.R) / np.sqrt((1. / (stadatar.rayp[i] ** 2. * (dep_mod.R / dep_mod.vs) ** -2)) - 1))
# x_p[i] = np.cumsum((dep_mod.dz / dep_mod.R) / np.sqrt((1. / (stadatar.rayp[i] ** 2. * (dep_mod.R / dep_mod.vp) ** -2)) - 1))
Tpds[i] = np.cumsum((np.sqrt((dep_mod.R / dep_mod.vs) ** 2 - stadatar.rayp[i] ** 2) -
np.sqrt((dep_mod.R / dep_mod.vp) ** 2 - stadatar.rayp[i] ** 2)) * (dep_mod.dz / dep_mod.R))
Tpds_ref = np.cumsum((np.sqrt((dep_mod.R / dep_mod.vs) ** 2 - raypref ** 2) -
np.sqrt((dep_mod.R / dep_mod.vp) ** 2 - raypref ** 2)) * (dep_mod.dz / dep_mod.R))
Tpds[i] = tpds(dep_mod, stadatar.rayp[i], stadatar.rayp[i])
Tpds_ref = tpds(dep_mod, raypref, raypref)
Newdatar = np.zeros([stadatar.ev_num, stadatar.rflength])
EndIndex = np.zeros(stadatar.ev_num)
......@@ -296,8 +272,7 @@ def psrf2depth(stadatar, YAxisRange, sampling, shift, velmod='iasp91', srayp=Non
for i in range(stadatar.ev_num):
x_s[i] = np.cumsum((dep_mod.dz / dep_mod.R) / np.sqrt((1. / (stadatar.rayp[i] ** 2. * (dep_mod.R / dep_mod.vs) ** -2)) - 1))
x_p[i] = np.cumsum((dep_mod.dz / dep_mod.R) / np.sqrt((1. / (stadatar.rayp[i] ** 2. * (dep_mod.R / dep_mod.vp) ** -2)) - 1))
Tpds[i] = np.cumsum((np.sqrt((dep_mod.R / dep_mod.vs) ** 2 - stadatar.rayp[i] ** 2) -
np.sqrt((dep_mod.R / dep_mod.vp) ** 2 - stadatar.rayp[i] ** 2)) * (dep_mod.dz / dep_mod.R))
Tpds[i] = tpds(dep_mod, stadatar.rayp[i], stadatar.rayp[i])
elif isinstance(srayp, str) or isinstance(srayp, np.lib.npyio.NpzFile):
if isinstance(srayp, str):
if not exists(srayp):
......@@ -311,8 +286,7 @@ def psrf2depth(stadatar, YAxisRange, sampling, shift, velmod='iasp91', srayp=Non
rayp = skm2srad(sdeg2skm(rayp))
x_s[i] = np.cumsum((dep_mod.dz / dep_mod.R) / np.sqrt((1. / (rayp ** 2. * (dep_mod.R / dep_mod.vs) ** -2)) - 1))
x_p[i] = np.cumsum((dep_mod.dz / dep_mod.R) / np.sqrt((1. / (stadatar.rayp[i] ** 2. * (dep_mod.R / dep_mod.vp) ** -2)) - 1))
Tpds[i] = np.cumsum((np.sqrt((dep_mod.R / dep_mod.vs) ** 2 - rayp ** 2) -
np.sqrt((dep_mod.R / dep_mod.vp) ** 2 - stadatar.rayp[i] ** 2)) * (dep_mod.dz / dep_mod.R))
Tpds[i] = tpds(dep_mod, rayp, stadatar.rayp[i])
else:
raise TypeError('srayp should be path to Ps rayp lib')
......@@ -360,8 +334,7 @@ def psrf_1D_raytracing(stadatar, YAxisRange, velmod='iasp91', srayp=None):
raylength_s[i] = (dep_mod.dz * dep_mod.R) / (np.sqrt(((dep_mod.R / dep_mod.vs) ** 2) - (stadatar.rayp[i] ** 2)) * dep_mod.vs)
x_p = np.cumsum((dep_mod.dz / dep_mod.R) / np.sqrt((1. / (stadatar.rayp[i] ** 2. * (dep_mod.R / dep_mod.vp) ** -2)) - 1))
raylength_p[i] = (dep_mod.dz * dep_mod.R) / (np.sqrt(((dep_mod.R / dep_mod.vp) ** 2) - (stadatar.rayp[i] ** 2)) * dep_mod.vp)
Tpds[i] = np.cumsum((np.sqrt((dep_mod.R / dep_mod.vs) ** 2 - stadatar.rayp[i] ** 2) -
np.sqrt((dep_mod.R / dep_mod.vp) ** 2 - stadatar.rayp[i] ** 2)) * (dep_mod.dz / dep_mod.R))
Tpds[i] = tpds(dep_mod, stadatar.rayp[i], stadatar.rayp[i])
pplat_s[i], pplon_s[i] = latlon_from(stadatar.stla, stadatar.stlo, stadatar.bazi[i], rad2deg(x_s))
pplat_p[i], pplon_p[i] = latlon_from(stadatar.stla, stadatar.stlo, stadatar.bazi[i], rad2deg(x_p))
elif isinstance(srayp, str) or isinstance(srayp, np.lib.npyio.NpzFile):
......@@ -379,8 +352,7 @@ def psrf_1D_raytracing(stadatar, YAxisRange, velmod='iasp91', srayp=None):
raylength_s[i] = (dep_mod.dz * dep_mod.R) / (np.sqrt(((dep_mod.R / dep_mod.vs) ** 2) - (rayp ** 2)) * dep_mod.vs)
x_p = np.cumsum((dep_mod.dz / dep_mod.R) / np.sqrt((1. / (stadatar.rayp[i] ** 2. * (dep_mod.R / dep_mod.vp) ** -2)) - 1))
raylength_p[i] = (dep_mod.dz * dep_mod.R) / (np.sqrt(((dep_mod.R / dep_mod.vp) ** 2) - (stadatar.rayp[i] ** 2)) * dep_mod.vp)
Tpds[i] = np.cumsum((np.sqrt((dep_mod.R / dep_mod.vs) ** 2 - rayp ** 2) -
np.sqrt((dep_mod.R / dep_mod.vp) ** 2 - stadatar.rayp[i] ** 2)) * (dep_mod.dz / dep_mod.R))
Tpds[i] = tpds(dep_mod, rayp, stadatar.rayp[i])
x_s = _imag2nan(x_s)
x_p = _imag2nan(x_p)
pplat_s[i], pplon_s[i] = latlon_from(stadatar.stla, stadatar.stlo, stadatar.bazi[i], rad2deg(x_s))
......@@ -476,30 +448,6 @@ def interp_depth_model(model, lat, lon, new_dep):
return vp, vs
class Mod3DPerturbation:
def __init__(self, modpath, YAxisRange, velmod='iasp91'):
dep_mod = DepModel(YAxisRange, velmod=velmod)
self.model = np.load(modpath)
new1dvp = interp1d(dep_mod.depthsraw, dep_mod.vpraw)(self.model['dep'])
new1dvs = interp1d(dep_mod.depthsraw, dep_mod.vsraw)(self.model['dep'])
new1dvp, _, _ = np.meshgrid(new1dvp, self.model['lat'], self.model['lon'], indexing='ij')
new1dvs, _, _ = np.meshgrid(new1dvs, self.model['lat'], self.model['lon'], indexing='ij')
self.dvp = (self.model['vp'] - new1dvp) / new1dvp
self.dvs = (self.model['vs'] - new1dvs) / new1dvs
self.cvp = dep_mod.vp
self.cvs = dep_mod.vs
def interpdvp(self, points):
dvp = interpn((self.model['dep'], self.model['lat'], self.model['lon']), self.dvp, points,
bounds_error=False, fill_value=None)
return dvp
def interpdvs(self, points):
dvs = interpn((self.model['dep'], self.model['lat'], self.model['lon']), self.dvs, points,
bounds_error=False, fill_value=None)
return dvs
def psrf_3D_migration(pplat_s, pplon_s, pplat_p, pplon_p, raylength_s, raylength_p, Tpds, YAxisRange, mod3d):
ev_num, _ = raylength_p.shape
timecorrections = np.zeros_like(raylength_p)
......
......@@ -21,8 +21,8 @@ class SlantStack():
self.ref_dis = 65
self.rayp_range = np.arange(-0.35, 0.35, 0.01)
self.tau_range = np.arange(0, 100)
self.tau = []
self.drayp = []
self.syn_tau = np.array([])
self.syn_drayp = np.array([])
def stack(self, ref_dis=None, rayp_range=None, tau_range=None):
if ref_dis is not None and isinstance(ref_dis, (int, float)):
......@@ -51,24 +51,18 @@ class SlantStack():
self.stack_amp += interp1d(self.time_axis, self.datar[i, :], fill_value='extrapolate')(tps)
self.stack_amp /= ev_num
# for j in range(self.rayp_range.shape[0]):
# for i in range(ev_num):
# self.datar[i, :] = self.datar[i, :] / np.max(np.abs(self.datar[i, :]))
# tps = self.tau_range - self.rayp_range[j] * (self.dis[i] - self.ref_dis)
# tmp[i, :] = interp1d(self.time_axis, self.datar[i, :], fill_value='extrapolate')(tps)
# self.stack_amp[j, :] = np.mean(tmp, axis=0)
def syn_tps(self, phase_list, model='iasp91'):
model = TauPyModel(model=model)
def syn_tps(self, phase_list, velmodel='iasp91', focal_dep=10):
model = TauPyModel(model=velmodel)
phase_list.insert(0, 'P')
arrs = model.get_travel_times(10, self.ref_dis, phase_list=phase_list)
arrs = model.get_travel_times(focal_dep, self.ref_dis, phase_list=phase_list)
p_arr = arrs[0].time
p_rayp = skm2sdeg(srad2skm(arrs[0].ray_param))
self.tau = [arr.time - p_arr for arr in arrs[1:]]
self.drayp = [p_rayp - skm2sdeg(srad2skm(arr.ray_param))for arr in arrs[1:]]
self.syn_tau = [arr.time - p_arr for arr in arrs[1:]]
self.syn_drayp = [p_rayp - skm2sdeg(srad2skm(arr.ray_param))for arr in arrs[1:]]
def plot(self, cmap='jet', xlim=None, vmin=None, vmax=None, figpath=None):
self.fig = plt.figure()
plt.style.use("bmh")
self.fig = plt.figure(figsize=(8,5))
self.ax = self.fig.add_subplot()
if vmin is None and vmax is None:
vmax = np.max(np.abs(self.stack_amp))
......@@ -83,8 +77,9 @@ class SlantStack():
raise TypeError('vmin and vmax must be in int or in float type')
im = self.ax.pcolor(self.tau_range, self.rayp_range, self.stack_amp,
vmax=vmax, vmin=vmin, cmap=cmap)
self.fig.colorbar(im, ax=self.ax)
self.ax.scatter(self.tau, self.drayp, color='k', marker='x')
# self.fig.colorbar(im, ax=self.ax)
self.ax.grid()
self.ax.scatter(self.syn_tau, self.syn_drayp, color='k', marker='x')
if xlim is not None and isinstance(xlim, (list, np.ndarray)):
self.ax.set_xlim(xlim)
self.ax.set_xlabel('Time (s)')
......
......@@ -2,6 +2,7 @@ from os.path import join, dirname, exists
from scipy.io import loadmat
from matplotlib.colors import ListedColormap
import numpy as np
from scipy.interpolate import interp1d, interpn
def load_cyan_map():
......@@ -22,6 +23,7 @@ def check_stack_val(stack_val, dep_val):
else:
raise ValueError('stack_val must be a multiple of dep_val')
def read_rfdep(path):
try:
return np.load(path, allow_pickle=True)
......@@ -29,4 +31,62 @@ def read_rfdep(path):
try:
return np.load(path+'.npy', allow_pickle=True)
except Exception as e:
raise FileNotFoundError('Cannot open file of {}'.format(path))
\ No newline at end of file
raise FileNotFoundError('Cannot open file of {}'.format(path))
class DepModel(object):
def __init__(self, YAxisRange, velmod='iasp91'):
VelocityModel = np.loadtxt(from_file(velmod))
self.depthsraw = VelocityModel[:, 0]
self.vpraw = VelocityModel[:, 1]
self.vsraw = VelocityModel[:, 2]
self.vp = interp1d(self.depthsraw, self.vpraw, bounds_error=False, fill_value='extrapolate')(YAxisRange)
self.vs = interp1d(self.depthsraw, self.vsraw, bounds_error=False, fill_value='extrapolate')(YAxisRange)
self.depths = YAxisRange
self.dz = np.append(0, np.diff(self.depths))
self.R = 6371 - self.depths
def from_file(mode_name):
if exists(mode_name):
filename = mode_name
elif exists(join(dirname(__file__), 'data', mode_name.lower()+'.vel')):
filename = join(dirname(__file__), 'data', mode_name.lower()+'.vel')
else:
raise ValueError('No such velocity mode')
return filename
def tpds(dep_mod, rayps, raypp):
return np.cumsum((np.sqrt((dep_mod.R / dep_mod.vs) ** 2 - rayps ** 2) -
np.sqrt((dep_mod.R / dep_mod.vp) ** 2 - raypp ** 2)) *
(dep_mod.dz / dep_mod.R))
def radius_s(dep_mod, rayps):
return np.cumsum((dep_mod.dz / dep_mod.R) /
np.sqrt((1. / (rayps ** 2. * (dep_mod.R / dep_mod.vs) ** -2)) - 1))
class Mod3DPerturbation:
def __init__(self, modpath, YAxisRange, velmod='iasp91'):
dep_mod = DepModel(YAxisRange, velmod=velmod)
self.model = np.load(modpath)
new1dvp = interp1d(dep_mod.depthsraw, dep_mod.vpraw)(self.model['dep'])
new1dvs = interp1d(dep_mod.depthsraw, dep_mod.vsraw)(self.model['dep'])
new1dvp, _, _ = np.meshgrid(new1dvp, self.model['lat'], self.model['lon'], indexing='ij')
new1dvs, _, _ = np.meshgrid(new1dvs, self.model['lat'], self.model['lon'], indexing='ij')
self.dvp = (self.model['vp'] - new1dvp) / new1dvp
self.dvs = (self.model['vs'] - new1dvs) / new1dvs
self.cvp = dep_mod.vp
self.cvs = dep_mod.vs
def interpdvp(self, points):
dvp = interpn((self.model['dep'], self.model['lat'], self.model['lon']), self.dvp, points,
bounds_error=False, fill_value=None)
return dvp
def interpdvs(self, points):
dvs = interpn((self.model['dep'], self.model['lat'], self.model['lon']), self.dvs, points,
bounds_error=False, fill_value=None)
return dvs
\ No newline at end of file
支持 Markdown
0% .
You are about to add 0 people to the discussion. Proceed with caution.
先完成此消息的编辑!
想要评论请 注册