import sys
sys.path.insert(0,"/home/ebiondi/research/projects/LongValley/DAS-taptest-widgets/Python")
import TapTestWdgts
import numpy as np
import pandas as pd
import os
from scipy.interpolate import interp1d
from scipy.signal import hilbert
from tqdm import tqdm
import datetime, pytz
import dateutil.parser
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
params = {
'image.interpolation': 'nearest',
'image.cmap': 'gray',
'savefig.dpi': 300, # to adjust notebook inline plot size
'axes.labelsize': 12, # fontsize for x and y labels (was 10)
'axes.titlesize': 12,
'font.size': 12,
'legend.fontsize': 12,
'xtick.labelsize': 12,
'ytick.labelsize': 12,
'text.usetex':False
}
matplotlib.rcParams.update(params)
import gpxpy
time_format = "%Y-%m-%dT%H:%M:%S+00:00"
# Getting event list from data time range
from obspy.clients.fdsn.client import Client
client = Client("IRIS")
import obspy
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
#Adding pyDAS location
pyDAS_path = "/home/elibird/research/packages/DAS-utilities/build/"
import os
os.environ['LD_LIBRARY_PATH'] = pyDAS_path
import sys
sys.path.insert(0,pyDAS_path)
sys.path.insert(0,'/home/elibird/research/packages/DAS-utilities/python')
import DASutils
# import importlib
# importlib.reload(DASutils)
# importlib.reload(TapTestWdgts)
#Reading channel locations, we are only using the south array in this work
array_df_north = pd.read_csv('/kuafu/DASdata/DASinfo/DAS_ChannelLocation/Mammoth_north_100km.csv')
array_df_south = pd.read_csv('/kuafu/DASdata/DASinfo/DAS_ChannelLocation/Mammoth_south_100km.csv')
good_chN=np.where(array_df_north['status']=='good')[0]
good_chS=np.where(array_df_south['status']=='good')[0]
good_arr2=array_df_north[array_df_north['status']=='good']
ch_lat2=good_arr2['latitude']
ch_lon2=good_arr2['longitude']
good_arr=array_df_south[array_df_south['status']=='good']
ch_lat=good_arr['latitude']
ch_lon=good_arr['longitude']
good_ch=np.where(array_df_south['status']=='good')[0]
import geopandas as gpd
# Reading DEM from npz file
filename = "data/LongValley-30m-DEM.npz"
DEM_npz = np.load(filename)
DEM_array = DEM_npz['data']
coord_box = DEM_npz['boundbox']
minLatdem = coord_box[3,1]
maxLatdem = coord_box[0,1]
minLondem = coord_box[0,0]
maxLondem = coord_box[1,0]
# Latitude axis 0
# Longitude axis 1
# first element top-left corner
DEM_array.shape
#DEM axes
lonaxdem=np.arange(minLondem,maxLondem,(maxLondem-minLondem)/DEM_array.shape[1])
lataxdem=np.arange(minLatdem,maxLatdem,(maxLatdem-minLatdem)/DEM_array.shape[0])[::-1]
import utm
#read catalog of events, get relevant parameters for said events
cat=pd.read_csv("data/catalog_data.csv")
_,_,zone_n,zone_lt = utm.from_latlon(0.5*(minLatdem+maxLatdem),0.5*(minLondem+maxLondem))
x,y,_,_=utm.from_latlon(np.array(cat['latitude']),np.array(cat['longitude']),zone_n,zone_lt)
z=np.array(cat['depth_km'])*1000
mags=np.array(cat['magnitude'])
eqLats=np.array(cat['latitude'])
eqLons=np.array(cat['longitude'])
x_p,y_p,_,_=utm.from_latlon(37.4,-118.4,zone_n,zone_lt)
dists=np.sqrt((x-x_p)**2+(y-y_p)**2)/1000
#make smaller catalog of useful events close to the array
housecat=cat[((cat['magnitude']>2.0)&(dists<30))|((cat['magnitude']>2.5)&(dists<50))]
np.save('data/eqinfo.npy',housecat)
events=np.array(housecat['event_id'])
eqSels=np.load('data/eqinfo.npy',allow_pickle=True)
if False:
#make bandpassed event data file, this event file is saved
events=np.array(housecat['event_id'])
count=-1
event_data=np.zeros((2147,2200,40))
for eventID in tqdm.tqdm(events):
count+=1
eqSel = cat[cat["event_id"] == eventID]
eqTime = dateutil.parser.parse(eqSel["event_time"].item())
data_folder="/kuafu/DASEventData/mammoth_south_100km/data/"
EventFile = data_folder+"%s.h5"%eventID
if os.path.exists(EventFile):
################################################
# Reading Event data
data, info = DASutils.read_das_event_data_h5(EventFile, ichan=good_ch, freqmin=1, freqmax=8)
dCh = info["dx"]
ot = info["begTime"]
dt = info["dt"]
fs = 1.0/info["dt"]
otAx = ot.timestamp() - eqTime.timestamp()
nt = data.shape[1]
nCh = data.shape[0]
timeAx = np.linspace(otAx, otAx+(nt-1)/fs, nt)
min_t = -2.0
max_t = 20.0
it_min = int((min_t-timeAx[0])/dt+0.5)
it_max = int((max_t-timeAx[0])/dt+0.5)
event_data[:,:,count]=data[7000:,it_min:it_max]
np.save('data/convphase_data_BP1-8Hz.npy',event_data)
import tqdm
#function to normalize times series
def divbysd(x):
return (x-np.mean(x))/np.std(x)
#function for plotting
def forceAspect(ax,aspect=1):
im = ax.get_images()
extent = im[0].get_extent()
ax.set_aspect(abs((extent[1]-extent[0])/(extent[3]-extent[2]))/aspect)
#function to read in picks from Phasenet-DAS
def read_TT(pick_df, good_ch, phase, eqTime):
"""Function to read traveltime from dataframe"""
nCh = good_ch.shape[0]
TT = np.zeros(nCh)
W = np.zeros(nCh)
perc_picked = 0.0
pick_idx = pick_df[pick_df["phase_type"]==phase]["channel_index"].to_numpy()
pick_wei = pick_df[pick_df["phase_type"]==phase]["phase_score"].to_numpy()
pick_TT = [dateutil.parser.parse(itime).timestamp()- eqTime.timestamp()
for itime in pick_df[pick_df["phase_type"]==phase]["phase_time"]]
pick_TT = np.array(pick_TT)
idxTemp = np.array([[idx,np.where(ch_id==pick_idx)[0][0]]
for idx, ch_id in enumerate(good_ch) if ch_id in pick_idx])
if len(idxTemp) > 0:
pick_IdxGood = idxTemp[:,0]
pick_IdxAll = idxTemp[:,1]
TT[pick_IdxGood] = pick_TT[pick_IdxAll]
W[pick_IdxGood] = pick_df[pick_df["phase_type"]==phase]["phase_score"].to_numpy()[pick_IdxAll]
perc_picked = float(np.where(W > 0.0)[0].shape[0])/nCh*100.0
return TT, W, perc_picked
#define moving average function
def moving_average(a, n=50):
ret = np.cumsum(a, dtype=float)
ret[n:] = ret[n:] - ret[:-n]
return ret[n - 1:] / n
#based on moving average, define stlt function, with short window of 12 samples, long window of 30 samples
def stlt(a,ns=12,nl=30):
st=moving_average(a, n=ns)[nl:]
lt=moving_average(a, n=nl)[:-ns]
return st/lt
#load event data
event_data=np.load('data/convphase_data_BP1-8Hz.npy')
evids=list()
for eqSel in eqSels:
evids.append(eqSel[0])
#get channel locations
latmean=np.mean(ch_lat[7000:])
lonmean=np.mean(ch_lon[7000:])
x_mid,y_mid,zn,zlt=utm.from_latlon(latmean,lonmean)
x_end,y_end,_,_=utm.from_latlon(ch_lat[7000:],ch_lon[7000:],zn,zlt)
x_end=np.array(x_end)/1000
y_end=np.array(y_end)/1000
eqSels=np.load('data/eqinfo.npy',allow_pickle=True)
#Make converted phase picks (amplitude and timing) and plot event with picks
if True:
from matplotlib.patches import Polygon
nbb=np.load('data/nbb.npy')
nbb=np.concatenate((nbb,[nbb[0,:]]))
eqLons=eqSels[:,2]
eqLats=eqSels[:,3]
import warnings
conv_picks=np.zeros((40, 2098, 5))
warnings.filterwarnings('ignore')
event_data=np.load('data/convphase_data_BP1-8Hz.npy')
eqSels=np.load('data/eqinfo.npy',allow_pickle=True)
for evno in tqdm.tqdm(range(40)):
#calculate sta/lta for an event and read in the P arrival times.
evid=eqSels[evno][0]
pick_df=pd.read_csv('/kuafu/users/elibird/PhaseNetPicks/OwensRiverConvPhase/picks_phasenet_das_raw/%s.csv'%evid)
eqSel=eqSels[np.where(np.array(evids)==evid)][0]
eqTime=dateutil.parser.parse(eqSel[1])
P_TT,WP,perc_pickedP=read_TT(pick_df, good_ch, 'P', eqTime)
P_TT[P_TT==0]=np.nan
P_TT=moving_average(P_TT)
data=event_data[:,:,evno]
datanorm=np.apply_along_axis(divbysd,1,data)
ev=abs(hilbert(datanorm))
ev_ma=np.apply_along_axis(moving_average,0,ev)
evstlt=np.apply_along_axis(stlt,1,ev)
evstlt_ma=np.apply_along_axis(moving_average,0,evstlt)
eqLon=eqSel[2]
eqLat=eqSel[3]
eq_x,eq_y,_,_=utm.from_latlon(eqLat,eqLon,zn,zlt)
eq_x=eq_x/1000
eq_y=eq_y/1000
dists=np.sqrt((eq_x-x_end)**2+(eq_y-y_end)**2)
dists_ma=moving_average(dists)
eqDep=eqSel[4]
Magtype=eqSel[6]
eqMag=eqSel[5]
max_t=20-6/100
min_t=-2+36/100
#Find the delay time between the P and conv. phases with varying bounds on the sta/lta
uppers=np.linspace(2,4,20)
lags=np.zeros((evstlt_ma.shape[0],len(uppers)))
for chan in range(evstlt_ma.shape[0]):
for i in range(20):
try:
upper=uppers[i]
lower=2*upper/3
P_ind=int((P_TT[7000+chan]*100)-min_t*100+0.5)
x_findstart=evstlt_ma[chan,(P_ind+10):]
first_ind=P_ind+10+np.where(x_findstart<lower)[0][0]
conv_ind=np.where(evstlt_ma[chan,first_ind:]>upper)[0][0]+first_ind
lag=(conv_ind-P_ind)/100
if lag>1:
lag=np.nan
lags[chan,i]=lag
except:
lags[chan,i]=np.nan
#considering these variable bounds get error bounds on delay time and amplitude ratio
delaymean=np.nanmean(lags,axis=1)
delaymean[np.sum(lags==lags,axis=1)<5]=np.nan
delaystd=np.nanstd(lags,axis=1)
delaystd[np.sum(lags==lags,axis=1)<5]=np.nan
ampl_rats=np.empty(ev_ma.shape[0])
ampl_rats[:]=np.nan
min_t_amp=-2
max_t_amp=20
for chan in range(ev_ma.shape[0]):
if delaymean[chan]==delaymean[chan]:
P_ind_amp=int((P_TT[7000+chan]*100)-min_t_amp*100+0.5)
conv_ind_amp=P_ind_amp+int(delaymean[chan]*100)
firstlen=5
Pamp=np.max(ev_ma[chan,(P_ind_amp-5):(P_ind_amp+20)])
convamp=np.max(ev_ma[chan,(conv_ind_amp-5):(conv_ind_amp+20)])
ampl_rats[chan]=Pamp/convamp
convval=delaymean+P_TT[7000:]
convval[np.arange(7025,len(P_TT)+25)<7800]=np.nan
delaystd[np.arange(7025,len(P_TT)+25)<7800]=np.nan
p_val=P_TT[7000:]/dists_ma
conv_picks[evno,:,0]=delaymean
conv_picks[evno,:,1]=np.sqrt((delaystd**2)+(0.06**2))
conv_picks[evno,:,2]=p_val
conv_picks[evno,:,3]=ampl_rats
#plot the event
min_t_p=-2
max_t_p=20
plt.rcParams.update({'font.size': 16})
fig, ax = plt.subplots(1,figsize=(10,8))
clipVal = np.percentile(np.absolute(datanorm), 95)
ax.imshow(datanorm.T,
extent=[7000, 7000+datanorm.shape[0], max_t_p, min_t_p],
aspect='auto', vmin=-clipVal, vmax=clipVal, cmap='seismic')
ax.plot(np.arange(7025,len(P_TT)+25),P_TT[7000:],color='k',label='P')
ax.plot(np.arange(7025,len(P_TT)+25),convval,color='orange',label='Ps')
ax.fill_between(np.arange(7025,len(P_TT)+25),convval-delaystd*2,convval+delaystd*2,alpha=0.25,color='orange')
ax.set_ylabel("Time from earthquake [s]",fontsize=20)
ax.set_xlabel("Channel number",fontsize=20)
ax.set_ylim((max_t_p,0))
if (evid=='nc73746481')|(evid=='nc73753116'):
ax.set_ylim((12.5,0))
ax.grid()
ax.legend(fontsize=15)
plt.savefig('figs/%s_conv.png'%evid)
plt.close()
np.save('data/convpicks.npy',conv_picks)
conv_picks=np.load('data/convpicks.npy')
100%|███████████████████████████████████████████| 40/40 [05:44<00:00, 8.60s/it]
#Example figure with STA/LTA at varying channel/time for a single event
evid='nc73731161'
pick_df=pd.read_csv('/kuafu/users/elibird/PhaseNetPicks/OwensRiverConvPhase/picks_phasenet_das_raw/%s.csv'%evid)
eqSel=eqSels[np.where(np.array(evids)==evid)][0]
evind=np.where(np.array(evids)==evid)[0][0]
eqTime=dateutil.parser.parse(eqSel[1])
P_TT,WP,perc_pickedP=read_TT(pick_df, good_ch, 'P', eqTime)
S_TT,WS,perc_pickedS=read_TT(pick_df, good_ch, 'S', eqTime)
data=event_data[:,:,evind]
datanorm=np.apply_along_axis(divbysd,1,data)
ev=abs(hilbert(datanorm))
evstlt=np.apply_along_axis(stlt,1,ev)
evstlt_ma=np.apply_along_axis(moving_average,0,evstlt)
clipVal = np.percentile(evstlt_ma, 99)
max_t=20-6/100
min_t=-2+36/100
fig, ax= plt.subplots(1,figsize=(10,8))
stltplt=ax.imshow(evstlt_ma.T,
extent=[7025, 7025+evstlt_ma.shape[0], max_t, min_t],
aspect='auto', vmax=clipVal, cmap=plt.get_cmap('Spectral'))
fig.colorbar(stltplt,ax=ax).set_label('STA/LTA',fontsize=15)
ax.set_ylim(15,5)
ax.plot(np.arange(7025,len(P_TT)-25),P_TT[7025:-25],'r',label='P arrival')
ax.plot(np.arange(7025,len(S_TT)-25),S_TT[7025:-25],'r--',label='S arrival')
ax.fill_between(np.arange(7025,len(P_TT)-25),P_TT[7025:-25],P_TT[7025:-25]+1,alpha=0.5,label='Search Range')
ax.legend(fontsize=15)
ax.set_xlabel('Channel No.',fontsize=15)
ax.set_ylabel('Time After Event (s)',fontsize=15)
plt.savefig('figs/stlt.png')
#We have manually picked which earthquakes are
ampevbool=np.zeros(len(eqSels))
delevbool=np.zeros(len(eqSels))
ampevs=np.load('data/ampevs.npy')
delevs=np.load('data/delevs.npy')
for iev, ev in enumerate(eqSels[:,0]):
if ev in ampevs:
ampevbool[iev]=1
if ev in delevs:
delevbool[iev]=1
delevbool=delevbool.astype('bool')
ampevbool=ampevbool.astype('bool')
fig, ax=plt.subplots(1,figsize=(8,6))
ax.plot(moving_average(np.arange(7000,len(P_TT))),conv_picks[delevbool,:,0].T,color='k',alpha=0.2)
meanpick=np.nanmean(conv_picks[delevbool,:,0],axis=0)
ppick=np.empty(conv_picks.shape[1])
pstdpick=np.empty(conv_picks.shape[1])
ppick[:]=np.nan
pstdpick[:]=np.nan
for i in range(conv_picks.shape[1]):
p4ave=conv_picks[delevbool,i,2][conv_picks[delevbool,i,0]==conv_picks[delevbool,i,0]]
if len(p4ave>=1):
ppick[i]=np.mean(p4ave)
pstdpick[i]=np.std(p4ave)
numpicks=np.sum(conv_picks[delevbool,:,0]==conv_picks[delevbool,:,0],axis=0)
stdpick=np.nanstd(conv_picks[delevbool,:,0],axis=0)
meanratpick=np.nanmean(conv_picks[ampevbool,:,3],axis=0)
stdratpick=np.nanstd(conv_picks[ampevbool,:,3],axis=0)
meanpickalt=np.nanmean(conv_picks[:,:,0],axis=0)
numpicksalt=np.sum(conv_picks[:,:,0]==conv_picks[:,:,0],axis=0)
meanpickalt[numpicksalt<3]=np.nan
meanpick[numpicks<3]=np.nan
ppick[numpicks<3]=np.nan
stdpick[numpicks<3]=np.nan
meanratpick[numpicks<3]=np.nan
stdratpick[numpicks<3]=np.nan
ax.plot(moving_average(np.arange(7000,len(P_TT))),meanpick)
ax.fill_between(moving_average(np.arange(7000,len(P_TT))),meanpick-stdpick,meanpick+stdpick,alpha=0.5)
ax.set_xlim(7850,9150)
ax.set_xlabel('Channel No.')
ax.set_ylabel('P to Converted Phase Delay')
plt.savefig('figs/delay_curve.png')
np.save('data/meanpick.npy',meanpick)
np.save('data/stdpick.npy',stdpick)
np.save('data/meanratpick.npy',meanratpick)
np.save('data/stdratpick.npy',stdratpick)
np.save('data/ppick.npy',ppick)
meanratpick=np.nanmean(conv_picks[ampevbool,:,3])
stdratpick=np.nanstd(conv_picks[ampevbool,:,3])
p_overall=np.mean(conv_picks[delevbool,:,2][conv_picks[delevbool,:,0]==conv_picks[delevbool,:,0]])
ratdat=np.array((meanratpick,stdratpick,p_overall))
np.save('data/ratdat.npy',ratdat)
plt.rcParams.update({'font.size': 15})
plt.tight_layout()
matplotlib.rc('xtick', labelsize=18)
matplotlib.rc('ytick', labelsize=18)
fig, ax1= plt.subplots(1,figsize=(13.5,10))
dem=DEM_array[np.argmin(abs(lataxdem-39.00)):np.argmin(abs(lataxdem-36.75)),np.argmin(abs(lonaxdem+120)):np.argmin(abs(lonaxdem+117.5))]
dem_small=DEM_array[np.argmin(abs(lataxdem-37.65)):np.argmin(abs(lataxdem-37.18)),np.argmin(abs(lonaxdem+118.72)):np.argmin(abs(lonaxdem+118.18))]
shading_style = {'vert_exag': 0.1,
'vmin': -3000,
'vmax': np.max(dem_small),
'azdeg': 315,
'altdeg': 25,
'cmap': plt.get_cmap('gray'),
'blend_mode': 'overlay'}
ls = matplotlib.colors.LightSource()
dem_rgb = ls.shade(dem, cmap=shading_style['cmap'], blend_mode=shading_style['blend_mode'], vmin=shading_style['vmin'], vmax=shading_style['vmax'])
demplot=ax1.imshow(dem_rgb,
extent=[-120,-117.5,36.75,39])
eqLons,eqLats=eqSels[:,2],eqSels[:,3]
ax1.plot(eqLons[delevbool],eqLats[delevbool],'r*',alpha=0.5,label='All Earthquakes Considered')
for ev_i in range(2):
evid=['nc73746481','nc73753116'][ev_i]
coluse=['violet','turquoise'][ev_i]
eqSel=eqSels[np.where(eqSels[:,0]==evid)[0][0]]
ax1.plot(eqSel[2], eqSel[3],'*', markersize=18,label='Event %s'%(ev_i+1),color='k',markerfacecolor=coluse)
ax1.plot(ch_lon,ch_lat,color='orange',label='Fiber')
ax1.plot(ch_lon[5999:],ch_lat[5999:],color='purple',label='Modeled Section')
nbbp=Polygon(nbb,alpha=0.5,label='NBB (Stevens et al., 2013)')
ax1.add_patch(nbbp)
ax1.set_xlim(-118.72,-118.18)
ax1.set_ylim(37.18,37.65)
#ax1.plot(np.concatenate((nbb[:,0],[nbb[0,0]])),np.concatenate((nbb[:,1],[nbb[0,1]])))
for c in np.arange(0,len(ch_lon),500):
if c>=5999:
ax1.plot(np.array(ch_lon)[c],np.array(ch_lat)[c],'o',color='purple')
else:
ax1.plot(np.array(ch_lon)[c],np.array(ch_lat)[c],'o',color='orange')
for c in np.arange(6000,len(ch_lon),1000):
ax1.text(np.array(ch_lon)[c]+0.005,np.array(ch_lat)[c]+0.005,c,color='purple')
ax1.grid()
ax1.set_xlabel('Longitude', fontsize=18)
ax1.set_ylabel('Latitude', fontsize=18)
plt.plot(-118.3952,37.3614,'s',color='darkblue',markersize=10)
plt.text(-118.3952-.085,37.3614-.007,'Bishop, CA',color='darkblue')
ax1.legend(fontsize=15)
clipVal = np.percentile(np.absolute(datanorm), 95)
import cartopy.crs as ccrs
import cartopy.feature as cfeature
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='50m',
facecolor='none')
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='50m',
facecolor='none')
extent_ax = [-118.72, -118.18, 37.18, 37.65]
box_lon = [extent_ax[i] for i in [0,1,1,0,0]]
box_lat = [extent_ax[i] for i in [2,2,3,3,2]]
ax2 = fig.add_axes([0.23, 0.15, 0.22, 0.3],projection=ccrs.PlateCarree())
ax2.set_extent([-125,-112,32,42.5])
ax2.plot(box_lon, box_lat, 'r-', lw=2)
ax2.add_feature(cfeature.LAND)
ax2.add_feature(cfeature.COASTLINE)
ax2.add_feature(cfeature.BORDERS)
ax2.add_feature(states_provinces, edgecolor='gray')
#ax2.outline_patch.set_linewidth(2)
aspect_ratio = 110.574/(111.320*np.cos(37.8/180.0*np.pi))
ax2.set_aspect(aspect_ratio)
#ax2.set_aspect(aspect_ratio)
plt.savefig('figs/map.png')
<Figure size 640x480 with 0 Axes>