Visualizing 3D Seismic Volumes Made Easy with Python and Mayavi

Visualization of migrated, post-stack seismic volumes is a very crucial component of interpretation workflows, be it to pick salt domes, interpret horizons, identify fault planes, or classify rock facies.

While there are already professional tools like OpendTect and Petrel available to experienced seismic interpreters, setting up a project from scratch in either of these software tools can be overwhelming for beginners and non-specialists, considering all of the various parameters that need to be set before one can visualize the seismic volume of interest. This is not to mention the inconvenience of working with different data types for processing and loading data for visualization, especially in the context of Data Science projects with seismic volumes – repeatedly converting seismic data back and forth between different data types at different stages of processing can be very time consuming and inconvenient for users.

I recently came across a very handy, open-source visualization tool, Mayavi, that allows seamless integration of visualization of 3D data like seismic volumes into data science workflows with Python. It is made all the better owing to its very interactive capabilities, where one may drag inline/crossline/depth slices throughout the volume at different angles to examine the features of interest. An example screenshot of such a 3-D visualization is shown below.

This post will walk you through some basic code to get you up and running with using Mayavi for visualizing seismic volumes along with their various interpreted features by seismic interpreters.

Installation

Detailed instructions for downloading Mayavi and setting it up for use with Python may be found directly on its homepage. Refer here for more details.

Importing Mayavi into Python Environment

For our purposes, we will be concerned with using the mlab API for using Mayavi as a 3-D plotting tool in python scripts. Assuming everything went smoothly in the last step, use the following code snippet to import mlab:

import numpy as np         # import numpy for data loading
from mayavi import mlab    # import mlab

Loading Seismic Data into the Python Session

Seismic datasets usually come in the popular .segy (pronounced: seg-y) format. In addition to storing amplitude data, segy files contain a variety of header information, like the crossline and inline numbers associated with each trace, time sample information etc. However, popular machine learning and data science tools usually ask for data to be in the form of numpy arrays. If your data happens to be in segy format, you may use the popular python package, Segyio, to convert it to a numpy array. See here, for more details, instructions, and examples on using Segyio, refer to their GitHub page.

For the rest of this post, we assume you are working with the Netherlands Offshore F3 Block’s seismic data and its interpreted labels, as developed by the OLIVES lab at Georgia Tech. For more information on this work, refer to the Github page here. Notice that converting the data to numpy format will make it lose tertiary information contained in the headers of the original segy file. But this is OK, since we are usually interested just in the amplitudes and labels for the specific purpose of visualization. Use the following code snippet to load both the raw seismic amplitudes and the associated rock facies labels as 3-D arrays.

seismic_path = 'data/train/train_seismic.npy'  
label_path = 'data/train/train_labels.npy'     

seismic = np.load(seismic_path)     # load 3-D seismic
labels = np.load(label_path)        # load 3-D labels

Setting up a Figure Object

Creating a visualization instance with Mayavi consists of two components: setting up a figure environment and secondly, populating this figure with views of the data one is interested in. This concept is somewhat similar to how one performs visualization with the popular Matplotlib library in Python. Set up the figure object using the code below:

fig = mlab.figure(figure='seismic', bgcolor=(1, 1, 1), fgcolor=(0, 0, 0))

Viewing Inline/Crossline/Depth Slices in the Seismic Volume

This is where we add views of data to the figure object we just created. To view orthogonal slices in 3-D data arrays, the function in mlab API we are interested in is _volumeslice. It takes as input the data array, the index of the slice one is interested in, the orientation of the slice, and the figure object to populate with this information. You may add as many slices in either of the three orientations you want. The slice index specifies the slice you will be looking at the first time the figure pops into view, but being interactive, you may then drag it to the position you want. See the code snippet below for visualizing slices along each of the crossline, inline, and depth orientations in the seismic volume we loaded above.

scalars = seismic   # specifying the data array

mlab.volume_slice(scalars, slice_index=0, plane_orientation='x_axes', figure=fig)   # crossline slice

mlab.volume_slice(scalars, slice_index=0,  plane_orientation='y_axes', figure=fig)   # inline slice

mlab.volume_slice(scalars, slice_index=0, plane_orientation='z_axes', figure=fig)   # depth slice

mlab.show()

This produces the visualization shown below:

You may want to explore different colormaps to suit best your application. One may also add information to the axes in this figure by a slightly modified version of the code snippet above:

scalars = seismic   # specifying the data array

mlab.volume_slice(scalars, slice_index=0, plane_orientation='x_axes', figure=fig)   # crossline slice

mlab.volume_slice(scalars, slice_index=0,  plane_orientation='y_axes', figure=fig)   # inline slice

mlab.volume_slice(scalars, slice_index=0, plane_orientation='z_axes', figure=fig)   # depth slice

mlab.axes(xlabel='Inline', ylabel='Crossline', zlabel='Depth', nb_labels=10)                             # Add axes labels 

mlab.show()

This produces the figure below:

Rendering Multiple Data Arrays on the Same Figure

Finally, one may render multiple data arrays onto the same figure object for an even more enriching experience with visualization. The interpretations of the seismic volume we used label all voxels in salt domes as the integer ‘3’. This is done by invoking the mlab API’s function, contour3d, that plots iso-surfaces in a 3-D volume given a list of integer(s) associated with the iso-surface(s). The following snippet of code overlays the interpreted salt dome in the migrated dataset on to the raw seismic amplitudes, producing the figure shown underneath:

scalars = seismic   # specifying the data array

mlab.volume_slice(scalars, slice_index=0, plane_orientation='x_axes', figure=fig)   # crossline slice

mlab.volume_slice(scalars, slice_index=0,  plane_orientation='y_axes', figure=fig)   # inline slice

mlab.volume_slice(scalars, slice_index=0, plane_orientation='z_axes', figure=fig)   # depth slice

mlab.contour3d(labels, contours=[3], figure=fig)  # plot isosurface   

mlab.show()

Summary

Visualization is a crucial tool in any interpreter’s toolbox. Conventional seismic visualization packages offer several inconveniences when it comes to visualizing migrated seismic volumes, especially for data science workflows. This is where the popular visualization package, Mayavi comes in, providing several handy features for quickly and efficiently setting up visualizations of migrated seismic volumes and their interpretations. We walked you through using some of these features in a Python environment. We hope this post proves beneficial to interpreters and non-experts in their interpretation workflows involving Python.