In [1]:
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import numpy as np
import sys
sys.path.append('../zephyr/')
from pygeo.segyread import SEGYFile
from pygeo.analysis import energyratio
from zephyr.middleware import FullwvDatastore
%matplotlib inline
matplotlib.rcParams['savefig.dpi'] = 300
In [2]:
# Set parameters
projnm = 'toyani' # Project name
pickOffset = -5
# Parameters for STA/LTA
windowLength = 20 # How long a window to use
zeroStart = 50 # How far back to blank in the first pass
zeroRel = -10 # How far back from the initial pick to blank in the second pass
# How far back to shift picks from peak STA/LTA
erPickOffset = pickOffset - windowLength
# Parameters for inflection picker
inflectionSign = +1 # Final pick on + or - inflection
inflHalfLength = 10 # Half-width of search window for inflection picker
# How far back to shift picks from inflection
inflPickOffset = pickOffset
unit = 1e3 # 1000 ms / s
outputFormat = '%(index)8d %(tag)4d %(sx)0.6e %(sz)0.6e %(rx)0.6e %(rz)0.6e %(pick)e\n'
banner = """Picks output for %s by simplepicker
https://github.com/uwoseis/simplepicker/
"""%(projnm,)
In [3]:
fds = FullwvDatastore(projnm)
systemConfig = fds.systemConfig
infile = '%s.ppre'%(projnm,)
outfile = '%s.til'%(projnm,)
nsrc = fds.ini['ns']
nrec = fds.ini['nr']
sf = SEGYFile(infile)
data = sf[:]
er = energyratio(data, windowLength)
er0 = er.copy()
er0[:,:zeroStart] = 0
picks0 = np.argmax(er0, axis=1) + erPickOffset
del(er0)
er1 = er.copy()
for itrace in xrange(sf.ntr):
er1[itrace,:picks0[itrace]+zeroRel] = 0
picks1 = np.argmax(er1, axis=1) + erPickOffset
del(er1)
der2 = np.diff(data)
gathers = np.array_split(der2, nsrc)
subpicks = np.array_split(picks1, nsrc)
inflector = np.argmax if inflectionSign > 0 else np.argmin
findInflection = lambda traces, picks: [pick + inflector(trace[pick-inflHalfLength:pick+inflHalfLength+1])-inflHalfLength for pick, trace in zip(picks, traces)]
picks2 = np.concatenate(map(findInflection, gathers, subpicks)) + inflPickOffset
clipScale = 0.5*abs(data).max()
def getInfo(isrc):
trStart = isrc * nrec
trEnd = trStart + nrec
return data[trStart:trEnd], er[trStart:trEnd], picks0[trStart:trEnd], picks1[trStart:trEnd], picks2[trStart:trEnd]
In [4]:
shotGather = 10
# 72a-[esc]
fig = plt.figure()
plotOpts = {
}
def makePlot(title, picks):
plt.imshow(picks.reshape((nsrc, nrec)), **plotOpts)
plt.hlines(shotGather, 0, nrec, color='w')
plt.xlabel('Receiver No.')
plt.ylabel('Source No.')
plt.title(title)
ax = plt.subplot(1,3,1)
makePlot('0: STA/LTA', picks0)
plt.subplot(1,3,2)
makePlot('1: STA/LTA', picks1)
plt.subplot(1,3,3)
makePlot('2: Inflection', picks2)
fig.tight_layout()
# 72a-[esc]
dataGather, erGather, picks0Gather, picks1Gather, picks2Gather = getInfo(shotGather)
fig = plt.figure()
plt.subplot(1,2,1)
plotOpts = {
'vmin': -clipScale,
'vmax': clipScale,
'cmap': cm.bwr,
}
plt.imshow(dataGather.T, **plotOpts)
axis0 = plt.axis()
plt.plot(picks0Gather, 'm-')
plt.plot(picks1Gather, 'g-')
plt.plot(picks2Gather, 'k-')
plt.axis(axis0)
plt.subplot(1,2,2)
plotOpts = {
'cmap': cm.jet,
}
plt.imshow(erGather.T, **plotOpts)
plt.plot(picks0Gather, 'm-')
plt.plot(picks1Gather, 'g-')
plt.plot(picks2Gather, 'w-')
axis0 = plt.axis(axis0)
In [5]:
timePicks = picks2 * fds.ini['deltatt'] / unit
if fds.ini['reduce']:
reducedTime = True
tbegin = fds.ini['tbegin']
vel = fds.ini['vmin']
else:
reducedTime = False
sx, sz, sw = fds.ini['srcs'].T
rx, rz, rw = fds.ini['recs'].T
offsets = np.sqrt((sx.reshape((nsrc,1)) - rx.reshape((1,nrec)))**2 + (sz.reshape((nsrc,1)) - rz.reshape((1,nrec)))**2).ravel()
shifts = tbegin + offsets / vel
with open(outfile, 'w') as fp:
fp.write(banner)
fp.write('%d\n%d\n'%(len(timePicks),7))
for isrc in xrange(nsrc):
for irec in xrange(nrec):
index = isrc*nrec + irec
if reducedTime:
pick = unit * (timePicks[index] + shifts[index])
else:
pick = unit * timePicks[index]
data = {
'index': index + 1,
'tag': 0,
'sx': sx[isrc],
'sz': sz[isrc],
'rx': rx[irec],
'rz': rz[irec],
'pick': pick,
}
line = outputFormat%data
fp.write(line)
In [ ]: