I want to display amplitude together with semblance. We'll try combining them into a single display, and making a 2D colourmap to go with the display.
This notebook goes with the blog post, Corendering attributes and 2D colourmaps.
First, the usual preliminaries.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
% matplotlib inline
Let's load two time slices from disk. The first one is a time slice of the seismic amplitudes, the second is a time slice from a continuity attribute. Continuity is also sometimes referred to as similarity, variance, or coherence. These are edge-detection type attributes that are useful for looking for faults, fractures, and other discontinuities in the data.
A = np.load('data/amp_slice.npy')
C = np.load('data/cont.npy')
The arrays should be exactly the same size:
A.shape == C.shape
True
This is just a function to scale the amplitudes between 0-1:
vmin = -0.75*A.ptp() # make axis 50% larger than max/min of trace
vmax = 0.75*A.ptp()
Anorm = 0.5 + (A / vmax)
Let's use the 'RdBu' colormap for seismic.
img_array = plt.get_cmap('RdBu')(Anorm)
Let's plot the data. In the first panel, we'll plot the amplitudes. In the second panel, we'll plot the continuity attribute, in the third panel we'll set the continuity image to 75% opacity, and lay it ontop of the amplitudes to see both data in the same space.
Inline and crossline spacing for these data is 25 m.
w = 0.025 # km
import matplotlib.pyplot as plt
from matplotlib.path import Path
import matplotlib.patches as patches
def add_scale(ax, xo=10., yo=-10., dx=40., dy=2.5, spacing=25., units='km'):
distance = int(spacing * dx / 1000.0)
verts = [
(xo, yo), # left, bottom
(xo, yo+dy), # left, top
(xo+dx, yo+dy), # right, top
(xo+dx, yo), # right, bottom
(0., 0.), # ignored
]
codes = [Path.MOVETO,
Path.LINETO,
Path.LINETO,
Path.LINETO,
Path.CLOSEPOLY,
]
path = Path(verts, codes)
patch = patches.PathPatch(path, facecolor='w', lw=1, alpha=0.5)
ptch = ax.add_patch(patch)
patch.set_clip_on(False) # to draw scale outside axes boundary
ax.text(xo+dx, yo, ' ' + str(distance) + ' ' + units, horizontalalignment='left', verticalalignment='center')
return
First, we can try just making one attribute slightly transparent, and overlaying on the other. We can do this by just using imshow
twice, once with alpha
(opacity) set to something less than 1.
fig = plt.figure(figsize=(15,5))
ax1 = fig.add_subplot(131)
ax2 = fig.add_subplot(132)
ax3 = fig.add_subplot(133)
ax1.set_title('amplitude')
ax1.imshow(img_array,origin='bottom')
ax1.set_yticklabels([])
ax1.set_xticklabels([])
add_scale(ax1)
ax2.set_title('continuity')
ax2.imshow(C, cmap='Greys_r', origin='bottom')
ax2.set_xticklabels([])
ax2.set_yticklabels([])
ax3.set_title('amplitude with continuity')
ax3.imshow(img_array, cmap='coolwarm_r', origin='bottom')
ax3.imshow(C, cmap='Greys_r', alpha=0.75, origin='bottom')
ax3.set_xticklabels([])
ax3.set_yticklabels([])
plt.show()
#fig.savefig('images/Amp_Cont_3panel.png')
But this presentation isn't ideal. It looks washed out. Let's try modulating the lightness in the seismic amplitude image by using the continuity.
In order to change the lightness of the colourmap, we need to convert the red, green, blue (RGB) triplets of the image into hue, saturation, and value (HSV) triplets first. Then we can multiply the "value" channel by the continuity (which already fortuitously has a range between 0 and 1.0).
In order to convert an image of RGB triplets into HSV, and vice versa, we need to import two new functions from the matplotlib.colors module,
from matplotlib.colors import rgb_to_hsv, hsv_to_rgb
Modulate lightness with continuity...
hsv = rgb_to_hsv(img_array[:,:,:3])
hsv[:, :, 2] *= C
corendered = hsv_to_rgb(hsv)
Let's remake our plot
fig = plt.figure(figsize=(15,5))
ax1 = fig.add_subplot(131)
ax2 = fig.add_subplot(132)
ax3 = fig.add_subplot(133)
ax1.set_title('amplitude')
ax1.imshow(img_array, origin='bottom')
ax1.set_xticklabels([])
ax1.set_yticklabels([])
add_scale(ax1)
ax2.set_title('continuity')
ax2.imshow(C, cmap='Greys_r', origin='bottom')
ax2.set_xticklabels([])
ax2.set_yticklabels([])
ax3.set_title('amplitude and continuity')
ax3.imshow(corendered, origin='bottom')
ax3.set_xticklabels([])
ax3.set_yticklabels([])
plt.show()
#fig.savefig('images/Amp_Cont_3panel_better.png')
Not bad!
In order to make a 2D colourmap to place as our legend, we need to import the makeMappingArray function from the matplotlib.colors module. We'll also import all the colormaps (called cm) that matplotlib has to offer, so you can try this out with any one of the built in colour maps.
For a list of all the default colormaps, see http://wiki.scipy.org/Cookbook/Matplotlib/Show_colormaps (Remember, to reverse the direction of any of these colourmaps, simply append an "_r" to the name when you call it.
from matplotlib import cm
from matplotlib.colors import makeMappingArray
Choose a colour map from Matplotlib's list of standard colourmaps.
seisbar = cm.get_cmap('RdBu')
Let's make an colourmap array that has 32 colours in it. This number can be any number you like. We probably wouldn't see much difference for any more than 64 colours, if the colour map is smoothly varying.
ncolours = 32
seis_array = makeMappingArray(ncolours, seisbar)
seis_array
is an RGBA representation of our colourmap. But it has the shape,
np.shape(seis_array)
(32, 4)
So we need to change the shape of the array, by adding an dimension to it. We can also drop the A (alpha) channel, because we aren't going to be using that,
color_arr = seis_array[np.newaxis, :]
color_arr = color_arr[:,:,:-1]
Now, we can build out a square image for plotting (32 x 32 x 3) using NumPy's tile function,
colour_roll = np.rollaxis(color_arr, 1)
seis_rgb_mtx = np.tile(colour_roll, (32,1))
np.shape(seis_rgb_mtx)
plt.imshow(seis_rgb_mtx)
plt.show()
Since matplotlib
uses RGB to diplay graphics, passing the HSV representation to the plot looks incorrect. But, now that we have an HSV representation of our colourmap, we are ready to modulate the lightness values across the image
seis_hsv = rgb_to_hsv(seis_rgb_mtx)
Create a ramp of lightness that goes from 0 to 1 along the horizontal axis of the colormap, then modulate "values" in last slice.
hues, lightness = np.mgrid[0:1:32j, 0:1.0:32j]
seis_hsv[:,:,2] *= lightness
plt.imshow(seis_hsv)
plt.show()
But in order to render out modified image we need to conver the triplets back into RGB format,
RGB = hsv_to_rgb(seis_hsv)
Make a quick function to label the plots...
def plot_annot(ax, text1='attribute 1', text2='attribute 2'):
"""A helper function to label the axes, to save typing."""
ax.set_ylabel(text1, rotation=0, ha='left', va='center')
ax.yaxis.set_label_coords(1.025, 0.5)
ax.yaxis.tick_right()
ax.set_xlabel(text2)
ax.set_xticklabels(['', 'low', '', '', 'high', ''])
ax.set_yticklabels(['', 'neg.', '', '', 'pos.', ''])
return
fig = plt.figure(figsize = (4,5))
ax = fig.add_subplot(111,axisbg='k')
ax.imshow(RGB, origin="lower", extent=[0, 1, 0, 1])
plot_annot(ax, 'amplitude', 'continuity')
plt.show()
#fig.savefig('Amp_Cont_2D_colourmap.png')
fig = plt.figure(figsize=(15,5), facecolor='white')
ax = fig.add_axes([0.1,0.1,0.5,0.8])
ax.set_title('amplitude and continuity')
ax.imshow(corendered, origin='bottom')
ax.set_xticklabels([])
ax.set_yticklabels([])
add_scale(ax)
cax = fig.add_axes([0.2,0.1,0.9,0.8], axisbg='k')
cax.imshow(RGB, origin="lower", extent=[0, 1, 0, 1])
plot_annot(cax, 'amplitude', 'continuity')
plt.show()
#fig.savefig('images/Data_Amp_Cont_2D_colourmap.png')
Now every value in the image space represents an specific value in the x-y plane of the 2-D colourmap legend.
Instead of using a linear opacity function, we might want to try a nonlinear function — that way we can emphasize the biggest discontinuities.
We'll define a sigmoid curve:
def opaque_transform(x, mu, s):
"""
Using a cumulative distribution of logistic function:
http://en.wikipedia.org/wiki/Logistic_distribution
x: the x values
mu: std deviation
s: horizontal shift
"""
y = 0.5 + 0.5*(np.tanh((x-mu)/s))
return y
And we need a way to plot it:
def opacity_plot(ax, data):
"""
ax: axes to plot in
data: data to make histogram
"""
# Histogram
hy, hx = np.histogram(np.ravel(data), 100)
denom = np.amax(np.percentile(hy,98)) # scaling correction for plot
hy = hy / (1.1*float(denom))
ax.fill_between(hx[1:-1], hy[1:], 0, color='k', alpha=0.5)
# Opacity curve
x = np.linspace(0,1.0,100)
y = opaque_transform(x, mu, s)
mid = int(len(x)/2.0)
ax.plot(x,y, '-', lw=3, alpha = 0.5)
xbar, ybar = mu, opaque_transform(mu, mu, s)
ax.scatter(xbar, ybar, s=50, alpha=1.0)
ax.text(xbar-0.025, ybar+0.05, s='mu: '+'%.2f' %mu, alpha=1.0, ha='right', va='center',
fontsize=10)
ax.text(xbar-0.025, ybar-0.05, s='s: '+'%.2f' %s, alpha=1.0, ha='right', va='center',
fontsize=10)
ax.set_ylabel('transparency')
ax.yaxis.set_label_position("right")
ax.yaxis.set_ticks_position("right")
ax.set_xlabel('continuity')
ax.set_xlim([0,1.0])
ax.set_ylim([0,1.0])
ax.grid()
return ax
Here's what that does... the left-hand plot is just a map (with imshow
), the right-hand plot is the histogram-plus-sigmoid-function curve we jsut defined.
# Scale the opacity
mu = 0.4 # lateral shift
s = 0.2 # steepness of sigmoid (small number = fast)
Cscaled = opaque_transform(C, mu, s)
fig = plt.figure(figsize=(14,5))
# Left-hand plot
ax = fig.add_subplot(121)
plt.imshow(Cscaled, cmap='Greys_r', origin='bottom', vmin=0.0, vmax=1.0)
plt.colorbar(shrink=0.5)
plt.title('continuity')
# Right-hand plot
ax2 = fig.add_subplot(122)
ax2 = opacity_plot(ax2, C)
plt.show()
I'll also make a way to display the 2D colourbar, as before.
hues, sats = np.mgrid[0:1:32j, 0:1:32j]
Cscbar = opaque_transform(sats, mu, s)
opacity = np.expand_dims(Cscbar, -1)
RGBAmap = np.dstack((seis_rgb_mtx, opacity))
fig = plt.figure()
ax = fig.add_axes([0.1,0.1,0.5,0.8], axisbg='k')
ax.imshow(RGBAmap)
plt.show()
Now we can apply our new opacity map to the actual data, creating a new version of the amplitude map.
img_array = plt.get_cmap('RdBu')(Anorm)
img_array = img_array[:,:,:3]
opacity2 = np.expand_dims(Cscaled, -1)
RGBA = np.dstack((img_array, opacity2))
And plot this alongside the 2D colourmap and the histogram with the opacity curve.
fig = plt.figure(figsize=(12,6), facecolor='w')
# Main map
ax = fig.add_axes([0.1,0.1,0.5,0.8], axisbg='k')
ax.set_title('amplitude with continuity')
ax.imshow(RGBA, origin='bottom')
ax.set_xticklabels([])
ax.set_yticklabels([])
add_scale(ax)
# Colourbar
cax = fig.add_axes([0.5,0.55,0.35,0.35],axisbg='k')
cax.imshow(RGBAmap, origin="lower", extent=[0, 1, 0, 1])
plot_annot(cax, 'amplitude', 'continuity')
# Opeacity curve
grax = fig.add_axes([0.588,0.1,0.174,0.35], axisbg='w') # this is a bit clunky, adjusting manually to line up
opacity_plot(grax, C)
plt.show()
#fig.savefig('images/Data_Amp_Cont_Sculplted.png')
Nice!
© 2015 Agile Geoscience — CC-BY — Have fun!