# 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 → io → core → analysis →
   manager`. A layer imports only from below — never sideways, never up.
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.

---

## 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)
  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 · min_curvature · statistics · filters · arithmetic
  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            # -> Stats(mean,min,max,p10,p50,p90,std) over the Brent interval
  w.brent.phie.mean
  w.xyz(2450)            # interpolated position at MD 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 green-room SEG-Y.**
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, MIT/Apache, `-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 green-room) 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.
