Python function to plot inline, crossline or horizontal slice from 3D seismic volume. Also included example showing how to load a 3D cube and display it.

example_plot3dseis.py

import numpy as np
import matplotlib.pyplot as plt
import time
from obspy.io.segy.segy import _read_segy

# sample 3d cube from GOM, download here:
# https://walrus.wr.usgs.gov/namss/survey/b-05-88-la/
filename='G3D1304_005_L88_056_LN_1664-3428.sgy'

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
t0=time.time()
segy = _read_segy(filename)
print('--> data read in {:.1f} sec'.format(time.time()-t0))
ntraces=len(segy.traces)
sr=segy.binary_file_header.sample_interval_in_microseconds/1e3 # sr is 30 ft
ns=segy.binary_file_header.number_of_samples_per_data_trace
size_kb=segy.traces[0].data.nbytes/1024.*ntraces
z=np.arange(0,ns*sr,sr)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
inl = np.arange(1664,3429,4)
crl = np.arange(29280,32871,10)
seis = np.vstack([xx.data for xx in segy.traces]).T
seis = seis.reshape(ns,inl.size,crl.size)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
plot3dseis(seis,z,inl,crl,d1_sel=None,d2_sel=2736,d3_sel=None,id1='depth [ft]',id2='IL',id3='XL')
plt.ylim(15000,0),plt.clim(-150,150)

plot3dseis.py

def plot3dseis(seis,d1,d2,d3,d1_sel=None,d2_sel=None,d3_sel=None,id1='Z',id2='IL',id3='XL',cm='RdBu'):
    '''
    plot3dseis (C) aadm 2016
    Plot inline, crossline or horizontal (time/depth) slice from 3D seismic volume.

    INPUT
    seis:        3D seismic cube numpy array, shape (a,b,c);
                 usually this will be (twt.size, inl.size, crl.size)
    d1, d2, d3:  numpy array containing the range of first, second and third dimension in the cube;
                 usually d1=twt, d2=inl, d3=crl are defined beforehand,
                 e.g. inl=np.arange(1000,2000,4)
    d1_sel, ...: to select one horizontal slice, one inline or on crossline
                 e.g., d2_sel=1500 to plot IL 1500 (if d2=inl); only set one, the others must be None
    id1, ...:    names for first, second, third dimensions
    cmap:        colormap to use (choose one of the matplotlib colormaps)
    '''

    if (d1_sel is None) & (d2_sel is None) & (d3_sel is None):
        d2_sel = d2[int(d2.size/2)]
    if d1_sel is not None:
        ssplot=seis[d1.tolist().index(d1_sel),:,:]
        nome1,nome2,nomex,nomey = id1,str(d1_sel),id2,id3
        xleft,xrite = d2.min(),d2.max()
        ytop,ybot   = d3.min(),d3.max()
    if d2_sel is not None:
        ssplot=seis[:,d2.tolist().index(d2_sel),:]
        nome1,nome2,nomex,nomey = id2,str(d2_sel),id3,id1
        xleft,xrite = d3.min(),d3.max()
        ytop,ybot   = d1.min(),d1.max()
    if d3_sel is not None:
        ssplot=seis[:,:,d3.tolist().index(d3_sel)]
        nome1,nome2,nomex,nomey = id3,str(d3_sel),id2,id1
        xleft,xrite = d2.min(),d2.max()
        ytop,ybot   = d1.min(),d1.max()

    print('input 3D dimensions: {}'.format(seis.shape))
    print('{0}: {1}-{2}'.format(id1, d1.min(), d1.max()))
    print('{0}: {1}-{2}'.format(id2, d2.min(), d2.max()))
    print('{0}: {1}-{2}'.format(id3, d3.min(), d3.max()))
    print('dataset to plot dimensions: {}'.format(ssplot.shape))

    plt.figure(figsize=(10,8))
    plt.imshow(ssplot,extent=[xleft,xrite,ybot,ytop],aspect='auto',cmap=cm)
    plt.colorbar(shrink=0.75), plt.grid()
    plt.tick_params(axis='both', which='major', labelsize=8)
    plt.title(nome1+' '+nome2)
    plt.xlabel(nomex), plt.ylabel(nomey)