seismic amplitude extraction¶

filename='3d_nearstack.sgy'
scan = segy_header_scan(filename, max_traces_scan=0)
with pd.option_context('display.max_rows', 89):
    display(scan[scan['mean'] > 0])
#     display(scan)

bytelocs = dict(iline=41, xline=21, cdpx=73, cdpy=77)
near = segy_loader(filename, **bytelocs)
near = near.assign_coords(iline=near.iline/1000, twt=near.twt+1500)

print('Inlines: {:.0f}, min={:.0f}, max={:.0f}'.format(near.iline.size, near.iline.values.min(), near.iline.values.max()))
print('Crosslines: {:.0f}, min={:.0f}, max={:.0f}'.format(near.xline.size, near.xline.values.min(), near.xline.values.max()))
print('TWT: {:.0f}, min={:.0f}, max={:.0f}'.format(near.twt.size, near.twt.values.min(), near.twt.values.max()))
print('Total number of traces: {}'.format(near.iline.size*near.xline.size))

display(ui, tool)

val = 1376

iline = near.data.sel(iline=val, twt=slice(1800, 2400))
opt = dict(y='twt', add_colorbar=True,
        interpolation='spline16', robust=True,
        yincrease=False, cmap='Greys')

f, ax = plt.subplots(figsize=(10,6))
iline.plot.imshow(ax=ax, x='xline', **opt)
x, t = hrz[hrz.il==val]['xl'], hrz[hrz.il==val]['twt']
# x, t = hrz.sel(iline=val).xline, hrz.sel(iline=val).twt
ax.plot(x, t, color='r')

seis = near.data

@interact(ii=(1300,1500,24),cc=(1500,2000,48),zz=(1800,2496,48))
def seisview3d(ii=1372, cc=1788, zz=2136):
    clip=abs(np.percentile(seis, 0.80))   

    ii_idx=inl.tolist().index(ii)
    cc_idx=crl.tolist().index(cc)
    zz_idx=twt.tolist().index(zz)
    
    f, ax = plt.subplots(1,3, figsize=(14,6))
    ax[0].imshow(seis[:,ii_idx,:],extent=[crl.min(),crl.max(),twt.max(),twt.min()],cmap='Greys',vmax=clip,vmin=-clip,aspect='auto')
    ax[0].plot(hrz[hrz['il']==ii]['xl'],hrz[hrz['il']==ii]['z'],color='r') 
    ax[0].axhline(zz,color='y',lw=2)
    ax[0].axvline(cc,color='y',lw=2)   
    ax[0].set_ylim(twt.max(),1800)
    ax[0].set_title('IL {0}'.format(ii))
    ax[0].set_xlabel('XL'), ax[0].set_ylabel('TWT [ms]')

    ax[1].imshow(seis[:,:,cc_idx],extent=[inl.min(),inl.max(),twt.max(),twt.min()],cmap='Greys',vmax=clip,vmin=-clip,aspect='auto')
    ax[1].plot(hrz[hrz['xl']==cc]['il'],hrz[hrz['xl']==cc]['z'],color='r')
    ax[1].axhline(zz,color='y',lw=2)
    ax[1].axvline(ii,color='y',lw=2)
    ax[1].set_ylim(twt.max(),1800)
    ax[1].set_title('XL {0}'.format(cc))
    ax[1].set_xlabel('IL'), ax[1].set_ylabel('TWT [ms]')

    ax[2].imshow(seis[zz_idx,:,:],extent=[inl.min(),inl.max(),crl.min(),crl.max()],cmap='Greys',vmax=clip,vmin=-clip,aspect='auto')
    ax[2].axhline(cc,color='y',lw=2)
    ax[2].axvline(ii,color='y',lw=2)
    ax[2].set_title('TWT {0}'.format(zz))
    ax[2].set_xlabel('IL'), ax[2].set_ylabel('XL')
    f.subplots_adjust(wspace=.4)

from ipywidgets import interact
    
@interact(ii=(1300,1500,24),cc=(1500,2000,48),zz=(1800,2496,48))
def seisview3d(ii=1372,cc=1788,zz=2136):
    clip=abs(np.percentile(seis, 0.80))   

    ii_idx=inl.tolist().index(ii)
    cc_idx=crl.tolist().index(cc)
    zz_idx=twt.tolist().index(zz)
    
    f, ax = plt.subplots(1,3, figsize=(14,6))

    ax[0].imshow(seis[:,ii_idx,:],extent=[crl.min(),crl.max(),twt.max(),twt.min()],cmap='Greys',vmax=clip,vmin=-clip,aspect='auto')
    ax[0].plot(hrz[hrz['il']==ii]['xl'],hrz[hrz['il']==ii]['z'],color='r') 
    ax[0].axhline(zz,color='y',lw=2)
    ax[0].axvline(cc,color='y',lw=2)   
    ax[0].set_ylim(twt.max(),1800)
    ax[0].set_title('IL {0}'.format(ii))
    ax[0].set_xlabel('XL'), ax[0].set_ylabel('TWT [ms]')

    ax[1].imshow(seis[:,:,cc_idx],extent=[inl.min(),inl.max(),twt.max(),twt.min()],cmap='Greys',vmax=clip,vmin=-clip,aspect='auto')
    ax[1].plot(hrz[hrz['xl']==cc]['il'],hrz[hrz['xl']==cc]['z'],color='r')
    ax[1].axhline(zz,color='y',lw=2)
    ax[1].axvline(ii,color='y',lw=2)
    ax[1].set_ylim(twt.max(),1800)
    ax[1].set_title('XL {0}'.format(cc))
    ax[1].set_xlabel('IL'), ax[1].set_ylabel('TWT [ms]')

    ax[2].imshow(seis[zz_idx,:,:],extent=[inl.min(),inl.max(),crl.min(),crl.max()],cmap='Greys',vmax=clip,vmin=-clip,aspect='auto')
    ax[2].axhline(cc,color='y',lw=2)
    ax[2].axvline(ii,color='y',lw=2)
    ax[2].set_title('TWT {0}'.format(zz))
    ax[2].set_xlabel('IL'), ax[2].set_ylabel('XL')
    f.subplots_adjust(wspace=.4)

import time
from scipy.interpolate import interp1d
from scipy.interpolate import splev, splrep

hrz_sub=hrz[hrz['il']==1376]
twt_finer=np.arange(hrz['z'].min(), hrz['z'].max(),0.1) # creates twt scale at 0.1 ms scale

amp0=np.zeros((hrz_sub.size,2))
start = time.time()
for i in range(hrz_sub.size):
    ii_idx=inl.tolist().index(hrz_sub['il'][i])
    cc_idx=crl.tolist().index(hrz_sub['xl'][i])  
    zz_idx = np.abs(twt-hrz_sub['z'][i]).argmin()
    amp0[i,0] = hrz_sub['xl'][i]
    amp0[i,1] = seis[zz_idx, ii_idx, cc_idx]
end = time.time()
elapsed_raw = end - start

amp1=np.zeros((hrz_sub.size,2))
start = time.time()
for i in range(hrz_sub.size):
    ii_idx=inl.tolist().index(hrz_sub['il'][i])
    cc_idx=crl.tolist().index(hrz_sub['xl'][i])  
    trace = seis[:, ii_idx, cc_idx].flatten()
    ff = interp1d(twt, trace, kind='cubic')
    amp1[i,0] = hrz_sub['xl'][i]
    amp1[i,1] = ff(hrz_sub['z'][i])
end = time.time()
elapsed_interp1d = end - start

amp2=np.zeros((hrz_sub.size,2))
start = time.time()
for i in range(hrz_sub.size):
    ii_idx=inl.tolist().index(hrz_sub['il'][i])
    cc_idx=crl.tolist().index(hrz_sub['xl'][i])  
    trace = seis[:, ii_idx, cc_idx].flatten()
    ff = splrep(twt, trace)
    amp2[i,0] = hrz_sub['xl'][i]
    amp2[i,1] = splev(hrz_sub['z'][i], ff)
end = time.time()
elapsed_splev = end - start

print('Time to extract amplitudes at closest sample: {:.2f} seconds.'.format(elapsed_raw))
print('Time to extract amplitudes by interpolation with interp1d: {:.2f} seconds.'.format(elapsed_interp1d))
print('Time to extract amplitudes by interpolation with splrep/splev: {:.2f} seconds.'.format(elapsed_splev))

diff = amp1[:,1] - amp2[:,1]
print('Average difference between interp1d and splrep/splev: {0}'.format(np.mean(diff)))

plt.figure(figsize=(15,5))
plt.plot(amp0[:,0], amp0[:,1], '.-k', lw=2, label='raw')
plt.plot(amp2[:,0], amp2[:,1], '.-r', lw=2, alpha=0.5, label='interpolated')
plt.legend()
plt.xlabel('XL')
plt.ylabel('amplitude')
plt.xlim(1600,1900)

hrz_extr = np.zeros((hrz.size,6))
twt_finer = np.arange(hrz['z'].min(), hrz['z'].max(),0.1) # creates twt scale at 0.1 ms scale
window = 8 # 8 samples = 32 ms

start = time.time()
for i in range(hrz.size):
    ii_idx=inl.tolist().index(hrz['il'][i])
    cc_idx=crl.tolist().index(hrz['xl'][i])  
    zz_idx = np.abs(twt-hrz['z'][i]).argmin()
    amp_raw = seis[zz_idx, ii_idx, cc_idx].flatten()
    amp_rms = np.sqrt(np.mean((np.square(seis[zz_idx-int(window/2):zz_idx+int(window/2), ii_idx, cc_idx]))))
    
    trace = seis[:, ii_idx, cc_idx].flatten()
    ff = splrep(twt, trace)
    amp_splev = splev(hrz['z'][i], ff)

    hrz_extr[i,0] = hrz['il'][i]
    hrz_extr[i,1] = hrz['xl'][i]
    hrz_extr[i,2] = hrz['z'][i]
    hrz_extr[i,3] = amp_raw
    hrz_extr[i,4] = amp_splev
    hrz_extr[i,5] = amp_rms
end = time.time()
elapsed = end - start

print('Time to extract amplitudes on 3D cube: {:.2f} seconds.'.format(elapsed))

clip_max=np.max(hrz_extr[:,3:5])
clip_min=np.min(hrz_extr[:,3:5])

f, ax = plt.subplots(1,3,figsize=(15,5))
map0 = ax[0].pcolormesh(X, Y, Araw, vmin=clip_min, vmax=clip_max, cmap='cubehelix')
map1 = ax[1].pcolormesh(X, Y, Aspl, vmin=clip_min, vmax=clip_max, cmap='cubehelix')
map2 = ax[2].pcolormesh(X, Y, Araw-Aspl, vmin=clip_min, vmax=clip_max, cmap='cubehelix')
ax[0].set_title('raw amplitude')
ax[1].set_title('interpolated amplitude')
ax[2].set_title('difference raw-interpolated')
for aa in ax.flatten():
    aa.set_xlabel('IL'), aa.set_ylabel('XL')
f.subplots_adjust(wspace=.3)
cax = f.add_axes([0.96, 0.2, 0.025, 0.6])
f.colorbar(map0, cax=cax, orientation='vertical')

plt.figure()
plt.pcolormesh(X, Y, Araw, vmin=clip_min, vmax=clip_max, cmap='cubehelix', alpha=0.5)
for r,val in enumerate(inl[10::20]):
    plt.axvline(val, color='k', ls='--', lw=2)
plt.colorbar()
plt.axes().set_aspect(0.4)
plt.xlabel('IL'), plt.ylabel('XL')
plt.title('selected inlines')
    
f, ax = plt.subplots(inl[10::20].size,1,figsize=(15,15))
for r,val in enumerate(inl[10::20]):
    crosslines = hrz_extr[hrz_extr[:,0]==val][:,1]
    amp_raw    = hrz_extr[hrz_extr[:,0]==val][:,3]
    amp_splev  = hrz_extr[hrz_extr[:,0]==val][:,4]

    ax[r].plot(crosslines, amp_raw, '-k', lw=2, label='raw')
    ax[r].plot(crosslines, amp_splev, '-r', lw=2, alpha=0.5, label='interpolated')

    ax[r].set_ylim(clip_min,clip_max)
    ax[r].set_yticklabels([])
    ax[r].set_ylabel('amplitude')
    ax[r].text(0.98,0.96,'IL {0}'.format(val),size='x-large',
    ha='right',va='top',transform=ax[r].transAxes)
ax[0].legend(loc='lower left')
ax[r].set_xlabel('XL')

f, ax = plt.subplots(1,2,figsize=(15,5))
map0 = ax[0].pcolormesh(X, Y, Aspl, cmap='cubehelix')
map1 = ax[1].pcolormesh(X, Y, Arms, cmap='cubehelix')
plt.colorbar(map0, ax=ax[0])
plt.colorbar(map1, ax=ax[1])
ax[0].set_title('amplitude')
ax[1].set_title('RMS amplitude (32ms window)')
for aa in ax.flatten():
    aa.set_xlabel('IL'), aa.set_ylabel('XL')
    CS=aa.contour(X,Y,Zf,10,colors='w')
    plt.clabel(CS, inline=1, fontsize=10, fmt='%.0f')
f.subplots_adjust(wspace=.4)

names=['raw','interpolated','rms 16ms', 'rms 48ms', 'max', 'min']
for i,aa in enumerate([Araw, Aspl, Arms16, Arms48, Amax, Amin]):
    clipmax=np.nanmax(aa)
    clipmin=np.nanmin(aa)
    print('{0}: max={1}, min={2}'.format(names[i], clipmax, clipmin))
    plt.figure(figsize=(6,6))
    plt.pcolormesh(X, Y, aa, cmap='cubehelix')
    plt.colorbar()
    plt.title(names[i])
    plt.xlabel('IL'), plt.ylabel('XL')