This is a short tutorial on how to load a matrix of numbers from a textfile into a 2-dimensional array, and then reshape that array into 3-dimensional array and make some annotated plots along each of the dimensions of a seismic volume.
This notebook is associated with the blog post, Slicing seismic arrays.
First we load the NumPy
package for number crunching and array manipulation. Load the pyplot
module from matplotlib
for plotting. It is conventional to give these packages the nicknames np
and plt
respectively.
import numpy as np
import matplotlib.pyplot as plt
%pylab inline
Populating the interactive namespace from numpy and matplotlib
The seismic data in this example has been made ultra-simple for illustration purposes. It a simple text file comprised of columns of numbers, with no header information at the start of the file. We can load this text file into our workspace by using the numpy loadtxt
command.
data = np.loadtxt('data/slicing_seismic.txt')
The np.loadtxt
function returns an array. Arrays have an intrinsic attribute called shape
that we can call to inspect the size of the array,
print data.shape
(10000, 500)
data.shape()
returns a tuple of the number of elements within each dimension of the array. This says that data
is a 2-dimensional array, with 100000 elements along the first dimension (traces), and 500 elements along the second dimensions (samples per trace).
Since data.shape()
returns a tuple, we can assign two new variables for describing the size of each of the dimensions in this textfile. The size of the first dimension is equal to the number of inlines, nIL
multiplied by the the number of crosslines, nXL
,
nILnXL, nt = data.shape
In our case, and by the way I've designed this example, the number of inlines equals the number of crosslines (100).
nIL = np.sqrt(nILnXL)
nXL = np.sqrt(nILnXL)
Now use the the variables, nIL
, nXL
, and nt
to np.reshape()
our 2D array into its proper 3D shape:
new_data = np.reshape(data, (nIL, nXL, nt))
print new_data.shape
(100, 100, 500)
Index into the time slice (0 to 499), an inline (0 to 99), and a crossline (0 to 99)
z = 122
inline = 30
xline = 60
Create a figure with three panels: timeslice, inline and crossline
f = plt.figure(figsize=(12,4), facecolor='white')
# Plot timeslice
ax = f.add_axes([0.1, 0.05, 0.35, 0.85])
ax.imshow(np.transpose(new_data[:,::-1,z]),
cmap = 'Greys',
vmin = -4000, vmax = 4000,
origin='lower')
ax.set_title('time slice at %d (ms)' % z, fontsize = 12)
# Plot inline
ax2 = f.add_axes([0.525, 0.05, 0.2, 0.85])
ax2.imshow(np.transpose(new_data[:,inline,:200]),
cmap = 'Greys',
vmin = -4000, vmax = 4000,
origin='upper')
ax2.set_title('inline %d' % inline, fontsize = 12)
# Plot xline
ax3 = f.add_axes([0.75, 0.05, 0.2, 0.85])
ax3.imshow(np.transpose(new_data[xline,:,:200]),
cmap = 'Greys',
vmin = -4000, vmax = 4000,
origin='upper')
ax3.set_title('crossline %d' % xline, fontsize = 12)
plt.show()
Let's add some code to make a number of plotting annotations and decorations. Admittedly it is a rather verbose.
f = plt.figure(figsize=(12,4), facecolor='white')
ax = f.add_axes([0.1, 0.05, 0.35, 0.85])
ax.imshow(np.transpose(new_data[:,::-1,z]),
cmap = 'Greys',
vmin = -4000, vmax = 4000,
origin='lower')
ax.set_title('time slice at %d (ms)' % z, fontsize = 12)
# Plot inline
ax2 = f.add_axes([0.525, 0.05, 0.2, 0.85])
ax2.imshow(np.transpose(new_data[:,inline,:200]),
cmap = 'Greys',
vmin = -4000, vmax = 4000,
origin='upper')
ax2.set_title('inline %d' % inline, fontsize = 12)
# Plot xline
ax3 = f.add_axes([0.75, 0.05, 0.2, 0.85])
ax3.imshow(np.transpose(new_data[xline,:,:200]),
cmap = 'Greys',
vmin = -4000, vmax = 4000,
origin='upper')
ax3.set_title('crossline %d' % xline, fontsize = 12)
ax3.hlines(z, 0,100, 'r', lw=3, alpha = 0.5)
ax2.hlines(z, 0,100, 'r', lw=3, alpha = 0.5)
ax.hlines(inline, 0,100,'g', lw=3, alpha = 0.5)
ax.vlines(xline, 0,100, 'b', lw=3, alpha = 0.5)
ax2.vlines(xline,0,200,'b', lw=3, alpha = 0.5)
ax3.vlines(inline,0,200,'g', lw=3, alpha = 0.5)
# Color and line weight of time slice
ax.spines['bottom'].set_color('red')
ax.spines['bottom'].set_lw('3')
ax.spines['top'].set_color('red')
ax.spines['top'].set_lw('3')
ax.spines['right'].set_color('red')
ax.spines['right'].set_lw('3')
ax.spines['left'].set_color('red')
ax.spines['left'].set_lw('3')
# Color and line weight of inline
ax2.spines['bottom'].set_color('green')
ax2.spines['bottom'].set_lw('3')
ax2.spines['top'].set_color('green')
ax2.spines['top'].set_lw('3')
ax2.spines['right'].set_color('green')
ax2.spines['right'].set_lw('3')
ax2.spines['left'].set_color('green')
ax2.spines['left'].set_lw('3')
# Color and line weight of xline
ax3.spines['bottom'].set_color('blue')
ax3.spines['bottom'].set_lw('3')
ax3.spines['top'].set_color('blue')
ax3.spines['top'].set_lw('3')
ax3.spines['right'].set_color('blue')
ax3.spines['right'].set_lw('3')
ax3.spines['left'].set_color('blue')
ax3.spines['left'].set_lw('3')
# Set the colors on the axis labels to match
ax.tick_params(axis='x', colors='blue')
ax.tick_params(axis='y', colors='green')
ax2.tick_params(axis='x', colors='blue')
ax3.tick_params(axis='x', colors='green')
ax.set_ylabel('inline', fontsize=10)
ax.set_xlabel('crossline', fontsize=10)
ax2.set_ylabel('twt [ms]', fontsize=10)
ax2.set_xlabel('crossline', fontsize=10)
ax3.set_ylabel('twt [ms]', fontsize=10)
ax3.set_xlabel('inline', fontsize=10)
ax3.axes.get_yaxis().set_visible(False)
plt.show()