Source code for seispy.ccpprofile

import numpy as np
import seispy
from seispy.geo import km2deg, deg2km, latlon_from, geoproject, sind, rad2deg, skm2srad
from seispy.setuplog import setuplog
from seispy.distaz import distaz
from seispy.rfcorrect 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, create_center_bin_profile
from os.path import exists, dirname, basename, join


[docs]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
[docs]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
[docs]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=None): 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 = []
[docs] 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.width is None and self.cpara.shape == 'circle': 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) self._write_sta() else: if self.cpara.width is not None: self.logger.CCPlog.error('Width of profile was set, the bin shape of {} is invalid'.format(self.cpara.shape)) raise ValueError('Width of profile was set, the bin shape of {} is invalid'.format(self.cpara.shape)) else: self.logger.CCPlog.error('Width of profile was not set, the bin shape of {} is invalid'.format(self.cpara.shape)) raise ValueError('Width of profile was not set, the bin shape of {} is invalid'.format(self.cpara.shape)) 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])) self._pierce_project(self.rfdep[idx]) 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)) raise ValueError('Satisfied stations not found') 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.width is None and self.cpara.shape == 'circle': 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. where 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)