Plot SH 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 pickle
import scipy.ndimage.filters
from scipy.ndimage.filters import gaussian_filter
from matplotlib.ticker import ScalarFormatter, FormatStrFormatter
Activate different Post-Processing Options
In [15]:
GAUSSIAN=1;
Import Colormap
In [16]:
fp = open('cmap_cm.pkl', 'rb')
my_cmap_cm = pickle.load(fp)
fp.close()
FD grid dimensions
In [17]:
DH = 0.05;
NX = 512;
NY = 256;
vpmin=800.0
Define fonts
In [18]:
FSize = 20
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 [19]:
f = open ('../inv_result_new_2016_profile_1/Single_vel_corr/PCG_04_12_2016_vs_rho_Qs_15/modelTest_vs_stage_4.bin')
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 [20]:
f1 = open ('SH_image_profile_1/07_12_2016_FWI_smooth_70_100_Hz/RTM_profile_1_FWI_smooth.bin')
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
Gaussian filter
In [21]:
if(GAUSSIAN==1):
RTM = gaussian_filter(RTM, sigma=[1,10])
Define Axis
In [22]:
x = np.arange(0.0, DH*NX, DH)
y = np.arange(0.0, DH*NY, DH)
Scale RTM result with depth
In [23]:
RTM_scale = np.zeros((NX,NY))
RTM_scale += np.flipud(y)**0
RTM*=RTM_scale.transpose()
Plot SubPlots
In [24]:
#plt.close('all')
#plt.figure()
#do_plot(1, 'inferno', "Inferno")
#do_plot(2, 'viridis', "Viridis")
#do_plot(3, my_cmap_cm, "My cool Colormap")
#plt.savefig('test.png', format='png', dpi=100)
#plt.savefig('test.pdf', format='pdf')
#plt.show()
In [25]:
def do_plot(n, an, title):
ax=plt.subplot(2, 1, n)
extent = [DH,NX*DH,DH,NY*DH]
cmax=1e-4
cmin=-cmax
plt.rc('font',**{'family':'sans-serif','sans-serif':['Helvetica']})
plt.rc('text', usetex=True)
im1 = plt.imshow(RTM, cmap=plt.cm.gray, interpolation='none',
extent=extent, vmin=cmin, vmax=cmax)
plt.hold(True)
if(n==2):
im2 = plt.imshow(vp, cmap=plt.cm.jet, alpha=.7, interpolation='bicubic',
extent=extent)
a = gca()
a.set_xticklabels(a.get_xticks(), font)
a.set_yticklabels(a.get_yticks(), font)
a.xaxis.set_major_formatter(FormatStrFormatter('%0.0f'))
a.yaxis.set_major_formatter(FormatStrFormatter('%0.0f'))
#plt.axis('scaled')
plt.title(title, fontdict=font)
plt.ylabel('Depth [m]', fontdict=font)
if(n==2):
plt.xlabel('Distance [m]', fontdict=font)
plt.gca().invert_yaxis()
#cbar=plt.colorbar()
#cbar.set_label('Vp[m/s]', fontdict=font, labelpad=1)
In [26]:
plt.close('all')
plt.figure()
do_plot(1, '(a)', r"SH-RTM (smooth FWI result)")
do_plot(2, '(b)', r" ")
#plt.savefig('test.png', format='png', dpi=100)
plt.savefig('test.pdf', bbox_inches='tight', format='pdf')
plt.tight_layout()
plt.show()
In [ ]: