#!/usr/bin/env python
# coding: utf-8

# # Reference
# 
# Mousa, W.A. and Al-Shuhail A.A., 2011, Processing of Seismic Reflection Data Using MATLAB, Morgan & Claypool

# In[1]:


get_ipython().run_line_magic('pylab', 'inline')


# In[2]:


from pkprocess import *


# In[3]:


trace=read_su('data/Book_Seismic_Data.su')


# In[4]:


get_ns(trace),get_ntr(trace),get_dt(trace)


# In[5]:


print_range(trace)


# # Header information

# In[6]:


plot(ntr_per_shot(trace),'o')
xlabel('Shot gather numbers')
ylabel('Number of traces/shot gather')
grid()


# In[7]:


plot(get_key(trace,'sx'))
xlabel('Trace number')
ylabel('Sources x-axis locations (ft)')
grid()


# In[8]:


plot(get_key(trace,'selev'))
xlabel('Trace number')
ylabel('Sources elevation (ft)')
grid()


# In[9]:


plot(get_key(trace,'gelev'))
xlabel('Trace number')
ylabel('Receivers elevation (ft)')
grid()


# In[10]:


stacking_chart(trace)


# # Displaying seismic data

# In[11]:


csg8=window(trace,'fldr',8)


# In[12]:


plot_wiggle(csg8,perc=99)


# In[13]:


plot_image(csg8,perc=99)


# In[14]:


plot_image(csg8,perc=99,cmap='bwr')


# # Trace editing

# In[15]:


plot_image(trace,figsize=[15,10],perc=99)


# In[16]:


for shot in trace_split(trace,'fldr'):
    plot_wiggle(shot)


# In[17]:


csg12=window(trace,'fldr',12)
plot_wiggle(csg12,perc=99)


# In[18]:


csg12.data[21,:]=0.
plot_wiggle(csg12,perc=99)


# In[19]:


csg16=window(trace,'fldr',16)
plot_wiggle(csg16,perc=99)


# In[20]:


csg16.data[30,:]=0.
plot_wiggle(csg16,perc=99)


# In[21]:


trace.data[ntr_per_shot(trace)*(12-1)+21,:]=0.
trace.data[ntr_per_shot(trace)*(16-1)+30,:]=0.
plot_image(trace,figsize=[15,10],perc=99)


# # Gain control

# In[22]:


csg8=window(trace,'fldr',8)
csg8_tpow=gain(csg8,tpow=2)
csg8_epow=gain(csg8,epow=2)
plot_comp((csg8_tpow,csg8_epow),plot='wiggle',perc=99)


# In[23]:


seis_env_dB(csg8,csg8_tpow)


# In[24]:


seis_env_dB(csg8,csg8_epow)


# In[25]:


trc_gained=gain(trace,tpow=2.)
plot_image(trc_gained,figsize=[15,10],perc=99)


# In[26]:


for ishot,shot in enumerate(trace_split(trc_gained,'fldr')):
    plot_wiggle(shot)
    xlabel('Trace number - shot %s'%ishot)


# # Frequency filtering

# In[27]:


specfx(csg8_tpow)


# In[28]:


specfk(csg8_tpow)


# In[29]:


csg8_bp=bpfilter(csg8_tpow,[15.,60.])
csg8_bp.print_log()
plot_comp((csg8_tpow,csg8_bp,csg8_tpow-csg8_bp),plot='wiggle',figsize=[12,10])


# In[30]:


plot_comp((csg8_tpow,csg8_bp,csg8_tpow-csg8_bp),plot='specfx',cmap='jet',perc=99,figsize=[12,10])


# In[31]:


ns=get_ns(csg8_tpow)
dt=get_dt(csg8_tpow)
tmax=(ns-1)*dt
df=1./tmax
nf=int(ns/2)
af=arange(nf)*df
t1=fft.fft(csg8_tpow.data[0])
t2=fft.fft(csg8_bp.data[0])
plot(af,abs(t1[:nf]),label='before')
plot(af,abs(t2[:nf]),label='after')
legend()
xlabel('Frequency (Hz)')
ylabel('Amplitude')
xlim([0,60])


# In[32]:


trc_bp=bpfilter(trc_gained,[15.,60.])
trc_bp.print_log()
#trc_bp.write('trace_gain_bpf.sd')


# In[33]:


trc=window(trc_bp,'fldr',[4,5,6])
plot_wiggle(trc,figsize=[10,10],perc=99)


# In[34]:


plot_image(trc_bp,figsize=[15,10],perc=99)


# In[35]:


for shot in trace_split(trc_bp,'fldr'):
    plot_wiggle(shot)


# # Autocorrelograms

# In[36]:


max_lag=0.2
auto=auto_correlation_map(trc,max_lag)
plot_wiggle(auto,figsize=[10,10],perc=95)
ylabel('Time lag (s)')


# # Spiking deconvolution

# In[37]:


mu=0.1
trc_sd=spiking_decon(trc_bp,max_lag,mu)


# In[38]:


plot_image(trc_sd,figsize=[15,10],perc=99)


# In[39]:


plot_wiggle(window(trc_sd,'fldr',[4,5,6]),figsize=[10,10],perc=99)


# # Automatic Gain Control

# In[40]:


trc_agc=gain(trc_sd,agc=True,agc_gate=0.5,norm='rms')


# In[41]:


plot_wiggle(window(trc_agc,'fldr',[4,5,6]),figsize=[10,10],perc=99)


# In[42]:


trc_agc.print_log()
#trc_agc.write('trc_gain_bpf_sd_agc.sd')


# In[43]:


for shot,sd in zip(trace_split(trace,'fldr'),trace_split(trc_agc,'fldr')):
    plot_comp((shot,sd),plot='wiggle',perc=99)


# In[44]:


plot_image(trc_agc,figsize=[15,10],perc=99)


# # CDP sorting

# In[45]:


#trc_agc=read('trc_gain_bpf_sd_agc.sd')
trc_sort=trace_sort(trc_agc,('+cdp','-offset'))
cmps=get_key(trc_agc,'cdp')
cmpu=get_key_unique(trc_agc,'cdp')
fold_num = [sum(icmp==cmps) for icmp in cmpu]

stem(cmpu,fold_num,'k-')
xlabel('CMP numbers')
ylabel('Fold')


# In[46]:


trc_sort.print_log()


# In[47]:


trc_sort.write("trc_gain_bpf_sd_agc_sort.sd")


# # Velocity analysis
# 
# 1. Open and edit `Par_velan_land.py` file.
# 2. Run `Run_velan_land.py`. (`$ python Run_velan_land.py`)
#     - Left click: Add a point
#     - Right click: Remove a point

# In[48]:


trc_nmo=read("trc_gain_bpf_sd_agc_sort_nmo.sd")
trc_nmo.print_log(nmo=True)


# In[49]:


plot_image(trc_nmo,figsize=[15,10],perc=99)


# # Stacking

# In[50]:


trc_stk=stack(trc_nmo)
trc_stk.print_log()


# In[51]:


plot_wiggle(trc_stk,perc=99,figsize=[10,10])


# In[52]:


plot_image(trc_stk,perc=99,figsize=[10,10])


# In[53]:


plot_image(trc_stk,perc=99,cmap='bwr',figsize=[10,10])


# # Surface-consistent static correction

# In[54]:


cmp_start=205
cmp_end=255
lags=20
trc_static=scr_static(trc_nmo,cmp_start,cmp_end,lags)


# In[55]:


trc_stk_static=stack(trc_static)
trc_stk_static.print_log()


# In[56]:


plot_comp((trc_stk,trc_stk_static,trc_stk_static-trc_stk),plot='wiggle',perc=99,figsize=[12,10])


# # Stolt migration

# In[57]:


# distance in feet
v=10000.
dx=110.
trc_mig=stolt_mig(trc_stk_static,v,dx)


# In[58]:


plot_wiggle(trc_mig,perc=99,figsize=[10,10])


# In[59]:


plot_wiggle(trc_mig-trc_stk_static,perc=99,figsize=[10,10])


# In[60]:


trc_mig.print_log()


# In[ ]: