seismic-interpretation

Seismic Interpretation Workflow

End-to-end pipeline for seismic data analysis, from loading SEG-Y files through signal processing, rock physics modelling, and 3D visualization.

Skill Chain

segyio          obspy           bruges          disba           pyvista
[SEG-Y I/O] --> [Signal Proc] --> [Rock Physics] --> [Dispersion] --> [3D Viz]
  |               |                |                 |                |
  Load traces     Filter/FFT       AVO modelling     Surface waves    Volume render
  Read headers    Instrument resp   Fluid sub         Phase velocity   Slice display
  3D geometry     Spectral anal    Synthetics        Inversion        Horizon pick

Decision Points

Step-by-Step Orchestration

Stage 1: Data Loading (segyio)

This loader handles a complete post-stack grid, including crossline-fast trace ordering. Supply the survey's projected CRS from authoritative metadata. Header coordinate units must be length units, and the binary measurement system must identify metres or feet. Angular coordinates need reprojection first. Confirm that the sample axis represents two-way time in milliseconds for this workflow.

import segyio
import numpy as np

def load_poststack(path, crs):
    if not isinstance(crs, str) or not crs.strip():
        raise ValueError('Supply the projected survey CRS; SEG-Y coordinates alone do not identify it')
    with segyio.open(path, 'r', iline=189, xline=193) as f:
        if len(f.offsets) != 1:
            raise ValueError('Select or stack offsets before loading a post-stack cube')
        time_ms = np.asarray(f.samples, dtype=float).copy()
        if len(time_ms) < 2 or not np.all(np.isfinite(time_ms)):
            raise ValueError('Require at least two finite time samples')
        dt_ms = time_ms[1] - time_ms[0]
        if dt_ms <= 0 or not np.allclose(np.diff(time_ms), dt_ms):
            raise ValueError('Resample inconsistent time sampling before building a cube')
        to_metres = {1: 1.0, 2: 0.3048}.get(f.bin[segyio.BinField.MeasurementSystem])
        if to_metres is None:
            raise ValueError('Confirm SEG-Y measurement system: metres or feet')
        ilines, xlines = np.array(f.ilines), np.array(f.xlines)
        ii = {int(value): index for index, value in enumerate(ilines)}
        jj = {int(value): index for index, value in enumerate(xlines)}
        shape = (len(ilines), len(xlines))
        cube = np.empty((*shape, len(time_ms)), dtype=np.float32)
        easting, northing = np.empty(shape), np.empty(shape)
        seen = np.zeros(shape, dtype=bool)
        for trace_index, header in enumerate(f.header):
            i, j = ii[header[189]], jj[header[193]]
            if seen[i, j]:
                raise ValueError('Duplicate inline/crossline location')
            if header[segyio.TraceField.CoordinateUnits] != 1:
                raise ValueError('Require projected coordinates in length units')
            if (not np.isclose(header[segyio.TraceField.TRACE_SAMPLE_INTERVAL] / 1000, dt_ms)
                    or not np.isclose(header[segyio.TraceField.DelayRecordingTime], time_ms[0])):
                raise ValueError('Trace time axes differ; align them before stacking into a cube')
            scalar = header[segyio.TraceField.SourceGroupScalar]
            factor = scalar if scalar > 0 else 1 / abs(scalar) if scalar < 0 else 1.0
            easting[i, j] = header[segyio.TraceField.CDP_X] * factor * to_metres
            northing[i, j] = header[segyio.TraceField.CDP_Y] * factor * to_metres
            cube[i, j] = f.trace[trace_index]
            seen[i, j] = True
        if not np.all(seen):
            raise ValueError('Incomplete grid: retain missing traces explicitly or regrid')
    return dict(amplitude=cube, easting_m=easting, northing_m=northing,
                time_ms=time_ms, crs=crs, ilines=ilines, xlines=xlines)