Skip to content

Instantly share code, notes, and snippets.

@aadm
Last active May 9, 2024 14:10
Show Gist options
  • Save aadm/424490cfbb4e7980f8d14ed9fd69dcb2 to your computer and use it in GitHub Desktop.
Save aadm/424490cfbb4e7980f8d14ed9fd69dcb2 to your computer and use it in GitHub Desktop.
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.
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)
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment