Read and Write SEGY files with pysegy

This tutorial demonstrates how to read existing SEGY data shipped with the repository and how to write new files using pysegy. We use one of the example files located in the data directory.

import pysegy as seg
import fsspec
import matplotlib.pyplot as plt

Read a SEGY file

We first load the file using segy_read without specifying a filesystem.

path = '../data/overthrust_2D_shot_1_20.segy'
block = seg.segy_read(path)
print('Samples per trace:', block.fileheader.bfh.ns)
print('Number of traces:', len(block.traceheaders))
Reading SEGY file ../data/overthrust_2D_shot_1_20.segy
Loaded header ns=751 dt=4000 from ../data/overthrust_2D_shot_1_20.segy
Samples per trace: 751
Number of traces: 3300

The first trace header contains a variety of fields. Here we display a few common values.

hdr = block.traceheaders[0]
print('SourceX:', hdr.SourceX)
print('GroupX:', hdr.GroupX)
print('ns:', hdr.ns)
SourceX: 400
GroupX: 100
ns: 751

Visualise a trace

The trace data are stored in a NumPy array. The following plot shows the amplitudes of the first trace.

plt.figure(figsize=(6,3))
plt.plot(block.data[:,0])
plt.title('Trace 0')
plt.xlabel('Sample')
plt.ylabel('Amplitude')
plt.tight_layout()
plt.show()

Visualise the whole gather

Using imshow helps see the shot as an image.

seg.plot_sdata(block.data[:, :800], (.004, 25))

Write a new file

We create a small subset containing the first 20 traces and write it to a temporary file.

subset = seg.SeisBlock(block.fileheader, block.traceheaders[:20], block.data[:, :20])
out_path = 'subset.segy'
seg.segy_write(out_path, subset)
# reload to verify
out = seg.segy_read(out_path)
print('Wrote and reloaded', out_path, 'with', len(out.traceheaders), 'traces')
Writing SEGY file subset.segy
Finished writing subset.segy
Reading SEGY file subset.segy
Loaded header ns=751 dt=4000 from subset.segy
Wrote and reloaded subset.segy with 20 traces

Using fsspec filesystems

The same functions work with any filesystem that follows the fsspec interface. Here we use the local file implementation.

fs = fsspec.filesystem('file')
block_fs = seg.segy_read(path, fs=fs)
print('Read with filesystem, ns =', block_fs.fileheader.bfh.ns)
fs_path = 'subset_fs.segy'
seg.segy_write(fs_path, subset, fs=fs)
reloaded = seg.segy_read(fs_path, fs=fs)
print('Reloaded via fs, traces =', len(reloaded.traceheaders))
Reading SEGY file ../data/overthrust_2D_shot_1_20.segy
Loaded header ns=751 dt=4000 from ../data/overthrust_2D_shot_1_20.segy
Read with filesystem, ns = 751
Writing SEGY file subset_fs.segy
Finished writing subset_fs.segy
Reading SEGY file subset_fs.segy
Loaded header ns=751 dt=4000 from subset_fs.segy
Reloaded via fs, traces = 20

Scan the full dataset

We can scan all the overthrust files and access each shot without loading everything at once.

scan = seg.segy_scan('../data', 'overthrust_2D_shot_*.segy')
print('Shots found:', len(scan))
Scanning 5 files in ../data with 8 threads
ThreadPoolExecutor-16_0 scanning file ../data/overthrust_2D_shot_1_20.segy
ThreadPoolExecutor-16_1 scanning file ../data/overthrust_2D_shot_21_40.segy
Header for ../data/overthrust_2D_shot_1_20.segy: ns=751 dt=4000
ThreadPoolExecutor-16_2 scanning file ../data/overthrust_2D_shot_41_60.segy
ThreadPoolExecutor-16_3 scanning file ../data/overthrust_2D_shot_61_80.segy
Header for ../data/overthrust_2D_shot_21_40.segy: ns=751 dt=4000
ThreadPoolExecutor-16_4 scanning file ../data/overthrust_2D_shot_81_97.segy
Header for ../data/overthrust_2D_shot_41_60.segy: ns=751 dt=4000
Header for ../data/overthrust_2D_shot_81_97.segy: ns=751 dt=4000
Header for ../data/overthrust_2D_shot_61_80.segy: ns=751 dt=4000
ThreadPoolExecutor-16_4 found 17 shots in ../data/overthrust_2D_shot_81_97.segyThreadPoolExecutor-16_0 found 20 shots in ../data/overthrust_2D_shot_1_20.segy

ThreadPoolExecutor-16_3 found 20 shots in ../data/overthrust_2D_shot_61_80.segy
ThreadPoolExecutor-16_2 found 20 shots in ../data/overthrust_2D_shot_41_60.segy
ThreadPoolExecutor-16_1 found 20 shots in ../data/overthrust_2D_shot_21_40.segy
Combined scan has 97 shots
Shots found: 97
first = scan[0]
seg.plot_sdata(first)

Back to top