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) # unpack_headers=True slows you down here
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()
Formatted header:
In [7]:
def chunk(string, width=80):
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 [8]:
section.binary_file_header
Out[8]:
In [9]:
section.traces[0].header
Out[9]:
In [10]:
len(section.traces[0].data)
Out[10]:
In [11]:
scaled = data / 1000
scaled[np.isnan(scaled)] = 0
In [12]:
scaled
Out[12]:
In [13]:
vm = np.percentile(scaled, 99)
plt.figure(figsize=(16,8))
plt.imshow(scaled.T, cmap="Greys", vmin=-vm, vmax=vm)
plt.colorbar(shrink=0.5)
plt.show()
In [14]:
from obspy.core import Trace, Stream, UTCDateTime
from obspy.io.segy.segy import SEGYTraceHeader
In [26]:
stream = Stream()
for i, trace in enumerate(scaled):
# Make the trace.
tr = Trace(trace)
# Add required data.
tr.stats.delta = 0.004
tr.stats.starttime = 0 # Not strictly required.
# Add yet more to the header (optional).
tr.stats.segy = {'trace_header': SEGYTraceHeader()}
tr.stats.segy.trace_header.trace_sequence_number_within_line = i + 1
tr.stats.segy.trace_header.receiver_group_elevation = 0
# Append the trace to the stream.
stream.append(tr)
In [27]:
stream
Out[27]:
In [28]:
stream.write('../data/out.sgy', format='SEGY', data_encoding=5) # encode 5 for IEEE
So far we only attached metadata to the traces, but we can do more by attaching some filewide metadata, like a textual header. A SEGY file normally has a file wide text header. This can be attached to the stream object.
If this header and the binary header are not set, they will be autocreated with defaults.
In [29]:
from obspy.core import AttribDict
from obspy.io.segy.segy import SEGYBinaryFileHeader
In [30]:
# Text header.
stream.stats = AttribDict()
stream.stats.textual_file_header = '{:80s}'.format('This is the textual header.').encode()
stream.stats.textual_file_header += '{:80s}'.format('This file contains seismic data.').encode()
# Binary header.
stream.stats.binary_file_header = SEGYBinaryFileHeader()
stream.stats.binary_file_header.trace_sorting_code = 4
stream.stats.binary_file_header.seg_y_format_revision_number = 0x0100
In [31]:
import sys
stream.write('../data/out.sgy', format='SEGY', data_encoding=5, byteorder=sys.byteorder)
© Agile Geoscience 2016