In [1]:

    

from __future__ import division
import matplotlib.pyplot as plt
import matplotlib as mpl
import numpy as np
from matplotlib.colors import LightSource, Normalize
from matplotlib.pyplot import gca
from pylab import rcParams
from matplotlib import rc
import scipy.ndimage.filters
import pickle


    

In [2]:

    

DH = 20.0;
NX = 500;
NY = 174;
N = NX*NY;


    

In [9]:

    

clip = 6e0         # wavefield clipping
NSNAP1 = 1         # first snapshot
NSNAP2 = 50        # last snapshot
DSNAP = 1          # snapshot increment
TSNAP1 = 0.002     # time of first snapshot
TSNAP2 = 3.0       # time of last snapshot
TSNAPINC = 0.06    # time increment between snapshots


    

In [4]:

    

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


    

In [5]:

    

f = open ('../start/marmousi_II_marine.vp')
data_type = np.dtype ('float32').newbyteorder ('<')
vp = np.fromfile (f, dtype=data_type)
vp = vp.reshape(NX,NY)
vp = np.transpose(vp)
vp = np.flipud(vp)


    

In [6]:

    

def load_snap(offset):
    f1 = open ('../snap/waveform_forward.bin.p')
    data_type1 = np.dtype ('float32').newbyteorder ('<')
    offset_snap = 4*NX*NY*(offset-1);
    snap = np.memmap(f1, dtype=data_type1, mode='c', shape=(NX,NY), offset=offset_snap)
    snap = np.transpose(snap)
    snap = np.flipud(snap)
    return snap


    

In [7]:

    

x = np.arange(0.0, DH*NX/1000.0, DH)
y = np.arange(0.0, DH*NY/1000.0, DH)


    

In [10]:

    

extent = [0.0,NX*DH/1000.0,0.0,NY*DH/1000.0]

fig = plt.figure(frameon=True)

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

i = NSNAP1
while i <= NSNAP2:
    
    offset = i;
    time = TSNAP1 + (i-1) * TSNAPINC

    snap = load_snap(offset);

    im1 = plt.imshow(vp, cmap=plt.cm.gray, interpolation='nearest', extent=extent)
    plt.hold(True)

    im2 = plt.imshow(snap, cmap=plt.cm.seismic, alpha=.75, 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.axis('scaled')
    time_str = str(time)
    name_title = "Pressure wavefield (time = " + time_str[0:5] + " s)"
    plt.title(name_title, fontdict=font)
    plt.ylabel('Depth [km]', fontdict=font)
    plt.xlabel('Distance [km]', fontdict=font)
    plt.gca().invert_yaxis()
    #cbar=plt.colorbar()
    #cbar.set_label('Vp[m/s]', fontdict=font, labelpad=1)
    plt.tight_layout()
    #plt.savefig('Marmousi_RTM.pdf', format='pdf')
    name_snap = "Marmousi_snap_" + "%0.*f" %(0,np.fix(offset)) + ".pdf"
    plt.savefig(name_snap, format='pdf')
    plt.hold(False)
    #plt.show()
    
    i = i + DSNAP;


    

In [ ]: