obspyBefore going any further, you might like to know, What is SEG-Y?. See also the articles in SubSurfWiki and Wikipedia.
We'll use the obspy seismology library to read and write SEGY data.
Technical SEG-Y documentation:
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
ls -l ../data/*.sgy
In [3]:
filename = '../data/HUN00-ALT-01_STK.sgy'
In [4]:
from obspy.io.segy.segy import _read_segy
section = _read_segy(filename)
# OPTIONS
# headonly=True — only reads the header info, then you can index in on-the-fly.
# unpack_headers=True — slows you down here and isn't really required.
In [5]:
data = np.vstack([t.data for t in section.traces])
In [6]:
plt.figure(figsize=(16,8))
plt.imshow(data.T, cmap="Greys")
plt.colorbar(shrink=0.5)
plt.show()
In [7]:
section.traces[0]
Out[7]:
In [8]:
section.textual_file_header
Out[8]:
Aargh...
OK, fine, we'll reformat this.
In [9]:
def chunk(string, width=80):
try:
# Make sure we don't have a ``bytes`` object.
string = string.decode()
except:
# String is already a string, carry on.
pass
lines = int(np.ceil(len(string) / width))
result = ''
for i in range(lines):
line = string[i*width:i*width+width]
result += line + (width-len(line))*' ' + '\n'
return result
s = section.textual_file_header.decode()
print(chunk(s))
In [10]:
section.traces[0]
Out[10]:
In [15]:
t = section.traces[0]
t.npts
Out[15]:
In [16]:
t.header
Out[16]:
Either use the small volume, or get the large dataset from Agile's S3 bucket
In [13]:
filename = '../data/F3_very_small.sgy'
# filename = '../data/Penobscot_0-1000ms.sgy'
In [14]:
from obspy.io.segy.segy import _read_segy
raw = _read_segy(filename)
In [15]:
data = np.vstack([t.data for t in raw.traces])
I happen to know that the shape of this dataset is 601 × 481.
In [17]:
_, t = data.shape
seismic = data.reshape((601, 481, t))
Note that we don't actually need to know the last dimension, if we already have two of the three dimensions. np.reshape() can compute it for us on the fly:
In [ ]:
seismic = data.reshape((601, 481, -1))
Plot the result...
In [18]:
clip = np.percentile(seismic, 99)
fig = plt.figure(figsize=(12,6))
ax = fig.add_subplot(111)
plt.imshow(seismic[100,:,:].T, cmap="Greys", vmin=-clip, vmax=clip)
plt.colorbar(label="Amplitude", shrink=0.8)
ax.set_xlabel("Trace number")
ax.set_ylabel("Time sample")
plt.show()
In [ ]:
© Agile Geoscience 2016