import pickle
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.ticker import FuncFormatter, MultipleLocator
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from obspy import Stream, UTCDateTime
from obspy.signal.rotate import rotate2zne
from obspy.clients.filesystem.sds import Client as sdsClient
from scipy.signal import correlation_lags, correlate
from scipy.signal import hilbert
def corr_phaseshift(zaux, zaur, period, srate):
corrs_ph = correlate(zaux, zaur)
lags = correlation_lags(len(zaux), len(zaur))
corrs_ph /= np.max(corrs_ph)
maxshift = round(period*srate/2)
minshift = -maxshift
boolsh = np.logical_and(lags>=minshift, lags<=maxshift)
newcorr = corrs_ph[boolsh]
newlag = lags[boolsh]
argmaxc = np.argmax(newcorr)
best_lag = newlag[argmaxc]
phase_shift = best_lag*360/(period*srate)
return phase_shift, best_lag
def calc_vars(radial, vertical, period, ref_time):
odict = {'ph_shift': [], 'char_funs': [],
'times_cf': [], 'hvs': [], 'ccs': []}
phs = np.arange(0, 182, 2)
lags = np.unique((period*vertical.stats.sampling_rate*phs/360).astype(int))
for lag in lags:
ph = lag*360/(period*vertical.stats.sampling_rate)
wwz = vertical.copy()
wwr = radial.copy()
wwz.data = wwz.data[lag:]
if lag!=0:
wwr.data = wwr.data[:-lag]
time_arr = wwr.times() + (wwr.stats.starttime - ref_time)
wlen = int(period)
step = int(round(period/5, 0))
if step<1:
step = 1
if wlen<1:
wlen = 1
t0win = []
corrs = []
aux_zr = Stream(traces=[wwz, wwr])
for window in aux_zr.slide(window_length=wlen, step=step):
auxz = window.select(component='Z')[0].data
auxr = window.select(component='R')[0].data
norm = (np.sum(auxz**2)*np.sum(auxr**2))**0.5
fcc = obcorrelate(auxz, auxr, shift=0, demean=True, normalize=None)/norm
shift, valuecc = xcorr_max(fcc)
t0win.append(window[0].stats.starttime-ref_time)
corrs.append(valuecc)
return(odict)
def grid_data(inx, inzdata):
common_x = inx[0]
length = common_x.size
in_x = inx.copy()
in_zdata = inzdata.copy()
for it, etime in enumerate(in_x):
thisl = len(etime)
in_x[it] = np.concatenate([etime, np.nan*np.ones(length-thisl)])
in_zdata[it] = np.concatenate([in_zdata[it], np.nan*np.ones(length-thisl)])
return(in_x, in_zdata)
def shift_unc(xvals, yvals, zgrid, per, evname, thres=0.8, save=False):
percs = [16, 50, 84]
rd = 0.01
xt, yph = np.meshgrid(xvals, yvals)
newxt, newyph = xt.flatten(), yph.flatten()
points = np.vstack((newxt, newyph)).T
fig, ax = plt.subplots()
CS = ax.contour(xt, yph, zgrid, levels=[0.0, 0.4, 0.6, 0.8],
colors='gray', alpha=0.7)
maxval = np.nanmax(zgrid)
thres_unc = 0.95*maxval
CSun = ax.contour(xt, yph, zgrid, levels=[thres_unc],
colors='black', alpha=0.8)
p_un = CSun.collections[0].get_paths()[-1]
bool_sel = p_un.contains_points(points, radius=rd)
ph_min, ph_med, ph_max = np.percentile(newyph[bool_sel], percs)
xt_min, xt_med, xt_max = np.percentile(newxt[bool_sel], percs)
if maxval>=thres:
CSthres = ax.contour(xt, yph, zgrid, levels=[thres],
colors='red', alpha=0.8, linestyles='dashed')
for c, contour in enumerate(CSthres.collections[0].get_paths()):
if any(contour.contains_points(p_un.vertices)) or\
any(p_un.contains_points(contour.vertices)):
break
if thres_unc<thres:
new_bool = contour.contains_points(points, radius=rd)
ph_min, ph_med, ph_max = np.percentile(newyph[new_bool], percs)
xt_min, xt_med, xt_max = np.percentile(newxt[new_bool], percs)
inph = np.logical_and(newyph<=ph_max, newyph>=ph_min)
incon = contour.contains_points(points, radius=rd)
selpoints = points[np.logical_and(inph, incon)]
minTmax = [ np.nanpercentile(selpoints[selpoints.T[1]==auxph].T[0], [0, 100])
for auxph in np.unique(selpoints.T[1]) ]
_, med_minT, _ = np.percentile(np.asarray(minTmax).T[0], percs)
_, med_maxT, _ = np.percentile(np.asarray(minTmax).T[1], percs)
else:
print(f'Max value not exceeding the threshold [{maxval}/{thres}]')
med_minT = np.nan
med_maxT = np.nan
unc_dict = dict(min_ph=ph_min, med_ph=ph_med, max_ph=ph_max,
minT_thres=med_minT, maxT_thres=med_maxT,
maxchf_mint=xt_min, maxchf_medt=xt_med, maxchf_maxt=xt_max)
if save:
ax.clabel(CS, CS.levels, inline=True,
fmt=FuncFormatter(lambda y,_: '{:g}'.format(y)), fontsize=10)
ax.set_xlabel('Time [s]')
ax.set_ylabel('Phase shift [$^o$]')
dict_unc = dict(color='red', alpha=0.6, ls='--')
for phs in [ph_min, ph_med, ph_max]:
ax.axhline(phs, **dict_unc)
for xts in [med_minT, med_maxT, xt_min, xt_med, xt_max]:
ax.axvline(xts, **dict_unc)
ax.yaxis.set_major_locator(MultipleLocator(30))
ax.yaxis.set_minor_locator(MultipleLocator(10))
ax.grid(which='major', ls='--', alpha=0.5, color='gray')
ax.text(0.01, 0.98, 'T$_0$ = ' + f'{per} s \n{evname}', ha='left',
va='top', fontweight='bold', fontsize=10, transform=ax.transAxes)
ax.text(0.01, ph_med, '$\phi^*$ = ' + f'{ph_med:.1f}' + '$^o$', ha='right',
va='bottom', fontweight='bold', color='red',
fontsize=10, transform=ax.get_yaxis_transform())
fig.tight_layout()
fig.savefig(f'{evname}_{per}s_chf.png', dpi=200)
plt.close(fig)
return(unc_dict)
def shift_waves(auxz, auxr, blag, time_ref):
if blag>=0:
newz = auxz.data[blag:]
if blag==0:
newr = auxr.data
arr_time = auxr.times() + (auxr.stats.starttime - time_ref)
else:
newr = auxr.data[:-blag]
arr_time = auxr.times()[:-blag] + (auxr.stats.starttime - time_ref)
else:
newz = auxz.data[:blag]
newr = auxr.data[-blag:]
arr_time = auxz.times()[:blag] + (auxz.stats.starttime - time_ref)
unsh_zz = auxz.copy()
auxz.data = newz
auxr.data = newr
return unsh_zz, arr_time
