petekIO — an ingest tour¶

petekIO is the subsurface data layer: files in, clean validated domain objects out — surfaces, wells (trajectories / tops / logs), points, polygons — behind one GeoData substrate. This notebook authors a small synthetic project on disk, loads it once, inspects surfaces / wells / tops / zones, and round-trips the whole thing through a single .pproj file.

All data here is hand-authored to Petrel / LAS / IRAP format spec — no real data. Point GeoData at your own export and the same calls apply.

Author a synthetic project tree¶

The helper below writes a two-well field to a throwaway temp directory: a multi-bore well 25/1-1 (bore A sampled fine, bore B coarse) and a single-bore 25/1-2, plus two depth surfaces and a scattered-point CSV. It uses only petekIO's supported formats (Petrel .wellpath + tops, LAS 2.0, IRAP-classic).

In [1]:

import math, random, tempfile
from pathlib import Path
import petekio


def _write(p, body):
    p = Path(p); p.parent.mkdir(parents=True, exist_ok=True); p.write_text(body); return p


def write_wellpath(p, head, kb, stations):
    """A positioned Petrel `.wellpath` (vertical here): rows are MD X Y Z TVD ..."""
    rows = "".join(f"{md} {x} {y} {z} {md} 0 0 0 0 0 0\n" for md, x, y, z in stations)
    return _write(p, ("# WELL TRACE FROM PETREL\n"
        f"# WELL HEAD X-COORDINATE: {head[0]} (m)\n# WELL HEAD Y-COORDINATE: {head[1]} (m)\n"
        f"# WELL DATUM (KB, Kelly bushing, from MSL): {kb} (m)\n# CRS: ED50 / UTM zone 31N\n=====\n"
        "MD X Y Z TVD DX DY AZIM_TN INCL DLS AZIM_GN\n" + rows))


def write_las(p, curves, samples):
    """A LAS 2.0 file: a depth index plus one column per (mnemonic, unit)."""
    strt, stop = samples[0][0], samples[-1][0]
    hdr = "".join(f" {m}.{u} :\n" for m, u in curves)
    rows = "".join(" ".join(f"{v:.4f}" for v in row) + "\n" for row in samples)
    return _write(p, ("~Version\n VERS. 2.0 :\n WRAP. NO :\n"
        f"~Well\n STRT.M {strt} :\n STOP.M {stop} :\n STEP.M 0 :\n NULL. -999.25 :\n"
        f"~Curve\n DEPT.M :\n{hdr}~ASCII\n" + rows))


def write_tops(p, rows):
    """A Petrel multi-well tops file. `rows` = (md, surface, well)."""
    hdr = ("# Petrel well tops\nVERSION 2\nBEGIN HEADER\n"
           "X\nY\nZ\nTWT\nTWT2\nage\nMD\nPVD\nType\nSurface\nWell\nEND HEADER\n")
    body = "".join(f'1 2 -1 -999 -999 -999 {md} -1 Horizon "{s}" "{w}"\n' for md, s, w in rows)
    return _write(p, hdr + body)


def write_irap(p, xori, yori, xinc, yinc, ncol, nrow, fn):
    """An IRAP-classic ASCII surface; `fn(i, j)` -> depth (NaN = undefined)."""
    undef = 9999900.0
    xmax, ymax = xori + (ncol - 1) * xinc, yori + (nrow - 1) * yinc
    lines = [f"-996 {nrow} {xinc} {yinc}", f"{xori} {xmax} {yori} {ymax}",
             f"{ncol} 0 {xori} {yori}", "0  0  0  0  0  0  0"]
    toks = [(f"{undef:.6f}" if (v := fn(i, j)) != v else f"{v:.6f}")
            for j in range(nrow) for i in range(ncol)]
    for k in range(0, len(toks), 6):
        lines.append(" ".join(toks[k:k + 6]))
    return _write(p, "\n".join(lines) + "\n")


# The intended column, shallow -> deep: (name, md-below-datum, phi, sw, vsh).
_COLUMN = [("Top Reservoir", 0, 0.06, 0.85, 0.42), ("Upper Sand", 40, 0.24, 0.28, 0.10),
           ("Mid Sand", 80, 0.27, 0.20, 0.08), ("Lower Sand", 120, 0.19, 0.45, 0.22),
           ("Base Reservoir", 170, 0.08, 0.90, 0.45)]


def _samples(base, step, seed, phi_bias=0.0):
    """PHIE / SW / VSH down a hole: each zone's modelled values + a little noise."""
    rng = random.Random(seed); md = base - 30; out = []
    tops = [(n, base + d) for n, d, *_ in _COLUMN]
    while md <= base + 190:
        phi, sw, vsh = 0.05, 0.95, 0.5
        for n, d, p, s, v in _COLUMN:
            if md >= base + d:
                phi, sw, vsh = p, s, v
        out.append((md, max(0.0, phi + phi_bias + rng.uniform(-0.01, 0.01)),
                    min(1.0, max(0.0, sw + rng.uniform(-0.03, 0.03))),
                    max(0.0, vsh + rng.uniform(-0.02, 0.02))))
        md += step
    return out, tops


def make_field(out_dir):
    """Write a small SYNTHETIC field under `out_dir`; return the key paths.

    Two wells: `25/1-1` is multi-bore (bore A sampled fine, bore B coarse and
    biased low, so thickness-weighting visibly matters); `25/1-2` is single-bore.
    Plus two depth surfaces (a gentle dome) and a scattered-point CSV.
    """
    out = Path(out_dir); C = [("PHIE", "m3/m3"), ("SW", "m3/m3"), ("VSH", "m3/m3")]
    top_rows = []
    for bore, step, seed, bias in [("A", 0.5, 1, 0.0), ("B", 1.5, 2, -0.05)]:
        base = 2000.0
        write_wellpath(out / "wells" / "25_1-1" / f"25_1-1_{bore}.wellpath", (1000.0, 5000.0),
                       25.0, [(base - 40, 1000.0, 5000.0, -(base - 40)),
                              (base + 200, 1000.0, 5000.0, -(base + 200))])
        s, tops = _samples(base, step, seed, bias)
        write_las(out / "wells" / "25_1-1" / f"25_1-1_{bore}_logs.las", C, s)
        top_rows += [(md, n, f"25/1-1 {bore}") for n, md in tops]
    base = 2120.0
    write_wellpath(out / "wells" / "25_1-2" / "25_1-2.wellpath", (3000.0, 5200.0), 30.0,
                   [(base - 40, 3000.0, 5200.0, -(base - 40)), (base + 200, 3000.0, 5200.0, -(base + 200))])
    s, tops = _samples(base, 0.5, 3)
    write_las(out / "wells" / "25_1-2" / "25_1-2_logs.las", C, s)
    top_rows += [(md, n, "25/1-2") for n, md in tops]
    tp = write_tops(out / "tops" / "field.tops", top_rows)
    top = write_irap(out / "surfaces" / "top_res.irap", 0.0, 0.0, 50.0, 50.0, 40, 30,
                     lambda i, j: 2000.0 - 60.0 * math.exp(-(((i - 20) / 9) ** 2 + ((j - 15) / 7) ** 2)))
    base_s = write_irap(out / "surfaces" / "base_res.irap", 0.0, 0.0, 50.0, 50.0, 40, 30,
                        lambda i, j: 2060.0 - 55.0 * math.exp(-(((i - 20) / 9) ** 2 + ((j - 15) / 7) ** 2)))
    rng = random.Random(7)
    pts = "x,y,depth,poro\n" + "".join(
        f"{rng.uniform(0, 1950):.1f},{rng.uniform(0, 1450):.1f},"
        f"{2000 - 55 * math.exp(-(((rng.random() - .5) * 3) ** 2)):.2f},{rng.uniform(0.05, 0.30):.3f}\n"
        for _ in range(60))
    _write(out / "points" / "scatter.csv", pts)
    return {"root": out, "wells_dir": out / "wells", "tops": tp,
            "top_surface": top, "base_surface": base_s, "points": out / "points" / "scatter.csv"}


# Write the field into a throwaway temp directory and confirm the tree.
field = make_field(tempfile.mkdtemp())
sorted(str(p.relative_to(field["root"])) for p in field["root"].rglob("*") if p.is_file())

import math, random, tempfile from pathlib import Path import petekio def _write(p, body): p = Path(p); p.parent.mkdir(parents=True, exist_ok=True); p.write_text(body); return p def write_wellpath(p, head, kb, stations): """A positioned Petrel `.wellpath` (vertical here): rows are MD X Y Z TVD ...""" rows = "".join(f"{md} {x} {y} {z} {md} 0 0 0 0 0 0\n" for md, x, y, z in stations) return _write(p, ("# WELL TRACE FROM PETREL\n" f"# WELL HEAD X-COORDINATE: {head[0]} (m)\n# WELL HEAD Y-COORDINATE: {head[1]} (m)\n" f"# WELL DATUM (KB, Kelly bushing, from MSL): {kb} (m)\n# CRS: ED50 / UTM zone 31N\n=====\n" "MD X Y Z TVD DX DY AZIM_TN INCL DLS AZIM_GN\n" + rows)) def write_las(p, curves, samples): """A LAS 2.0 file: a depth index plus one column per (mnemonic, unit).""" strt, stop = samples[0][0], samples[-1][0] hdr = "".join(f" {m}.{u} :\n" for m, u in curves) rows = "".join(" ".join(f"{v:.4f}" for v in row) + "\n" for row in samples) return _write(p, ("~Version\n VERS. 2.0 :\n WRAP. NO :\n" f"~Well\n STRT.M {strt} :\n STOP.M {stop} :\n STEP.M 0 :\n NULL. -999.25 :\n" f"~Curve\n DEPT.M :\n{hdr}~ASCII\n" + rows)) def write_tops(p, rows): """A Petrel multi-well tops file. `rows` = (md, surface, well).""" hdr = ("# Petrel well tops\nVERSION 2\nBEGIN HEADER\n" "X\nY\nZ\nTWT\nTWT2\nage\nMD\nPVD\nType\nSurface\nWell\nEND HEADER\n") body = "".join(f'1 2 -1 -999 -999 -999 {md} -1 Horizon "{s}" "{w}"\n' for md, s, w in rows) return _write(p, hdr + body) def write_irap(p, xori, yori, xinc, yinc, ncol, nrow, fn): """An IRAP-classic ASCII surface; `fn(i, j)` -> depth (NaN = undefined).""" undef = 9999900.0 xmax, ymax = xori + (ncol - 1) * xinc, yori + (nrow - 1) * yinc lines = [f"-996 {nrow} {xinc} {yinc}", f"{xori} {xmax} {yori} {ymax}", f"{ncol} 0 {xori} {yori}", "0 0 0 0 0 0 0"] toks = [(f"{undef:.6f}" if (v := fn(i, j)) != v else f"{v:.6f}") for j in range(nrow) for i in range(ncol)] for k in range(0, len(toks), 6): lines.append(" ".join(toks[k:k + 6])) return _write(p, "\n".join(lines) + "\n") # The intended column, shallow -> deep: (name, md-below-datum, phi, sw, vsh). _COLUMN = [("Top Reservoir", 0, 0.06, 0.85, 0.42), ("Upper Sand", 40, 0.24, 0.28, 0.10), ("Mid Sand", 80, 0.27, 0.20, 0.08), ("Lower Sand", 120, 0.19, 0.45, 0.22), ("Base Reservoir", 170, 0.08, 0.90, 0.45)] def _samples(base, step, seed, phi_bias=0.0): """PHIE / SW / VSH down a hole: each zone's modelled values + a little noise.""" rng = random.Random(seed); md = base - 30; out = [] tops = [(n, base + d) for n, d, *_ in _COLUMN] while md <= base + 190: phi, sw, vsh = 0.05, 0.95, 0.5 for n, d, p, s, v in _COLUMN: if md >= base + d: phi, sw, vsh = p, s, v out.append((md, max(0.0, phi + phi_bias + rng.uniform(-0.01, 0.01)), min(1.0, max(0.0, sw + rng.uniform(-0.03, 0.03))), max(0.0, vsh + rng.uniform(-0.02, 0.02)))) md += step return out, tops def make_field(out_dir): """Write a small SYNTHETIC field under `out_dir`; return the key paths. Two wells: `25/1-1` is multi-bore (bore A sampled fine, bore B coarse and biased low, so thickness-weighting visibly matters); `25/1-2` is single-bore. Plus two depth surfaces (a gentle dome) and a scattered-point CSV. """ out = Path(out_dir); C = [("PHIE", "m3/m3"), ("SW", "m3/m3"), ("VSH", "m3/m3")] top_rows = [] for bore, step, seed, bias in [("A", 0.5, 1, 0.0), ("B", 1.5, 2, -0.05)]: base = 2000.0 write_wellpath(out / "wells" / "25_1-1" / f"25_1-1_{bore}.wellpath", (1000.0, 5000.0), 25.0, [(base - 40, 1000.0, 5000.0, -(base - 40)), (base + 200, 1000.0, 5000.0, -(base + 200))]) s, tops = _samples(base, step, seed, bias) write_las(out / "wells" / "25_1-1" / f"25_1-1_{bore}_logs.las", C, s) top_rows += [(md, n, f"25/1-1 {bore}") for n, md in tops] base = 2120.0 write_wellpath(out / "wells" / "25_1-2" / "25_1-2.wellpath", (3000.0, 5200.0), 30.0, [(base - 40, 3000.0, 5200.0, -(base - 40)), (base + 200, 3000.0, 5200.0, -(base + 200))]) s, tops = _samples(base, 0.5, 3) write_las(out / "wells" / "25_1-2" / "25_1-2_logs.las", C, s) top_rows += [(md, n, "25/1-2") for n, md in tops] tp = write_tops(out / "tops" / "field.tops", top_rows) top = write_irap(out / "surfaces" / "top_res.irap", 0.0, 0.0, 50.0, 50.0, 40, 30, lambda i, j: 2000.0 - 60.0 * math.exp(-(((i - 20) / 9) ** 2 + ((j - 15) / 7) ** 2))) base_s = write_irap(out / "surfaces" / "base_res.irap", 0.0, 0.0, 50.0, 50.0, 40, 30, lambda i, j: 2060.0 - 55.0 * math.exp(-(((i - 20) / 9) ** 2 + ((j - 15) / 7) ** 2))) rng = random.Random(7) pts = "x,y,depth,poro\n" + "".join( f"{rng.uniform(0, 1950):.1f},{rng.uniform(0, 1450):.1f}," f"{2000 - 55 * math.exp(-(((rng.random() - .5) * 3) ** 2)):.2f},{rng.uniform(0.05, 0.30):.3f}\n" for _ in range(60)) _write(out / "points" / "scatter.csv", pts) return {"root": out, "wells_dir": out / "wells", "tops": tp, "top_surface": top, "base_surface": base_s, "points": out / "points" / "scatter.csv"} # Write the field into a throwaway temp directory and confirm the tree. field = make_field(tempfile.mkdtemp()) sorted(str(p.relative_to(field["root"])) for p in field["root"].rglob("*") if p.is_file())

Out[1]:

['points/scatter.csv',
 'surfaces/base_res.irap',
 'surfaces/top_res.irap',
 'tops/field.tops',
 'wells/25_1-1/25_1-1_A.wellpath',
 'wells/25_1-1/25_1-1_A_logs.las',
 'wells/25_1-1/25_1-1_B.wellpath',
 'wells/25_1-1/25_1-1_B_logs.las',
 'wells/25_1-2/25_1-2.wellpath',
 'wells/25_1-2/25_1-2_logs.las']

Load the project once¶

GeoData(unit=...) is the substrate — petekIO never guesses the length unit. A well folder ingests every .wellpath (one per bore) and every .las (logs), auto-routing sidetracks; load_well_tops reads the Petrel tops and derives the field-wide lithostratigraphic column.

In [2]:

geo = petekio.GeoData(unit="m")
geo.load_surface("top_res", str(field["top_surface"]))
geo.load_surface("base_res", str(field["base_surface"]))
for wid in ["25/1-1", "25/1-2"]:
    geo.load_well(wid, files=str(field["wells_dir"]))
geo.load_well_tops(str(field["tops"]))

print("unit:", geo.unit)
print("wells:", [w.id for w in geo.wells.iter()])
print("strat column:", geo.strat_order)

geo = petekio.GeoData(unit="m") geo.load_surface("top_res", str(field["top_surface"])) geo.load_surface("base_res", str(field["base_surface"])) for wid in ["25/1-1", "25/1-2"]: geo.load_well(wid, files=str(field["wells_dir"])) geo.load_well_tops(str(field["tops"])) print("unit:", geo.unit) print("wells:", [w.id for w in geo.wells.iter()]) print("strat column:", geo.strat_order)

unit: m
wells: ['25/1-1', '25/1-2']
strat column: ['Top Reservoir', 'Upper Sand', 'Mid Sand', 'Lower Sand', 'Base Reservoir']

Inspect wells & sidetracks¶

Per-bore access is first-class. 25/1-1 is multi-bore, so choose a bore (or set a default) before the top-level accessors; each bore is positioned by its own trajectory.

In [3]:

w = geo.well("25/1-1")
print("25/1-1 bores:", w.bores(), " multibore:", w.is_multibore)
print("25/1-2 bores:", geo.well("25/1-2").bores())

bore = w.sidetrack("A")
print("bore A curves:", bore.mnemonics(), " MD range:", bore.md_range())
print("position at MD 2100:", bore.xyz(2100.0), " TVD:", bore.tvd(2100.0))

w = geo.well("25/1-1") print("25/1-1 bores:", w.bores(), " multibore:", w.is_multibore) print("25/1-2 bores:", geo.well("25/1-2").bores()) bore = w.sidetrack("A") print("bore A curves:", bore.mnemonics(), " MD range:", bore.md_range()) print("position at MD 2100:", bore.xyz(2100.0), " TVD:", bore.tvd(2100.0))

25/1-1 bores: ['', 'A', 'B']  multibore: True
25/1-2 bores: ['']
bore A curves: ['PHIE', 'SW', 'VSH']  MD range: (1960.0, 2200.0)
position at MD 2100: (1000.0, 5000.0, -2075.0)  TVD: 2075.0

Surfaces — sample, volumetrics, thickness¶

A Surface is a regular gridded layer; ops return new surfaces (immutable). sample is a NaN-aware bilinear read; area_below(depth) sums the cell-area shallower-or-equal to a depth (the GRV-style query); Surface.thickness is the isochore between two horizons.

In [4]:

top = geo.surface("top_res")
base = geo.surface("base_res")
s = top.stats()
print(f"top_res: {top.ncol}x{top.nrow} nodes  mean={s.mean:.1f}  p10={s.p10:.1f}  p90={s.p90:.1f}")
print("sample at (1000, 750):", top.sample(1000.0, 750.0))
print("area shallower than 1990 m:", top.area_below(1990.0))

thick = petekio.Surface.thickness(top, base, clamp_zero=True)
print("mean isochore thickness:", round(thick.stats().mean, 2))

top = geo.surface("top_res") base = geo.surface("base_res") s = top.stats() print(f"top_res: {top.ncol}x{top.nrow} nodes mean={s.mean:.1f} p10={s.p10:.1f} p90={s.p90:.1f}") print("sample at (1000, 750):", top.sample(1000.0, 750.0)) print("area shallower than 1990 m:", top.area_below(1990.0)) thick = petekio.Surface.thickness(top, base, clamp_zero=True) print("mean isochore thickness:", round(thick.stats().mean, 2))

top_res: 40x30 nodes  mean=1990.1  p10=1967.8  p90=1999.9
sample at (1000, 750): 1940.0
area shallower than 1990 m: 897500.0
mean isochore thickness: 60.82

Grid scattered points into a surface¶

Load scattered (x, y, z) data and grid it onto a lattice — "minimum_curvature" (Briggs biharmonic) honours the points as hard constraints.

In [5]:

pts = petekio.PointSet.load_csv(str(field["points"]), x="x", y="y", z="depth")
pb = pts.bbox()
print(f"{len(pts)} points; x[{pb.xmin:.0f}, {pb.xmax:.0f}] y[{pb.ymin:.0f}, {pb.ymax:.0f}]")

nx = ny = 25
geom = petekio.GridGeometry(
    xori=pb.xmin, yori=pb.ymin,
    xinc=(pb.xmax - pb.xmin) / (nx - 1), yinc=(pb.ymax - pb.ymin) / (ny - 1),
    ncol=nx, nrow=ny,
)
gridded = pts.to_surface(geom, "minimum_curvature")
print("gridded surface mean depth:", round(gridded.stats().mean, 1))

pts = petekio.PointSet.load_csv(str(field["points"]), x="x", y="y", z="depth") pb = pts.bbox() print(f"{len(pts)} points; x[{pb.xmin:.0f}, {pb.xmax:.0f}] y[{pb.ymin:.0f}, {pb.ymax:.0f}]") nx = ny = 25 geom = petekio.GridGeometry( xori=pb.xmin, yori=pb.ymin, xinc=(pb.xmax - pb.xmin) / (nx - 1), yinc=(pb.ymax - pb.ymin) / (ny - 1), ncol=nx, nrow=ny, ) gridded = pts.to_surface(geom, "minimum_curvature") print("gridded surface mean depth:", round(gridded.stats().mean, 1))

60 points; x[45, 1904] y[39, 1440]
gridded surface mean depth: 1973.1

Per-zone stats, in stratigraphic order¶

A top resolves an interval [top_md, base_md). zone_stats returns per-zone Stats in lithostratigraphic order (not measured-depth order).

In [6]:

for zone, st in bore.zone_stats("PHIE"):
    print(f"  {zone:16s}  mean PHIE = {st.mean:.3f}  (n={st.count})")

for zone, st in bore.zone_stats("PHIE"): print(f" {zone:16s} mean PHIE = {st.mean:.3f} (n={st.count})")

  Top Reservoir     mean PHIE = 0.060  (n=80)
  Upper Sand        mean PHIE = 0.241  (n=80)
  Mid Sand          mean PHIE = 0.272  (n=80)
  Lower Sand        mean PHIE = 0.191  (n=100)
  Base Reservoir    mean PHIE = 0.081  (n=41)

Persist the project — one .pproj file¶

The whole project serialises to a single structured file: atomic to write, inspectable without a full load, and re-openable into an equivalent GeoData.

In [7]:

import os, tempfile
pproj = os.path.join(tempfile.mkdtemp(), "field.pproj")
geo.save(pproj)

info = petekio.GeoData.inspect(pproj)          # manifest only — no full decode
print("saved:", os.path.basename(pproj), f"({os.path.getsize(pproj)} bytes)")
print("unit:", info["unit"], " elements:", info["elements"])

reopened = petekio.GeoData.open(pproj)
print("re-opened wells:", [w.id for w in reopened.wells.iter()])

import os, tempfile pproj = os.path.join(tempfile.mkdtemp(), "field.pproj") geo.save(pproj) info = petekio.GeoData.inspect(pproj) # manifest only — no full decode print("saved:", os.path.basename(pproj), f"({os.path.getsize(pproj)} bytes)") print("unit:", info["unit"], " elements:", info["elements"]) reopened = petekio.GeoData.open(pproj) print("re-opened wells:", [w.id for w in reopened.wells.iter()])

saved: field.pproj (21861 bytes)
unit: Metres  elements: [('surface', 'top_res'), ('surface', 'base_res'), ('well', '25/1-1'), ('well', '25/1-2')]
re-opened wells: ['25/1-1', '25/1-2']

Recap¶

One GeoData, loaded once, gave us surfaces (sample / volumetrics / thickness), scattered-point gridding, multi-bore wells positioned along their trajectories, per-zone stats in stratigraphic order, and single-file persistence — with no per-file parsing in our own code. The companion notebook 02_well_analysis goes deeper on zone tables, net-cutoff sweeps, and the log-correlation view.