import numpy as np
from seispy.geo import km2deg, deg2km, latlon_from, \
geoproject, sind, rad2deg, skm2srad, \
geo2sph, sph2geo
from seispy.setuplog import SetupLog
from seispy.distaz import distaz
from seispy.core.depmodel import DepModel
from seispy.rf2depth_makedata import Station
from seispy.ccppara import ccppara, CCPPara
from scikits.bootstrap import ci
from seispy.ccp3d import boot_bin_stack
from seispy.utils import check_stack_val, read_rfdep
from scipy.interpolate import interp1d
from os.path import exists, dirname, basename, join
import sys
def _prof_range(lat, lon):
dis = [0]
for i in range(lat.size-1):
dis.append(distaz(lat[i], lon[i], lat[i+1], lon[i+1]).degreesToKilometers())
return np.cumsum(dis)
[docs]
def create_center_bin_profile(stations, val=5, method='linear'):
"""Create bins along a profile with given stations
:param stations: Stations along a profile
:type stations: :class:`seispy.rf2depth_makedata.Station`
:param val: The interval between two points in km
:type val: float
:param method: Method for interpolation
:type method: str
:return: The location of bins (bin_loca), and length between each bin and the start point (profile_range)
:rtype: (numpy.array, numpy.array)
"""
if not isinstance(stations, Station):
raise TypeError('Stations should be seispy.rf2depth_makedata.Station')
dis_sta = _prof_range(stations.stla, stations.stlo)
dis_inter = np.append(np.arange(0, dis_sta[-1], val), dis_sta[-1])
r, theta, phi = geo2sph(np.zeros(stations.stla.size), stations.stla, stations.stlo)
# t_po = np.arange(stations.stla.size)
# ip_t_po = np.linspace(0, stations.stla.size, bin_num)
theta_i = interp1d(dis_sta, theta, kind=method, bounds_error=False, fill_value='extrapolate')(dis_inter)
phi_i = interp1d(dis_sta, phi, kind=method, bounds_error=False, fill_value='extrapolate')(dis_inter)
_, lat, lon = sph2geo(r, theta_i, phi_i)
return lat, lon, dis_inter
def _line_proj(lat1, lon1, lat2, lon2):
daz = distaz(lat1, lon1, lat2, lon2)
az1_begin = (daz.baz - 90) % 360
az2_begin = (daz.baz + 90) % 360
az1_end = (daz.az - 90) % 360
az2_end = (daz.az + 90) % 360
return az1_begin, az2_begin, az1_end, az2_end, daz
def _fix_filename(filename, typ='dat'):
dname = dirname(filename)
if not exists(dname) and dname != '':
raise FileExistsError('internal error')
bname = basename(filename)
sp_name = bname.split('.')
if len(sp_name) == 1:
return filename + '.' + typ
if sp_name[-1] != typ:
nname = '.'.join(sp_name[0:-1]) + '.' + typ
return join(dirname(filename), nname)
else:
return filename
def _bin_shape(cpara):
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 init_profile(lat1, lon1, lat2, lon2, val):
""" Initial bins along a profile with given position of two points.
:param lat1: The latitude of the start point
:type lat1: float
:param lon1: The lontitude of the start point
:type lon1: float
:param lat2: The latitude of the end point
:type lat2: float
:param lon2: The lontitude of the end point
:type lon2: float
:param val: The interval between two points in km
:type val: float
:return: The location of bins (bin_loca), and length between each bin and the start point (profile_range)
The bin_loca is positions of bins with a numpy.array with two column. The profile_range is distance between bin center and the start point with an 1D numpy.array.
:rtype: (numpy.array, numpy.array)
"""
azi = distaz(lat1, lon1, lat2, lon2).baz
dis = distaz(lat1, lon1, lat2, lon2).delta
profile_range = np.arange(0, deg2km(dis), val)
lat_loca, lon_loca = latlon_from(lat1, lon1, azi, km2deg(profile_range))
bin_loca = np.zeros([lat_loca.shape[0], 2])
bin_loca = np.vstack((lat_loca, lon_loca)).T
return bin_loca, profile_range
[docs]
class CCPProfile():
def __init__(self, cfg_file=None, log:SetupLog=SetupLog()):
"""CCPProfile class for CCP stacking along a profile
:param cfg_file: Configure file for CCP stacking
:type cfg_file: str
:param log: Log class for logging
:type log: :class:`seispy.setuplog.SetupLog`
"""
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 = []
def _load_para(self, cfg_file):
try:
self.cpara = ccppara(cfg_file)
except Exception as e:
self.logger.CCPlog.error('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):
"""Read RFdepth file
:raises FileNotFoundError: Not Found RFdepth file
"""
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_profile(self):
"""Initialize bins of profile
"""
self.read_rfdep()
if exists(self.cpara.stack_sta_list):
self.stations = Station(self.cpara.stack_sta_list)
if self.cpara.adaptive:
if not exists(self.cpara.stack_sta_list):
raise ValueError('Adaptive binning depends on existence of {} order by stations'.format(
self.cpara.stack_sta_list))
lat, lon, self.profile_range = create_center_bin_profile(self.stations, self.cpara.slide_val)
self.bin_loca = np.vstack((lat, lon)).T
else:
self.bin_loca, self.profile_range = init_profile(*self.cpara.line, self.cpara.slide_val)
self.fzone = _bin_shape(self.cpara)
self._get_sta()
self._select_sta()
def _get_sta(self):
"""Read station info from rfdep
"""
self.staname = [sta['station'] for sta in self.rfdep]
self.stalst = np.array([[sta['stalat'], sta['stalon']] for sta in self.rfdep])
def _select_sta(self):
if exists(self.cpara.stack_sta_list):
self.logger.CCPlog.info('Use stacking stations in {}'.format(self.cpara.stack_sta_list))
self.idxs = []
for sta in self.stations.station:
try:
idx = self.staname.index(sta)
self.idxs.append(idx)
except ValueError:
self.logger.CCPlog.warning('{} does not in RFdepth structure'.format(sta))
if self.cpara.shape == 'rect' and self.cpara.adaptive == False:
self._pierce_project(self.rfdep[idx])
elif self.cpara.shape == 'circle':
if self.cpara.width is not None:
self.logger.CCPlog.warning('The \'width\' will be ignored, when circle bin was set')
dep_mod = DepModel(self.cpara.stack_range)
x_s = np.cumsum((dep_mod.dz / dep_mod.R) / np.sqrt((1. / (skm2srad(0.085) ** 2. * (dep_mod.R / dep_mod.vs) ** -2)) - 1))
dis = self.fzone[-1] + rad2deg(x_s[-1]) + 0.3
# self.idxs = self._proj_sta(dis)
self.idxs = []
for i, bin_info in enumerate(self.bin_loca):
self.idxs.append(np.where(distaz(bin_info[0], bin_info[1], self.stalst[:, 0], self.stalst[:, 1]).delta <= dis)[0])
elif self.cpara.width is not None and self.cpara.shape == 'rect':
self.logger.CCPlog.info('Select stations within {} km perpendicular to the profile'.format(self.cpara.width))
self.idxs = self._proj_sta(self.cpara.width)
for idx in self.idxs:
self._pierce_project(self.rfdep[idx])
if self.cpara.stack_sta_list != '':
self._write_sta()
else:
self.logger.CCPlog.error('Width of profile was not set, the bin shape of {} is invalid'.format(self.cpara.shape))
sys.exit(1)
def _write_sta(self):
with open(self.cpara.stack_sta_list, 'w') as f:
for idx in self.idxs:
f.write('{}\t{:.3f}\t{:.3f}\n'.format(self.staname[idx],
self.stalst[idx,0], self.stalst[idx, 1]))
def _proj_sta(self, width):
az1_begin, az2_begin, az1_end, az2_end, daz = _line_proj(*self.cpara.line)
az_sta_begin = distaz(self.stalst[:, 0], self.stalst[:, 1], self.cpara.line[0], self.cpara.line[1]).az
az_sta_end = distaz(self.stalst[:, 0], self.stalst[:, 1], self.cpara.line[2], self.cpara.line[3]).az
if 0 <= daz.baz < 90 or 270 <= daz.baz < 360:
tmp_idx_begin = np.where(((az1_begin < az_sta_begin)&(az_sta_begin < 360)) | ((0 < az_sta_begin)&(az_sta_begin < az2_begin)))[0]
else:
tmp_idx_begin = np.where((az1_begin < az_sta_begin)&(az_sta_begin < az2_begin))[0]
if 90 <= daz.az < 270:
tmp_idx_end = np.where(((az1_end < az_sta_end) & (az_sta_end < az2_end)))[0]
else:
tmp_idx_end = np.where(((az1_end < az_sta_end) & (az_sta_end < 360)) | ((0 < az_sta_end) & (az_sta_end < az2_end)))[0]
sta_daz = distaz(self.stalst[:, 0], self.stalst[:, 1], self.cpara.line[0], self.cpara.line[1])
sta_dis = sind(daz.baz-sta_daz.az) * sta_daz.degreesToKilometers()
proj_idx = np.where(np.abs(sta_dis) < width)[0]
tmp_idx = np.intersect1d(tmp_idx_begin, tmp_idx_end)
final_idx = np.intersect1d(tmp_idx, proj_idx).astype(int)
if not final_idx.any():
self.logger.CCPlog.error('No stations within the profile belt with width of {}'.format(width))
sys.exit(1)
return final_idx
def _pierce_project(self, rfsta):
rfsta['projlat'] = np.zeros_like(rfsta['piercelat'])
rfsta['projlon'] = np.zeros_like(rfsta['piercelon'])
for i, dep in enumerate(self.cpara.depth_axis):
rfsta['projlat'][:, i], rfsta['projlon'][:, i] = geoproject(rfsta['piercelat'][:, i], rfsta['piercelon'][:, i], *self.cpara.line)
[docs]
def stack(self):
"""Stack RFs in bins
"""
if self.cpara.shape == 'circle' or self.cpara.adaptive:
field_lat = 'piercelat'
field_lon = 'piercelon'
elif self.cpara.shape == 'rect':
field_lat = 'projlat'
field_lon = 'projlon'
else:
pass
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 from {:.2f} km at lat: {:.3f} lon: {:.3f}'.format(i + 1, self.bin_loca.shape[0], self.profile_range[i], bin_info[0], bin_info[1]))
if self.cpara.shape == 'circle' and not exists(self.cpara.stack_sta_list):
idxs = self.idxs[i]
else:
idxs = self.idxs
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][field_lat][stop_idx, idx], self.rfdep[k][field_lon][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['profile_dis'] = self.profile_range[i]
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, format='npz'):
"""If format is ``npz``, saving 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.
If format is ``dat`` the stacked data will be save into a txt file with 8 columns,
including ``bin_lat``, ``bin_lon``, ``profile_dis``, ``depth``, ``amp``, ``ci_low``, ``ci_high`` and ``count``.
- ``bin_lat`` and ``bin_lon`` represent the position of each bin;
- ``profile_dis`` represents the distance in km between each bin and the start point of the profile;
- ``depth`` represents depth of each bin;
- ``amp`` means the stacked amplitude;
- ``ci_low`` and ``ci_high`` mean confidence interval with bootstrap method;
- ``count`` represents stacking number of each bin.
:param format: Format for stacked data
:type format: str
"""
self.cpara.stackfile = _fix_filename(self.cpara.stackfile, format)
self.logger.CCPlog.info('Saving stacked data to {}'.format(self.cpara.stackfile))
if not isinstance(self.cpara.stackfile, str):
self.logger.CCPlog.error('fname should be in \'str\'')
raise ValueError('fname should be in \'str\'')
if format == 'npz':
np.savez(self.cpara.stackfile, cpara=self.cpara, stack_data=self.stack_data)
elif format == 'dat':
with open(self.cpara.stackfile, 'w') as f:
for i, bin in enumerate(self.stack_data):
for j, dep in enumerate(self.cpara.stack_range):
if dep == self.cpara.stack_range[-1]:
f.write(
'{:.4f}\t{:.4f}\t{:.4f}\t{:.2f} nan nan nan nan\n'.format(
bin['bin_lat'], bin['bin_lon'],
self.profile_range[i], dep
)
)
else:
f.write(
'{:.4f}\t{:.4f}\t{:.4f}\t{:.2f}\t{:.4f}\t{:.4f}\t{:.4f}\t{:d}\n'.format(
bin['bin_lat'], bin['bin_lon'],
self.profile_range[i],
dep, bin['mu'][j], bin['ci'][j, 0],
bin['ci'][j, 1], int(bin['count'][j])
)
)
if __name__ == '__main__':
bin_loca, _ = init_profile(27.5, 94, 36.5, 92, 5)
print(bin_loca.shape)