Plot Marmousi RTM Results

Daniel Köhn

Kiel, 02/04/2016

Import necessary packages

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

FD grid dimensions

In [2]:
DH = 20.0;
NX = 500;
NY = 174;

Define fonts

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

Read S-wave velocity model and RTM result

In [4]:
f = open ('../start/marmousi_II_smooth2.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 [5]:
f1 = open ('../jacobian/jacobian_Test_P_image')
RTM = np.fromfile (f1, dtype=data_type)
RTM = RTM.reshape(NX,NY)
RTM = np.transpose(RTM)
RTM = np.flipud(RTM)
RTM = scipy.ndimage.filters.laplace(RTM) # suppress low-wavenumber artifacts in image

Define Axis

In [6]:
x = np.arange(0.0, DH*NX/1000.0, DH)
y = np.arange(0.0, DH*NY/1000.0, DH)

Plot $\alpha$-Blending of Vp FWI result (Jet) and Laplace-filtered RTM result (Gray)

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

fig = plt.figure(frameon=False)

plt.rc('text', usetex=True)

im1 = plt.imshow(vp,, interpolation='nearest', extent=extent)

im2 = plt.imshow(RTM,, alpha=.70, interpolation='bicubic',
                 extent=extent, vmin=cmin, vmax=cmax)

#im2 = plt.imshow(RTM,, interpolation='bicubic',
#                 extent=extent, vmin=cmin, vmax=cmax)

a = gca()
a.set_xticklabels(a.get_xticks(), font)
a.set_yticklabels(a.get_yticks(), font)
plt.title('Marmousi-RTM (smooth2 model)', fontdict=font)
plt.ylabel('Depth [km]', fontdict=font)
plt.xlabel('Distance [km]', fontdict=font)
#cbar.set_label('Vp[m/s]', fontdict=font, labelpad=1)
#plt.savefig('Marmousi_RTM.pdf', format='pdf')
#plt.savefig('Marmousi_RTM.png', format='png')

In [ ]: