In [ ]:

    

import obspy
from obspy import read
import numpy as np
import matplotlib as plt
import matplotlib.pyplot as plt
from matplotlib.pyplot import gca
import matplotlib as mpl
from pylab import rcParams
from matplotlib import rc
from scipy import interpolate
#%matplotlib notebook


    

In [ ]:

    

CONV_STF = 0   # convolve spike response with source wavelet


    

In [ ]:

    

clip = 1e1      # data clip
shotno = 1      # shot number


    

In [ ]:

    

FSize = 15
font = {'color':  'black',
        'weight': 'bold',
        'size': FSize}
mpl.rc('xtick', labelsize=FSize) 
mpl.rc('ytick', labelsize=FSize) 
rcParams['figure.figsize'] = 15, 15


    

In [ ]:

    

S = 0.5


    

In [ ]:

    

str_shotno = "%0.*f" %(0,np.fix(shotno))    # shotnumber2str
filename = "../../../DENISE-Black-Edition/par/su/DENISE_MARMOUSI_p.su.shot" + str_shotno
data = read(filename, unpack_trace_headers=True)


    

In [ ]:

    

#data.filter('lowpass', freq=8.0, corners=6, zerophase=False)


    

In [ ]:

    

# number of traces, time samples and sample interval
ntr = len(data.traces)
ns = len(data.traces[0].data)                                   
dt = data.traces[0].stats.delta                                 

# x- and y-source coordinates from trace header
xsrc = data.traces[0].stats.su.trace_header.source_coordinate_x
ysrc = data.traces[0].stats.su.trace_header.source_coordinate_y

# allocate memory for traces and receiver positions
traces = np.zeros((ns, ntr))
xrec = np.zeros((1, ntr))
yrec = np.zeros((1, ntr))

i=0
for tr in data.traces:
    
    # extract traces
    traces[:, i] = tr.data[:]
    
    # x- and y-receiver coordinates from trace header
    xrec[0,i] = data.traces[i].stats.su.trace_header.group_coordinate_x
    yrec[0,i] = data.traces[i].stats.su.trace_header.group_coordinate_y
    
    i += 1

# flip traces
traces = np.flipud(traces)
    
# offset [km]    
offset = (xrec - xsrc) / 1e6


    

In [ ]:

    

# number of frequencies nf and frequency interval df
FC_low = 0.015
FC_high = 15.0

nf = 1000
df = (FC_high - FC_low) / nf

# define size arrays of real and imaginary parts
FD_real = np.zeros((nf, ntr))
FD_cmplx = np.zeros((nf, ntr))

# define binary data type
data_type = np.dtype ('float32').newbyteorder ('<')

for i in range (nf):

    str_nf = "%0.*f" %(0,np.fix(i+1))    # nf2str
    input_file = "../seis/marmousi_p_shot_1_stage_1_nfreq_" + str_nf + ".bin"
    
    fp = open(input_file)
    tmp = np.fromfile (fp, dtype=data_type)
    tmp1 = tmp.reshape(1,2*ntr)

    # collect FDFD data in arrays
    FD_real[i,:] = tmp[0:2*ntr:2]
    FD_cmplx[i,:] = tmp[1:2*ntr:2]
    
    # clean memory
    tmp = None
    tmp2 = None
    
# assemble complex data
FDFD = FD_real + 1j*FD_cmplx
tmp = np.concatenate((FDFD, np.zeros((nf, ntr)), np.zeros((nf, ntr)), np.flipud(FDFD)))
FDFD = None
FDFD = tmp

# transformation from FD to TD
FD2TD = np.fft.ifft(FDFD,axis=0)
FD2TD = np.real(FD2TD)


    

In [ ]:

    

#plt.imshow(np.real(tmp))
#plt.show()


    

In [ ]:

    

TmaxTD = dt * ns       # maximum time TD
TmaxFD = 0.25 / df     # maximum time FD
dtFD = TmaxFD / nf     # time sampling FD

# maximum time sample of FD2TD corresponding to TmaxTD
nmaxFD = np.int(TmaxTD / dtFD)

# extract FD2TD data up to TmaxTD
tmp1 = FD2TD[1:nmaxFD,:]
FD2DTD = None

# define time samples in time and frequency domain
tTD = np.arange(0.0, TmaxTD, dt)
tFD = np.arange(0.0, TmaxTD - 2*dtFD, dtFD)
offset_1 = np.arange(0, ntr, 1)

# interpolate tmp values at tTD
ftmp = interpolate.interp2d(offset_1, tFD, tmp1, kind='cubic')
tmp1 = ftmp(offset_1, tTD)
FD2TD = np.flipud(tmp1)

# clean memory
tmp = None
ftmp = None


    

In [ ]:

    

print(np.shape(tFD), np.shape(offset_1), np.shape(tmp1))


    

In [ ]:

    

if CONV_STF == 1:
    # load TD source wavelet
    filename = "../../../DENISE-Black-Edition/par/start/marmousi_II_marine_source_signal.0.su.shot" + str_shotno
    data_stf = read(filename, unpack_trace_headers=True)

    # extract TD source wavelet
    stf = np.zeros((ns, 1))

    i=0
    for tr in data_stf.traces:
        stf[:, i] = tr.data[:]
    
    # convolve source wavelet with FD2TD data
    for i in range (ntr): 
        u = np.convolve(tmp1[:,i], stf[:,0])[0:ns]
        tmp1[:,i] = u


    

In [ ]:

    

tmp2 = np.exp(-S * tTD)
tdamp = np.tile(tmp2,(ntr,1))
tdamp = np.flipud(tdamp.T)
traces = traces * tdamp

tmp2 = None
tdamp = None


    

In [ ]:

    

max_TD = np.max(traces)
max_FD = np.max(tmp1)
tmp1 = tmp1 * max_TD / max_FD
FD2TD = np.flipud(tmp1)

tmp1 = None


    

In [ ]:

    

def do_plot(n, model, cm, an, title, clip):
    
    ax=plt.subplot(1, 2, n)
    extent = [np.min(offset), np.max(offset), dt, dt*ns]

    plt.rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
    plt.rc('text', usetex=True)

    im = plt.imshow(model, cmap=cm, interpolation='bicubic', extent=extent, vmin=-clip, vmax=clip)

    a = gca()
    a.set_xticklabels(a.get_xticks(), font)
    a.set_yticklabels(a.get_yticks(), font)

    plt.title(title, fontdict=font)
    if n==1:
        plt.ylabel('Time [s]', fontdict=font)
    plt.xlabel('Offset [km]', fontdict=font)
    plt.gca().invert_yaxis()


    

In [ ]:

    

plt.close('all')
plt.figure()
do_plot(1, traces, 'gray', '(a)', r"Marmousi-2 TDFD (DENISE Black-Edition)", clip)
do_plot(2, FD2TD, 'gray', '(b)', r"Marmousi-2 FDFD (GERMAINE)", clip)
plt.tight_layout()
output_file = "Marmousi_shot_" + str_shotno + "_TDFD_FDFD.pdf"
plt.savefig(output_file, bbox_inches='tight', format='pdf')
plt.show()