import numpy as np
from seispy.geo import km2deg, latlon_from, cosd, extrema, skm2srad, rad2deg
from seispy import distaz
from seispy.rfcorrect import DepModel
from seispy.setuplog import setuplog
from scikits.bootstrap import ci
from seispy.ccppara import ccppara, CCPPara
from seispy.signal import smooth
from seispy.utils import check_stack_val, read_rfdep
from scipy.interpolate import interp1d
import warnings
import sys
[docs]def gen_center_bin(center_lat, center_lon, len_lat, len_lon, val):
"""
Create spaced grid point with coordinates of the center point in the area in spherical coordinates.
:param center_lat: Latitude of the center point.
:type center_lat: float
:param center_lon: Longitude of the center point.
:type center_lon: float
:param len_lat: Half length in degree along latitude axis.
:type len_lat: float
:param len_lon: Half length in degree along longitude axis.
:type len_lon: float
:param val: Interval in degree between adjacent grid point.
:type val: float
:return: Coordinates of Grid points.
:rtype: 2-D ndarray of floats with shape (n, 2), where n is the number of grid points.
"""
lats = np.arange(0, 2*len_lat, val)
lons = np.arange(0, 2*len_lon, val)
plat, plon = latlon_from(center_lat, center_lon, 0, 90)
da = distaz(plat, plon, center_lat, center_lon)
begx = -len_lon
begy = -len_lat
bin_loca = []
bin_mat = np.zeros([lats.size, lons.size, 2])
bin_map = np.zeros([lats.size, lons.size]).astype(int)
n = 0
for j in range(lats.size):
delyinc = j * val + begy
delt = da.delta + delyinc
for i in range(lons.size):
azim = da.az + (begx + i * val) / cosd(delyinc)
glat, glon = latlon_from(plat, plon, azim, delt)
if glon > 180:
glon -= 360
bin_loca.append([glat, glon])
bin_mat[j, i, 0] = glat
bin_mat[j, i, 1] = glon
bin_map[j, i] = n
n += 1
return np.array(bin_loca), bin_mat, bin_map
[docs]def bin_shape(cpara):
if cpara.shape == 'rect':
raise ValueError('The shape of bins must be set to \'circle\' in ccp3d mode.')
if cpara.bin_radius is None:
depmod = DepModel(cpara.stack_range)
fzone = km2deg(np.sqrt(0.5*cpara.domperiod*depmod.vs*cpara.stack_range))
else:
fzone = np.ones_like(cpara.stack_range) * km2deg(cpara.bin_radius)
return fzone
[docs]def boot_bin_stack(data_bin, n_samples=3000):
warnings.filterwarnings("ignore")
data_bin = data_bin[~np.isnan(data_bin)]
count = data_bin.shape[0]
if count > 1:
if n_samples is not None:
cci = ci(data_bin, n_samples=n_samples)
else:
cci = np.array([np.nan, np.nan])
mu = np.nanmean(data_bin)
else:
cci = np.array([np.nan, np.nan])
mu = np.nan
return mu, cci, count
def _get_sta(rfdep):
return np.array([[sta['stalat'], sta['stalon']] for sta in rfdep])
def _sta_val(stack_range, radius):
dep_mod = DepModel(stack_range)
x_s = np.cumsum((dep_mod.dz / dep_mod.R) / np.sqrt((1. / (skm2srad(0.08) ** 2. * (dep_mod.R / dep_mod.vs) ** -2)) - 1))
dis = radius + rad2deg(x_s[-1]) + 0.5
return dis
[docs]class CCP3D():
def __init__(self, cfg_file=None, log=None):
"""Class for 3-D CCP stacking, Usually used to study mantle transition zone structure.
:param cfg_file: Path to configure file. If not defined a instance of CCP3D.cpara will be initialed, defaults to None
:type cfg_file: str, optional
:param log: A logger instance. If not defined, seispy.sutuplog.logger will be initialed, defaults to None
:type log: seispy.sutuplog.logger , optional
"""
if log is None:
self.logger = setuplog()
else:
self.logger = log
if cfg_file is None:
self.cpara = CCPPara()
elif isinstance(cfg_file, str):
self.load_para(cfg_file)
else:
raise ValueError('cfg_file must be str format.')
self.stack_data = []
self.good_410_660 = np.array([])
self.good_depth = np.array([])
self.bin_loca = None
self.bin_mat = None
self.bin_map = None
[docs] def load_para(self, cfg_file):
try:
self.cpara = ccppara(cfg_file)
except Exception as e:
self.logger.CCPlog('Cannot open configure file {}'.format(cfg_file))
raise FileNotFoundError('{}'.format(e))
try:
self.stack_mul = check_stack_val(self.cpara.stack_val, self.cpara.dep_val)
except Exception as e:
self.logger.CCPlog.error('{}'.format(e))
raise ValueError('{}'.format(e))
[docs] def read_rfdep(self):
self.logger.CCPlog.info('Loading RFdepth data from {}'.format(self.cpara.depthdat))
try:
self.rfdep = read_rfdep(self.cpara.depthdat)
except FileNotFoundError as e:
self.logger.CCPlog.error('{}'.format(e))
raise FileNotFoundError('Cannot open file of {}'.format(self.cpara.depthdat))
[docs] def initial_grid(self):
self.read_rfdep()
self.bin_loca, self.bin_mat, self.bin_map = gen_center_bin(*self.cpara.center_bin)
self.fzone = bin_shape(self.cpara)
self.stalst = _get_sta(self.rfdep)
self.dismin = _sta_val(self.cpara.stack_range, self.fzone[-1])
def _select_sta(self, bin_lat, bin_lon):
return np.where(distaz(bin_lat, bin_lon, self.stalst[:, 0], self.stalst[:, 1]).delta <= self.dismin)[0]
[docs] def stack(self):
"""Search conversion points falling within a bin and stack them with bootstrap method.
"""
for i, bin_info in enumerate(self.bin_loca):
boot_stack = {}
bin_mu = np.zeros(self.cpara.stack_range.size)
bin_ci = np.zeros([self.cpara.stack_range.size, 2])
bin_count = np.zeros(self.cpara.stack_range.size)
self.logger.CCPlog.info('{}/{} bin at lat: {:.3f} lon: {:.3f}'.format(i + 1, self.bin_loca.shape[0],
bin_info[0], bin_info[1]))
idxs = self._select_sta(bin_info[0], bin_info[1])
for j, dep in enumerate(self.cpara.stack_range):
idx = int(j * self.stack_mul + self.cpara.stack_range[0]/self.cpara.dep_val)
bin_dep_amp = np.array([])
for k in idxs:
stop_idx = np.where(self.rfdep[k]['stopindex'] >= idx)[0]
fall_idx = np.where(distaz(self.rfdep[k]['piercelat'][stop_idx, idx], self.rfdep[k]['piercelon'][stop_idx, idx],
bin_info[0], bin_info[1]).delta < self.fzone[j])[0]
bin_dep_amp = np.append(bin_dep_amp, self.rfdep[k]['moveout_correct'][stop_idx[fall_idx], idx])
bin_mu[j], bin_ci[j], bin_count[j] = boot_bin_stack(bin_dep_amp, n_samples=self.cpara.boot_samples)
boot_stack['bin_lat'] = bin_info[0]
boot_stack['bin_lon'] = bin_info[1]
boot_stack['mu'] = bin_mu
boot_stack['ci'] = bin_ci
boot_stack['count'] = bin_count
self.stack_data.append(boot_stack)
[docs] def save_stack_data(self, fname):
"""Save stacked data and parameters to local as a npz file. To load the file, please use data = np.load(fname, allow_pickle=True).
data['cpara'] is the parameters when CCP stacking.
data['stack_data'] is the result of stacked data.
:param fname: file name of stacked data
:type fname: str
"""
if not isinstance(fname, str):
self.logger.CCPlog.error('fname should be in \'str\'')
raise ValueError('fname should be in \'str\'')
np.savez(fname, cpara=self.cpara, stack_data=self.stack_data)
def _search_peak(self, tr, peak_410_min=380, peak_410_max=440, peak_660_min=630, peak_660_max=690):
tr = smooth(tr, half_len=4)
idx_all = extrema(tr)
idx_ex = np.where(tr[idx_all] > 0)[0]
#idx = idx_all[idx_ex]
peak = tr[idx_all[idx_ex]]
peak_depth = idx_all[idx_ex] * self.cpara.stack_val + self.cpara.stack_range[0]
idx_410 = np.where((peak_depth>peak_410_min) & (peak_depth<peak_410_max))[0]
try:
idx_410_max = idx_410[np.nanargmax(peak[idx_410])]
dep_410 = peak_depth[idx_410_max]
except:
dep_410 = np.nan
idx_660 = np.where((peak_depth>peak_660_min) & (peak_depth<peak_660_max))[0]
try:
idx_660_max = idx_660[np.nanargmax(peak[idx_660])]
dep_660 = peak_depth[idx_660_max]
except:
dep_660 = np.nan
return dep_410, dep_660
[docs] def search_good_410_660(self, peak_410_min=380, peak_410_max=440, peak_660_min=630, peak_660_max=690):
self.good_410_660 = np.zeros_like(self.bin_loca)
for i, boot_stack in enumerate(self.stack_data):
self.good_410_660[i, 0], self.good_410_660[i, 1] = self._search_peak(boot_stack['mu'], peak_410_min, peak_410_max, peak_660_min, peak_660_max)
[docs] def save_good_410_660(self, fname):
with open(fname, 'w') as f:
for i, good_peak in enumerate(self.good_410_660):
if np.isnan(good_peak[0]):
ci_410 = np.array([np.nan, np.nan])
count_410 = np.nan
else:
idx = int((good_peak[0]-self.cpara.stack_range[0]) / self.cpara.stack_val)
ci_410 = self.stack_data[i]['ci'][idx]
count_410 = self.stack_data[i]['count'][idx]
if np.isnan(good_peak[1]):
ci_660 = np.array([np.nan, np.nan])
count_660 = np.nan
else:
idx = int((good_peak[1]-self.cpara.stack_range[0]) / self.cpara.stack_val)
ci_660 = self.stack_data[i]['ci'][idx]
count_660 = self.stack_data[i]['count'][idx]
f.write('{:.3f} {:.3f} {:.0f} {:.4f} {:.4f} {:.0f} {:.0f} {:.4f} {:.4f} {:.0f}\n'.format(
self.bin_loca[i, 0], self.bin_loca[i, 1], good_peak[0], ci_410[0], ci_410[1], count_410,
good_peak[1], ci_660[0], ci_660[1], count_660))
[docs] @classmethod
def read_stack_data(cls, stack_data_path, cfg_file=None, good_depth_path=None, ismtz=False):
ccp = cls(cfg_file)
data = np.load(stack_data_path, allow_pickle=True)
ccp.stack_data = data['stack_data']
ccp.cpara = data['cpara'].any()
ccp.bin_loca, ccp.bin_mat, ccp.bin_map = gen_center_bin(*ccp.cpara.center_bin)
if good_depth_path is not None:
if ismtz:
ccp.good_410_660[:, 0] = np.loadtxt(good_depth_path, usecols=[2])
ccp.good_410_660[:, 0] = np.loadtxt(good_depth_path, usecols=[6])
else:
ccp.good_depth = np.loadtxt(good_depth_path, usecols=[2])
return ccp
[docs] def get_depth_err(self, type='std'):
moho_err = np.zeros([self.bin_loca.shape[0], 2])
self.logger.CCPlog.info('Computing errors of selected depth')
if self.good_depth.size == 0:
self.logger.CCPlog.error('Please load good depths before.')
sys.exit(1)
if np.isnan(self.stack_data['ci']).all() and type == 'ci':
self.logger.CCPlog.warning('No confidence intervals in stack data, using standard division instead.')
type = 'std'
for i, _ in enumerate(self.bin_loca):
if np.isnan(self.good_depth[i]):
moho_err[i, 0], moho_err[i, 1] = np.nan, np.nan
else:
idx = np.nanargmin(np.abs(self.cpara.stack_range-self.good_depth[i]))
mu = self.stack_data[i]['mu']
min_idxes = extrema(mu, opt='min')
try:
low_idx = min_idxes[np.max(np.where((min_idxes - idx) < 0)[0])]
up_idx = min_idxes[np.min(np.where((min_idxes - idx) > 0)[0])]
except:
moho_err[i, 0], moho_err[i, 1] = np.nan, np.nan
continue
if type == 'std':
cvalue = mu[idx] - 1.645 * np.std(mu[low_idx:up_idx+1])/np.sqrt(up_idx-low_idx+1)
elif type == 'ci':
cvalue = self.stack_data[i]['ci'][idx, 0]
else:
self.logger.error('Reference type should be in \'std\' and \'ci\'')
sys.exit(1)
moho_err[i, 0], moho_err[i, 1] = self._get_err(mu[low_idx:up_idx+1],
self.cpara.stack_range[low_idx:up_idx+1], cvalue)
return moho_err
def _get_err(self, tr, dep, cvalue):
result = np.array([])
for i, amp in enumerate(tr[:-1]):
if (amp <= cvalue < tr[i+1]) or (amp > cvalue >= tr[i+1]):
result = np.append(result, interp1d([amp, tr[i+1]],
[dep[i], dep[i+1]])(cvalue))
if len(result) == 2:
return result[0], result[1]
else:
return np.nan, np.nan
if __name__ == '__main__':
bin_loca = gen_center_bin(48.5, 100, 5, 8, km2deg(55))
with open('/workspace/WMHG_MTZ/ccp_results/bin_loca.dat', 'w') as f:
for binin in bin_loca:
f.write('{} {}\n'.format(binin[1], binin[0]))