# petekIO — subsurface data model & IO library (build spec)
> Repo: `Koding/Rust/petekIO` · crate name: `petekio` (lowercase for
> crates.io — confirm availability). Sibling to `logSuite`. This is the handover
> spec for the build agent.
A standalone, reusable **Rust** library (with optional **PyO3** Python bindings)
that is the complete **input data model** for subsurface work: surfaces, wells
(trajectories/tops/logs), points, polygons — with loading, calculations,
interpolation, filters, and statistics built in. It fills a real gap (no Rust
crate does this; xtgeo/welly are Python-only) and is the data foundation that
apps like SimulatoRS consume so they do **zero** parsing/interpolation themselves.
**Litmus test for what belongs here:** it must be useful to anyone working with
subsurface data — a geologist, a seismic interpreter, a petrophysicist — **not
just a reservoir simulator.** Nothing reservoir-specific (wireframes, grids,
volumetrics) lives here; that's the consumer's job.
---
## Design constitution (adopted from logSuite)
1. **Strictly layered, one-way deps.** `foundation → algorithms → io → core →
analysis → manager`. A layer imports only from below — never sideways, never
up. (`algorithms` depends only on `foundation`; `io` and `algorithms` are
siblings above it.)
2. **A manager substrate.** Load once into a `GeoData` project; everything reads
from it. **No per-item loops** — operations broadcast across the collection.
3. **Domain objects carry their operations** (arithmetic, filters, interpolation,
stats) as methods/traits — fluent and chainable.
4. **Views** = read-only filtered subsets (`project.wells.filter(...)`).
5. **Open/closed.** Extend by adding new readers/operations/artifacts, not by
editing existing types.
6. **Compartmentalized — split the elephant.** One module/topic, one
type/responsibility. Soft limits: module ≲600 lines, type ≲300, method ≲50.
7. **Minimal public surface.** Re-export only what users need.
8. **Rust core + thin PyO3.** All logic in Rust; bindings only marshal. The
Python API mirrors the Rust API in the fluent logSuite style.
9. **Algorithms are isolated, QC-able, discipline-grouped kernels.** High-value
numeric / geostatistical routines live in **`algorithms/`, grouped by
discipline** (`wells`, `grids`, …) as **pure, type-light functions** —
primitives + `foundation` types in/out, no domain-object (`Surface`/`Well`/…)
or IO coupling. Domain types call into them; a kernel's math has **one home**
(never duplicated across call sites). Rationale: each kernel is trivial to
**QC in isolation** (analytic tests on raw numbers), and a kernel that proves
high-value is a cheap **lift-and-shift into the external `petekTools`
library** (the module mirrors its type-light boundary). Don't inline a formula
in a domain type.
---
## Architecture
A Cargo workspace (or one crate with layered modules — the dep arrow must stay
scannable). Composes the author's own crates + the Rust geo ecosystem; does NOT
reinvent them.
```
petekio/
foundation/ errors · units · geometry (Point2/Point3, BBox, GridGeometry)
algorithms/ pure numeric kernels, grouped by discipline: wells (min-curvature survey) · grids (gridding) — type-light, QC-able, petekTools-offloadable
io/ irap · zmap · csv · las (←las_rs) · excel (←calamine) · survey · tops · vector (←geozero)
core/ Surface · Well→Sidetrack→Trajectory · Log · Top · PointSet · PolygonSet (+ operation traits)
analysis/ resample · statistics · filters · arithmetic · model-ready-inputs pipeline
manager/ GeoData (the substrate) · Views
py/ PyO3 bindings (optional feature)
```
**Dependencies (compose, don't reinvent):**
- `las_rs` (own) — LAS logs. `calamine` — Excel/tabular (`.xlsx`/`.xls`).
- `geo` / `geozero` — polygon/point IO (GeoJSON/shapefile/WKT) + 2-D predicates.
- `ndarray` — the array backbone for surfaces (efficient, contiguous, BLAS-ready).
- `rstar` — spatial index (nearest cell / point-in-set) when needed.
- `serde` — (de)serialization of the model.
- *(Phase 3, cubes only)* `giga-segy` — SEG-Y reader (rev2, active); `ibmfloat`
fallback codec. ZGY/OpenVDS have no Rust → FFI, deferred.
**Greenfield (the genuinely missing parts):** the `Surface` type + IRAP/ZMAP
readers, the `Well/Sidetrack/Trajectory` model + minimum-curvature, surface ops,
resampling, the tops→interval→stats machinery, the `GeoData` substrate.
---
## Conventions
- **Coordinate system:** project (x=Easting, y=Northing, z=depth increasing
**downward**), feet OR metres — carried as a `unit` on the project; the library
never guesses, conversions live in `foundation/units`.
- **Undefined values:** `f64::NAN` in arrays (and a fast `is_defined` mask).
Arithmetic propagates NaN; stats skip NaN.
- **Storage precision:** `f64` default; a `f32` feature/typedef for memory-bound
large surfaces. Arrays are contiguous `ndarray::Array2`.
- **Errors:** one `GeoError` enum (`thiserror`); `Result<T, GeoError>` everywhere.
- **Immutability of ops:** arithmetic/resample return **new** objects; mutation is
explicit (`set_*`).
---
## Core data model
### `Stats` — the universal aggregation result
Returned by every reduction (a surface's values, a log over an interval, a point
attribute). The thing `well.brent.ntg` evaluates to.
```rust
pub struct Stats {
pub count: usize,
pub mean: f64, pub min: f64, pub max: f64, pub std: f64, pub sum: f64,
pub p10: f64, pub p50: f64, pub p90: f64, // default percentiles
}
impl Stats { pub fn percentile(&self, p: f64) -> f64; } // arbitrary p
// Built via a builder that supports WEIGHTING (by interval length, or by another
// curve, e.g. pore-volume-weighted Sw): Stats::weighted(values, weights).
```
Python: `s.mean`, `s.p50`, `s.percentile(0.25)`.
### `Surface` — a regular gridded surface (the workhorse)
A rotated regular grid (IRAP/RMS model) holding a primary value layer plus named
**attribute** layers (thickness, seismic, time, …) on the same geometry.
```rust
pub struct GridGeometry { // regular, rotatable
pub xori: f64, pub yori: f64, // origin
pub xinc: f64, pub yinc: f64, // node spacing
pub ncol: usize, pub nrow: usize,
pub rotation_deg: f64, // CCW from East
pub yflip: bool,
}
pub struct Surface {
pub geom: GridGeometry,
values: Array2<f64>, // primary (e.g. depth); NaN = undefined
attributes: IndexMap<String, Array2<f64>>, // thickness, seismic_amp, twt, ...
}
```
**Required behaviour (maps to the requested features):**
- **Load:** `Surface::load_irap_classic(path)` — **FIRST format to implement.**
(Then `load_irap_binary`, `load_zmap`, `load_xyz_csv`.)
- **Attributes:** `surface.attr("thickness") -> &Array2` ; `set_attr(name, arr)` ;
`as_attr_surface("seismic") -> Surface` (promote an attribute for ops).
- **Scalar arithmetic + math:** operator overloads (`&s + 100.0`, `&a * &b`) and
methods on the active layer: `.ln() .log10() .exp() .powf(n) .sqrt() .abs()
.clamp_min(x) .clamp(lo,hi)`. Each returns a new `Surface`.
- **Surface↔surface:** `a.minus(&b)`, `a.times(&b)`, … require equal geometry
(else `GeoError::GeometryMismatch` — caller resamples first). Convenience:
`Surface::thickness(top, base, clamp_zero: bool) -> Surface`.
- **Statistics:** `surface.stats() -> Stats` (over defined nodes).
- **Area / volume:** `surface.area_below(depth) -> f64` (areal extent of nodes
with value ≤ depth × cell area — the GRV-style query); `area_above`,
`volume_between(&other) -> f64`, `hypsometry()` (area-vs-depth curve).
- **Resample (native):** `surface.resample(&target: GridGeometry) -> Surface`
(bilinear; NaN-aware). `target` may come from another surface or a grid's areal
geometry. Also `sample(x, y) -> Option<f64>` (point query, bilinear).
Python ergonomics:
```python
top = geo.load_surface("top.irap")
base = geo.load_surface("base.irap")
thick = petekio.Surface.thickness(top, base, clamp_zero=True) # or (base - top).clamp_min(0)
trend = top.attr("seismic").ln()
top.stats.p50
top.area_below(8240) # ft² below the OWC
ongrid = top.resample(grid_geom) # bilinear onto a target geometry
```
### `Well` → `Sidetrack` → `Trajectory` (+ tops + logs)
```rust
pub struct Well {
pub id: String,
pub head: (f64, f64), // wellhead x, y
pub kb: f64, // KB elevation / air gap — the MD datum
sidetracks: IndexMap<String, Sidetrack>, // "" = main (default), "a","b",...
}
pub struct Sidetrack {
pub label: String,
trajectories: Vec<Trajectory>, // a sidetrack may hold several survey versions
active: usize, // exactly one is active
logs: IndexMap<String, Log>, // assigned to THIS sidetrack
tops: Vec<Top>, // assigned to THIS sidetrack
}
pub enum TrajectoryInput { // construction modes, all normalized internally
Xyz(Vec<[f64;3]>), // explicit positions
MdIncAzi(Vec<(f64,f64,f64)>), // survey stations → minimum-curvature
Stations(Vec<Station>), // raw survey stations
Hold { from: Station, to_md: f64 }, // constant inc/azi segment
Steer { from: Station, build: f64, turn: f64, to_md: f64 }, // build/turn rates
}
pub struct Trajectory { /* normalized */ path: PositionedPath } // md → (x,y,z,tvd)
pub struct Top { pub name: String, pub md: f64 } // entry MD; interval base = next top's MD (or TD)
pub struct Log { pub mnemonic: String, md: Array1<f64>, values: Array1<f64>, unit: String }
```
**Required behaviour:**
- **Construction:** `Well::new(id, head, kb)`; `well.sidetrack_mut(label)` (creates
`""`/main lazily); `st.add_trajectory(TrajectoryInput) -> &mut Trajectory`
(newest becomes active unless told otherwise); `st.set_active(i)`.
- **Trajectory types:** every `TrajectoryInput` variant normalizes to a
`PositionedPath` — `MdIncAzi`/`Stations` via **minimum-curvature** (dogleg ratio
factor), `Xyz` directly, `Hold`/`Steer` by integrating the segment. From the
wellhead `head` + datum `kb`.
- **Interpolation (native):** on the active trajectory, exposed on the well:
`well.xyz(md) -> [f64;3]`, `well.tvd(md) -> f64`, `well.position_at(md)`,
`well.md_at_tvd(tvd)`. Arc/linear interpolation between stations.
- **Logs (native):** `well.log("PHIE")` → a `LogView` with `.filter(pred)`,
`.resample(step)`, `.stats() -> Stats`, `.values()`, `.at_md(md)`. Positioned in
3-D via the active trajectory (`.xyz()`).
- **Tops → interval → stats (native, the key ergonomic):** a top defines a depth
**interval** (its MD → the next top's MD). `well.top("Brent")` → an `Interval`;
`interval.log("NTG").stats() -> Stats`. The headline Python ergonomic is dynamic
attribute access:
```python
w.brent.ntg w.brent.phie.mean
w.xyz(2450) ```
(Rust: `well.top("Brent")?.log("NTG")?.stats()`. PyO3 `__getattr__` resolves
`.<top>` then `.<log>`.)
### `PointSet` (simpler)
```rust
pub struct PointSet { coords: Array2<f64> /* N×3 */, attributes: IndexMap<String, Array1<f64>> }
```
- Load: xyz/CSV, IRAP points, GeoJSON (`geozero`). Ops: `filter(pred)`,
`attr(name)`, `stats(attr) -> Stats`, `bbox()`, `nearest(x,y)` (`rstar`),
`to_surface(geom, method)` (grid scattered points → a `Surface`).
### `PolygonSet` (simpler)
```rust
pub struct PolygonSet { polygons: Vec<Polygon /* geo::Polygon + optional z */> }
```
- Load: GeoJSON/shapefile/WKT (`geozero`), IRAP polygons. Ops: `contains(x,y)`,
`area()`, `clip(&surface)`, `bbox()`. Used for boundaries + (later) fault traces.
### `Cube` — 3-D regular volume (seismic & inversions) — *designed, deferred*
A seismic survey or an inversion result is a 3-D regular volume: a rotatable
(inline × crossline) lattice with a vertical axis (time or depth). An inversion
(acoustic impedance, Vp/Vs, a porosity cube) is the **same type** with a
different value — like `Surface` attributes but in 3-D. (Design distilled from
Equinor `xtgeo`/`segyio`/`open-zgy`/`OpenVDS` + the Rust `giga-segy` crate.)
**Geometry — xtgeo's model + one improvement.**
```rust
pub struct CubeGeometry {
pub areal: GridGeometry, // inline/xline (reuses the surface geometry: origin/inc/rotation/yflip)
pub zori: f64, pub zinc: f64, pub nz: usize,
pub ilines: Vec<i32>, pub xlines: Vec<i32>, // explicit line LABELS (segy access is by label, not index)
pub domain: Domain, // Time | Depth (xtgeo omits this — we keep it explicit)
}
pub struct Cube {
pub geom: CubeGeometry,
pub property: PropertyKind, // Amplitude | AcousticImpedance | VpVs | Porosity | ... + unit
store: CubeStore, // dense (small) OR brick-tiled (large) — see below
dead: Bitmask, // dead/undefined traces (xtgeo's traceidcodes)
}
```
One `Cube` carries seismic **and** inversion alike; the *meaning* lives in
`property`/unit metadata (the improvement over xtgeo, which leaves `values`
semantics-agnostic).
**Storage — the key architecture decision (do NOT copy xtgeo's whole-cube-in-RAM).**
xtgeo holds the entire `(ncol,nrow,nlay)` f32 array in memory — fine for small
cubes, the **ceiling to beat** for GB–TB volumes. ZGY / OpenVDS / seismic-zfp all
converge on the same answer, which we adopt:
> **Store the cube as a grid of fixed-size compressed bricks (~64³ samples) plus
> an offset lookup table — never a naïve contiguous N-D array.** A slice /
> sub-volume / horizon-extraction query then fetches only the few intersecting
> bricks (cheap random access, cloud-friendly ranged reads, per-brick
> compression, optional LOD pyramid).
`CubeStore` is an enum: `Dense(Array3<f32>)` for small cubes; `Bricked { bricks,
offset_table, .. }` for large. The `CubeGeometry` is the thin addressing layer
over the bricks. (Start with `Dense` + a seismic-zfp-style ZFP brick model;
full VDS/LOD is later.)
**IO — depend on `giga-segy`, do not reimplement SEG-Y from scratch.**
SEG-Y is a public standard (no IP upside to re-deriving) and quirk-heavy (IBM
floats, endianness, **configurable byte locations** — 189/193 default, never
assume — headers, rev1/rev2). The Rust crate **`giga-segy`** (GiGainfosystems,
rev2, active, permissively licensed, `-core`/`-in`/`-out` split that already matches our
layering) covers it. Wrap it behind a `SegyReader` trait in `io/`. **Caveat:
verify IBM-float (format code 1) round-trips on a real file before committing;**
`ibmfloat` crate as a fallback codec. Mirror segyio's two-layer split (lazy
trace reader + geometry/format detection) and **lazy/streaming by default**,
mmap opt-in. ZGY / OpenVDS have **no Rust** (C++/Python) → FFI, deferred.
**The link to the rest (the valuable op):** `surface.slice_cube(&cube, sampling)`
samples a cube along a horizon → a **surface attribute** (amplitude/impedance on
the top); `slice_cube_window(min|max|rms|mean over a z-window)` for windowed
attributes (xtgeo's `slice_cube_window`). This *is* the "seismic trend on a
surface" the MVP already supports as a surface attribute — the cube is just the
full-volume source. Also `cube.sample(x,y,z)`, `cube.inline(label)/zslice(z) ->
Surface`, `cube.stats()`, `cube.resample(&CubeGeometry)`.
**Domain conversion** (time↔depth via a velocity model) is heavier — deferred.
### `GeoData` — the manager substrate
Load everything once; named + collection access; views; broadcast.
```rust
pub struct GeoData { unit: Unit, surfaces: IndexMap<String, Surface>,
wells: IndexMap<String, Well>, points: ..., polygons: ... }
impl GeoData {
pub fn load_surface(&mut self, path) -> Result<&Surface>; // fluent (returns ref / self)
pub fn load_well(&mut self, dir_or_files, head, kb) -> Result<&Well>;
pub fn surface(&self, name) -> Option<&Surface>;
pub fn well(&self, id) -> Option<&Well>;
pub fn wells(&self) -> WellsView; // broadcastable + filterable
}
```
```python
geo = petekio.GeoData(unit="ft")
geo.load_surface("top.irap"); geo.load_well("wells/A1/", wellhead=(x,y), kb=82)
geo.wells.filter(field="Gullfaks").tops("Brent").ntg # broadcast → Stats per well
```
---
## IO formats — phasing
- **Phase 1 (MVP):** IRAP classic (ASCII) surfaces · LAS logs (`las_rs`) ·
deviation-survey CSV · tops CSV · scattered xyz/CSV points.
- **Phase 2:** IRAP binary (`.gri`, Fortran-record/byte-swapped — validate vs
`xtgeo`) · ZMAP+ surfaces · GeoJSON/shapefile polygons (`geozero`) · Excel
(`calamine`).
- **Phase 3 (3-D volumes):** SEG-Y (seismic & inversion cubes) via **`giga-segy`**
(depend, don't reinvent) behind a `SegyReader` trait; brick-tiled storage, not
load-all. Then OpenVDS / ZGY (no Rust → FFI, later).
- **Deferred:** GOCAD/SKUA · DLIS (bridge to `dlisio`) · RESQML · WITSML ·
time↔depth conversion.
Every reader validates on load (geometry sane, monotonic MD, units present) and
returns a typed error, not a panic.
---
## Build phasing (for the agent)
1. `foundation` (errors, units, `GridGeometry`, geometry primitives) + `Stats`.
2. `Surface` + IRAP-classic reader + scalar/math ops + `stats`/`area_below` +
`resample`/`sample`. (Golden tests: round-trip a known IRAP file; bilinear
resample vs hand calc; `area_below` vs analytic.)
3. `Trajectory` (minimum-curvature) + `Well`/`Sidetrack` + `well.xyz(md)`.
(Golden: a worked deviation survey; vertical-well degenerate case.)
4. `Log` + `Top`→interval + the `well.<top>.<log> -> Stats` ergonomic.
5. `PointSet` + `PolygonSet` (basic).
6. `GeoData` substrate + views/broadcast.
7. PyO3 bindings mirroring the above.
## Non-goals (keep it reusable)
- No reservoir concepts (wireframe, corner-point grid, volumetrics, MC) — those
belong to consumers (e.g. SimulatoRS).
- No plotting/visualization (a separate concern; data only).
- No format completeness for its own sake — add a reader when a real need exists.
## How SimulatoRS consumes it (the point)
SimulatoRS does **no** parsing or interpolation. Its `srs-data` adapter reads a
`GeoData`: surfaces → horizons, wells' tops/logs → control points + upscaled
properties, polygons → boundary/faults — and assembles the reservoir `Wireframe`.
The library owns the entire input-data lifecycle.