Metadata-Version: 2.4
Name: seiskit
Version: 0.1.0
Summary: Exploration (reflection) seismic processing toolkit: SEG-Y I/O, gain, filtering, deconvolution, velocity analysis, NMO, stacking, and post-stack migration.
Author: Hadi Azizpour Lindi
License: MIT
License-File: LICENSE
Keywords: geophysics,migration,segy,seismic,signal-processing
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.9
Requires-Dist: click>=8.1
Requires-Dist: matplotlib>=3.6
Requires-Dist: numpy>=1.23
Requires-Dist: pandas>=1.5
Requires-Dist: pyyaml>=6.0
Requires-Dist: scipy>=1.9
Requires-Dist: segyio>=1.9
Provides-Extra: dev
Requires-Dist: build>=1.0; extra == 'dev'
Requires-Dist: pytest>=7.0; extra == 'dev'
Requires-Dist: ruff>=0.4; extra == 'dev'
Description-Content-Type: text/markdown

# seiskit — exploration seismic processing toolkit

A Python library + CLI that runs the full **reflection (exploration) seismic**
processing chain, turning raw shot gathers into an imaged subsurface section:

```
SEG-Y I/O → geometry / CMP sort → gain → filter → deconvolution
          → velocity analysis → NMO → CMP stack → post-stack migration
```

Each stage is a pure function `stage(SeismicData, **params) -> SeismicData`, so
stages compose cleanly and are individually testable. A synthetic layered-earth
generator ships with the package, so everything is runnable without proprietary
field data.

## Install

```bash
cd seismic-toolkit
pip install -e ".[dev]"
```

Requires Python ≥ 3.9. Core deps: numpy, scipy, matplotlib, pandas, segyio,
pyyaml, click.

## Quick start (CLI, end-to-end)

```bash
# 1. make a synthetic survey from a model config
seiskit synth -c configs/model.yaml -o build/raw.sgy

# 2. inspect geometry / headers
seiskit info build/raw.sgy

# 3. run the whole pipeline (geometry → gain → filter → decon → stack → migration)
seiskit run configs/process.yaml
#   -> build/migrated.sgy  +  build/qc/00_raw_shot.png, build/qc/99_result.png
```

`configs/process.yaml` is the reproducible recipe — an ordered list of stages:

```yaml
input: build/raw.sgy
output: build/migrated.sgy
qc_dir: build/qc
steps:
  - geometry
  - { name: agc, window: 0.5 }
  - { name: bandpass, f1: 8, f2: 55 }
  - { name: spiking_decon, oplen: 0.12 }
  - { name: autostack, vmin: 1600, vmax: 3000, nv: 100 }
  - { name: stolt, v: 2300 }
```

### Master pipeline — one config, toggle everything, stack after each step

`configs/process_master.yaml` keeps every stage in a single list and lets you turn each on or
off with an `enabled:` flag, load the SEG-Y **and** geometry up front, and drop a stacked-section
PNG after every enabled stage so you can see what each step did:

```yaml
input:
  segy: build/raw.sgy
  geometry: { csv: configs/geom.csv }   # merge + assign offset/CDP/elevations up front
output: build/migrated_master.sgy
qc_dir: build/qc_master
qc_stack:
  every_step: true                            # write NN_after_<stage>.png after each stage
  velocity: 2300                              # constant fallback (brute stack)
  velocity_file: build/qc_master/vfield.npz   # use the real field once it exists
steps:
  - { name: trace_edit,    enabled: true,  min_rms_frac: 0.1, max_rms_frac: 8.0 }
  - { name: tpow,          enabled: true,  power: 2.0 }
  - { name: spiking_decon, enabled: true,  oplen: 0.12 }
  - { name: bandpass,      enabled: true,  f1: 8, f2: 55 }
  - { name: demultiple,    enabled: false, q_cut: 0.04 }   # toggled off
  - { name: agc,           enabled: true,  window: 0.5 }
  - { name: autostack,     enabled: true,  save_velocity: build/qc_master/vfield.npz }
  - { name: stolt,         enabled: true,  v: 2300 }
```

The per-step QC stack picks its velocity in this order: the **velocity file** when it is present
on disk (a real per-CMP field → accurate stack), otherwise a **constant-velocity brute stack**
(fast, and comparable across steps since the velocity is held fixed); a stage whose output is
already stacked/migrated is imaged directly. `autostack.save_velocity` writes that field mid-run,
so later stages — and every stage on a second run — stack with the real velocity.

## Loading geometry onto raw field data

Raw field SEG-Y often carries only the field-record number (`fldr`/FFID) and
channel (`tracf`); the real source/receiver coordinates live in separate survey
files. `seiskit` merges them on `(fldr, tracf)` from either a flat CSV
("observer log") or SPS S/R/X files, then assigns offset/CDP (2-D, any azimuth):

```bash
# demo: make a coordinate-stripped record + sidecar geometry files
seiskit synth -c configs/model.yaml -o build/raw.sgy --emit-geometry
#   -> build/raw.sgy (coords zeroed) + build/raw.geometry.csv + build/raw.s/.r/.x

seiskit geom build/raw.sgy --sps build/raw.s,build/raw.r,build/raw.x -o build/geo.sgy
seiskit geom build/raw.sgy --csv build/raw.geometry.csv            -o build/geo.sgy
```

Or as the first pipeline stage:

```yaml
steps:
  - { name: load_geometry, sps: ["raw.s", "raw.r", "raw.x"] }   # or: { csv: raw.geometry.csv }
  - { name: agc, window: 0.5 }
  # ... rest of the pipeline
```

Unmatched traces are reported and dropped (never silently kept with bad coords).
The flat CSV needs columns `fldr,tracf,sx,sy,gx,gy` (`selev,gelev` optional). SPS
column positions default to the common SEG rev 2.1 layout and are overridable via
`column_map` (field files vary by vendor/revision).

## Velocity picking (human-like, constrained DP)

Velocity auto-picking (`proc/velocity.py`) defaults to a **constrained dynamic
program** (`pick_dp`) rather than a per-time `argmax`: it finds the whole `v(t)` path
that rides the strong semblance peaks along a *smooth, generally-increasing* trend,
penalising roughness and — crucially — downward jumps. Because semblance is a
coherency measure, a multiple flattens (high semblance) at its own slow velocity just
like a primary; only this "velocity increases with depth" prior tells them apart,
exactly as a human interpreter reasons. `autostack` builds the per-CMP field, smooths
it **laterally** across CMPs (`lateral_smooth`), then NMO-stacks. On a noisy synthetic
this lifts reflector strength ~8 % over the old argmax picker; the legacy method is
still available via `analyze(..., method="argmax")`.

## Multiple attenuation (parabolic Radon)

Multiples (energy that reverberates more than once — e.g. the marine water-bottom
bounce) mimic reflectors at wrong times. After NMO with the *primary* velocity,
primaries are flat while multiples keep a residual parabolic moveout; a parabolic
Radon transform (`proc/radon.py`) separates events by that curvature (`q`), so the
multiple zone (`|q| ≥ q_cut`) can be modelled and subtracted:

```yaml
steps:
  - geometry
  - { name: demultiple, q_cut: 0.04, qmax: 0.30 }   # per-CMP parabolic Radon
  - { name: autostack, vmin: 1600, vmax: 3000 }
```

Validated on synthetic ground truth: a curved multiple is attenuated to <40 % energy
while a flat primary is preserved (>60 %). On real data it needs a **proper velocity
model** (per-CMP picked, not a single constant) and adequate fold — on the crude
constant-velocity, low-fold Viking subset it runs but adds striping in low-fold zones
rather than a clean uplift. This is the natural pairing with interactive velocity
picking.

## Real data (Mobil AVO Viking Graben Line 12)

The toolkit is validated end-to-end on a genuine public-domain field line (2-D
marine, North Sea; released by Mobil Oil for the 1994 SEG inversion workshop). The
raw file is format-1 **IBM float** with an **EBCDIC** textual header — both handled
automatically.

```bash
# fetch a 40-shot, 30 MB byte-range subset (the full line is 749 MB)
seiskit fetch --name viking --max-bytes 29955600 -o build/viking_subset.segy

seiskit scan build/viking_subset.segy          # IBM float, EBCDIC header, geometry
seiskit run  configs/process_real.yaml         # top-mute → gain → bandpass →
                                                # decon → AGC → velan+NMO+stack
```

`scan` reports the real headers (`format: IBM float`, 4 ms, 1500 samples, 120
channels/shot, CDP + offset present). The pipeline images continuous North Sea
reflectors in the high-fold CDPs; the near CDPs are low-fold edge (only 40 shots in
the subset — use more shots or the full line for full fold). The file already
carries `cdp`/`offset`, so the real-data config processes on its own geometry (no
`geometry`/`load_geometry` step needed).

## Statics (land near-surface corrections)

Topography and the weathering layer impose source/receiver time shifts that smear
the stack. Three corrections share one surface-consistent least-squares solver
(`proc/statics.py`):

- `elevation_statics` — datum correction from `selev`/`gelev` and a replacement velocity.
- `refraction_statics` — STA/LTA first-break picking → refractor velocity → per-source
  /receiver intercept-delay statics (single refractor).
- `residual_statics` — iterative pilot cross-correlation → surface-consistent
  source/receiver shifts (the biggest stack-quality lever).

Demo (injects known statics, then corrects them):

```bash
seiskit synth -c configs/model_statics.yaml -o build/raw_stat.sgy   # with statics
seiskit run   configs/process_statics.yaml                          # elevation + residual
```

On the shipped example, enabling the statics steps makes the stacked reflectors
~**3× stronger** (more coherent) than the identical pipeline without them.

## Quick start (library)

```python
from seiskit.synthetic import default_model, Survey, generate
from seiskit.geometry import assign_geometry
from seiskit.proc import gain, filters, decon, stack, migration
from seiskit.pipeline import autostack

data = assign_geometry(generate(default_model(), Survey()))
data = gain.agc(data, window=0.5)
data = filters.bandpass(data, 8, 55)
data = decon.spiking_decon(data, oplen=0.12)
section = autostack(data, vmin=1600, vmax=3000)   # velan + NMO + CMP stack
image = migration.stolt_migrate(section, v=2300)  # post-stack migration
```

## Pipeline stages

| Stage | Function | Purpose |
|-------|----------|---------|
| `load_geometry` | `geometry_load.apply_geometry` | merge coords from CSV/SPS onto traces |
| `geometry` | `geometry.assign_geometry` | 2-D offset + CDP-bin assignment |
| `trace_edit` | `proc.edit.kill_bad_traces` | drop/zero dead, hot, non-finite traces |
| `agc`, `tpow` | `proc.gain` | amplitude recovery |
| `bandpass`, `notch`, `fk_filter`, `top_mute` | `proc.filters` | band-limit / denoise / mute |
| `spiking_decon`, `predictive_decon` | `proc.decon` | wavelet whitening / dereverb |
| `demultiple` | `proc.radon` | parabolic Radon multiple attenuation |
| `elevation_statics`, `refraction_statics`, `residual_statics` | `proc.statics` | near-surface time corrections |
| `autostack` | `pipeline.autostack` | semblance velan → NMO → CMP stack |
| `stolt`, `kirchhoff` | `proc.migration` | post-stack time migration |
| `pstm` | `proc.migration` | Kirchhoff **pre-stack** time migration |

## Status

Built milestone-by-milestone (all green):

- **M1** — data model, SEG-Y I/O, synthetic generator ✅
- **M2** — geometry/CMP sort, gain (AGC), filters (bandpass/f-k), deconvolution ✅
- **M3** — velocity analysis, NMO, CMP stack ✅
- **M4** — post-stack migration (Stolt f-k, Kirchhoff) ✅
- **M5** — CLI + YAML pipeline runner + configs ✅
- **M7** — real geometry loading (SPS + CSV), 2-D binning ✅
- **M8** — statics: elevation + refraction + residual (surface-consistent) ✅
- **M9** — real-SEG-Y I/O hardening (IBM float, EBCDIC, scan), CI, real-data run ✅
- **M10** — Kirchhoff **pre-stack** time migration (PSTM) ✅
- **M11** — parabolic Radon multiple attenuation (demultiple) ✅
- **M12** — human-like velocity picking (constrained DP) + lateral smoothing ✅
- **M13** — trace editing / kill (dead, hot, non-finite) ✅

**Deferred to later:** pre-stack depth migration, surface-consistent decon,
high-resolution/hyperbolic Radon + SRME, learned (CNN) velocity picking, and any
GUI/web front-end.

## Tests & CI

```bash
pytest        # 29 checks (synthetic ground-truth + real-format SEG-Y)
ruff check src tests
```

GitHub Actions runs ruff + pytest + a wheel/sdist build on Python 3.9–3.12.
