01 - Stack Model From Scatter¶

This notebook shows the simplest petekStatic role-binding path: make or point at a project tree, load it, inspect Project.inventory(), then bind the project names to model roles before building a stack model through the peteksim facade.

To use your own export, edit only the data-source cell and the role-binding literals. The synthetic manifest checks are isolated in the final cell, so they naturally skip for real data.

1. Imports¶

Only the library imports live here.

In [1]:

from pathlib import Path
import tempfile

import petektools as pt
import peteksim as ps

from pathlib import Path import tempfile import petektools as pt import peteksim as ps

2. Generate Synthetic Data Or Point At Your Export¶

The runnable path creates a synthetic Petrel-style export. For your own data, set USE_SYNTHETIC = False and edit the three PROJECT_* values.

In [2]:

# Keep this True for the runnable example.
# Set it to False and edit PROJECT_ROOT / PROJECT_CRS / PROJECT_ALIASES for your export.
USE_SYNTHETIC = True

if USE_SYNTHETIC:
    synthetic_manifest = pt.synth_asset(
        tempfile.mkdtemp(prefix="petekstatic_stack_"),
        surfaces_as_points=True,
        seed=20260704,
        ncol=41,
    )
    PROJECT_ROOT = Path(synthetic_manifest["root"])
    PROJECT_CRS = synthetic_manifest["crs"]
    PROJECT_ALIASES = synthetic_manifest["aliases"]
else:
    synthetic_manifest = None
    PROJECT_ROOT = Path("/replace/with/your/export/root")
    PROJECT_CRS = "REPLACE WITH YOUR CRS LABEL"
    PROJECT_ALIASES = {
        "PHIE_2025": "PORO",
        "NTG_PhieLam_2025": "NTG",
        "SW_2025": "SW",
    }

print("project root:", PROJECT_ROOT)

# Keep this True for the runnable example. # Set it to False and edit PROJECT_ROOT / PROJECT_CRS / PROJECT_ALIASES for your export. USE_SYNTHETIC = True if USE_SYNTHETIC: synthetic_manifest = pt.synth_asset( tempfile.mkdtemp(prefix="petekstatic_stack_"), surfaces_as_points=True, seed=20260704, ncol=41, ) PROJECT_ROOT = Path(synthetic_manifest["root"]) PROJECT_CRS = synthetic_manifest["crs"] PROJECT_ALIASES = synthetic_manifest["aliases"] else: synthetic_manifest = None PROJECT_ROOT = Path("/replace/with/your/export/root") PROJECT_CRS = "REPLACE WITH YOUR CRS LABEL" PROJECT_ALIASES = { "PHIE_2025": "PORO", "NTG_PhieLam_2025": "NTG", "SW_2025": "SW", } print("project root:", PROJECT_ROOT)

project root: /var/folders/7h/n6cl9kt96g79hz0s16q5vfxc0000gn/T/petekstatic_stack_6zsycu8v

3. Load Data And Inspect Inventory¶

Loading is separate from synthetic generation. Inspect the inventory before changing the role-binding cell.

In [3]:

project = ps.Project.load(
    str(PROJECT_ROOT),
    settings=ps.LoadSettings(crs=PROJECT_CRS, aliases=PROJECT_ALIASES),
)

print("Inventory before role binding:")
print(project.inventory())

project = ps.Project.load( str(PROJECT_ROOT), settings=ps.LoadSettings(crs=PROJECT_CRS, aliases=PROJECT_ALIASES), ) print("Inventory before role binding:") print(project.inventory())

Inventory before role binding:
Inventory(surfaces=1, polygons=2, points=14, wells=8, tops=12, merged=0, skipped=0)

4. Bind Project Names To Model Roles¶

This is the user-facing replacement point: outline, horizon order, zones, subzones, contacts, and properties are explicit literals.

In [4]:

# These are the user-replaceable role bindings.
# For a real export, replace the string literals with names from Project.inventory().
OUTLINE = "ModelEdge"

HORIZON_NAMES = [
    "H0",
    "H1",
    "H2",
    "H3",
    "H4",
    "H5",
    "H6",
]

ZONE_NAMES = [
    "Z0",
    "Z1",
    "Z2",
    "Z3",
    "Z4",
    "Z5",
]

# The synthetic asset carries H3b as a tops-only internal split in Z3.
# Replace or empty this mapping for your own stratigraphy.
SUBZONE_BINDINGS = {
    "Z3": ["H3b"],
}

if synthetic_manifest is not None:
    CONTACTS_BY_ZONE = {
        "Z2": {"owc": synthetic_manifest["contacts"]["owc_z2"]},
        "Z4": {
            "goc": synthetic_manifest["contacts"]["goc_z4"],
            "fwl": synthetic_manifest["contacts"]["fwl_z4"],
        },
    }
else:
    # Replace these contact depths with your project contacts, positive-down metres.
    CONTACTS_BY_ZONE = {
        "Z2": {"owc": 2440.0},
        "Z4": {"goc": 2410.0, "fwl": 2455.0},
    }

horizons = ps.Horizons(
    *(ps.hz(name) for name in HORIZON_NAMES),
    zones=ZONE_NAMES,
    ties=ps.TieSettings(method="convergent"),
    gridding=ps.Gridding(collapse=True),
)
subzones = ps.Subzones({
    zone: ps.splits(*split_names, conformity="proportional")
    for zone, split_names in SUBZONE_BINDINGS.items()
})
layering = ps.Layering(nk=2)
contacts = ps.Contacts(CONTACTS_BY_ZONE)

props = ps.Props(
    ps.Prop(
        "PORO",
        net_only=True,
        propagate=ps.Propagate(
            variogram=ps.variogram("spherical", 1500.0),
            seed=11,
        ),
    ),
    ps.Prop(
        "NTG",
        propagate=ps.Propagate(
            variogram=ps.variogram("spherical", 1500.0),
            seed=12,
        ),
    ),
)

print("outline:", OUTLINE)
print("horizons:", " -> ".join(HORIZON_NAMES))
print("zones   :", " -> ".join(ZONE_NAMES))
print("contacts:")
for zone in ZONE_NAMES:
    print(f"  {zone}: {contacts.for_zone(zone) or '-'}")

# These are the user-replaceable role bindings. # For a real export, replace the string literals with names from Project.inventory(). OUTLINE = "ModelEdge" HORIZON_NAMES = [ "H0", "H1", "H2", "H3", "H4", "H5", "H6", ] ZONE_NAMES = [ "Z0", "Z1", "Z2", "Z3", "Z4", "Z5", ] # The synthetic asset carries H3b as a tops-only internal split in Z3. # Replace or empty this mapping for your own stratigraphy. SUBZONE_BINDINGS = { "Z3": ["H3b"], } if synthetic_manifest is not None: CONTACTS_BY_ZONE = { "Z2": {"owc": synthetic_manifest["contacts"]["owc_z2"]}, "Z4": { "goc": synthetic_manifest["contacts"]["goc_z4"], "fwl": synthetic_manifest["contacts"]["fwl_z4"], }, } else: # Replace these contact depths with your project contacts, positive-down metres. CONTACTS_BY_ZONE = { "Z2": {"owc": 2440.0}, "Z4": {"goc": 2410.0, "fwl": 2455.0}, } horizons = ps.Horizons( *(ps.hz(name) for name in HORIZON_NAMES), zones=ZONE_NAMES, ties=ps.TieSettings(method="convergent"), gridding=ps.Gridding(collapse=True), ) subzones = ps.Subzones({ zone: ps.splits(*split_names, conformity="proportional") for zone, split_names in SUBZONE_BINDINGS.items() }) layering = ps.Layering(nk=2) contacts = ps.Contacts(CONTACTS_BY_ZONE) props = ps.Props( ps.Prop( "PORO", net_only=True, propagate=ps.Propagate( variogram=ps.variogram("spherical", 1500.0), seed=11, ), ), ps.Prop( "NTG", propagate=ps.Propagate( variogram=ps.variogram("spherical", 1500.0), seed=12, ), ), ) print("outline:", OUTLINE) print("horizons:", " -> ".join(HORIZON_NAMES)) print("zones :", " -> ".join(ZONE_NAMES)) print("contacts:") for zone in ZONE_NAMES: print(f" {zone}: {contacts.for_zone(zone) or '-'}")

outline: ModelEdge
horizons: H0 -> H1 -> H2 -> H3 -> H4 -> H5 -> H6
zones   : Z0 -> Z1 -> Z2 -> Z3 -> Z4 -> Z5
contacts:
  Z0: -
  Z1: -
  Z2: {'owc': 1951.3975311345757}
  Z3: -
  Z4: {'goc': 1968.9450766064617, 'fwl': 1974.5450766064619}
  Z5: -

5. Build The Stack Model¶

This is the main notebook workflow: Project.load -> grid_geometry -> geom.build -> grid.model.

In [5]:

geom = project.grid_geometry(cell=(50.0, 50.0), orient=0.0)
grid = geom.build(
    horizons,
    subzones=subzones,
    layering=layering,
    collapse_negative=True,
    outline=OUTLINE,
    min_thickness_m=0.0,
)
model = grid.model(
    props,
    contacts,
    fluid="oil",
    fvf=1.30,
    gas_fvf=0.005,
)

print("model:", model)
print("properties:", model.property_names())

geom = project.grid_geometry(cell=(50.0, 50.0), orient=0.0) grid = geom.build( horizons, subzones=subzones, layering=layering, collapse_negative=True, outline=OUTLINE, min_thickness_m=0.0, ) model = grid.model( props, contacts, fluid="oil", fvf=1.30, gas_fvf=0.005, ) print("model:", model) print("properties:", model.property_names())

model: Model(zoned=True, props=['PORO', 'SW', 'NTG'])
properties: ['PORO', 'SW', 'NTG']

6. Inspect The Bound Framework¶

This prints the model roles from code, not hidden synthetic defaults.

In [6]:

print("outline:", OUTLINE)
print("horizons:", " -> ".join(HORIZON_NAMES))
print("zones   :", " -> ".join(ZONE_NAMES))
print("subzones:")
for zone in ZONE_NAMES:
    print(f"  {zone}: {SUBZONE_BINDINGS.get(zone, []) or '-'}")
print("contacts:")
for zone in ZONE_NAMES:
    print(f"  {zone}: {contacts.for_zone(zone) or '-'}")

print("outline:", OUTLINE) print("horizons:", " -> ".join(HORIZON_NAMES)) print("zones :", " -> ".join(ZONE_NAMES)) print("subzones:") for zone in ZONE_NAMES: print(f" {zone}: {SUBZONE_BINDINGS.get(zone, []) or '-'}") print("contacts:") for zone in ZONE_NAMES: print(f" {zone}: {contacts.for_zone(zone) or '-'}")

outline: ModelEdge
horizons: H0 -> H1 -> H2 -> H3 -> H4 -> H5 -> H6
zones   : Z0 -> Z1 -> Z2 -> Z3 -> Z4 -> Z5
subzones:
  Z0: -
  Z1: -
  Z2: -
  Z3: ['H3b']
  Z4: -
  Z5: -
contacts:
  Z0: -
  Z1: -
  Z2: {'owc': 1951.3975311345757}
  Z3: -
  Z4: {'goc': 1968.9450766064617, 'fwl': 1974.5450766064619}
  Z5: -

7. Per-Zone Volumes¶

in_place_by_zone() reads volumes from the populated StaticModel, using each zone's own contacts.

In [7]:

by_zone = model.in_place_by_zone()
rows = by_zone["zones"]
total = by_zone["total"]

hdr = f'{"zone":6} {"GRV[Mm3]":>10} {"HCPV[m3]":>14} {"STOIIP[MSm3]":>13} {"2-contact":>10}'
print(hdr)
print("-" * len(hdr))
for row in rows:
    print(
        f'{row["zone"]:6} {row["grv_mcm"]:>10.2f} {row["hcpv_m3"]:>14.1f} '
        f'{row["stoiip_msm3"]:>13.4f} {str(row["two_contact"]):>10}'
    )
print("-" * len(hdr))
print(f'{"TOTAL":6} {total["grv_mcm"]:>10.2f} {total["hcpv_m3"]:>14.1f} {total["stoiip_msm3"]:>13.4f}')

by_zone = model.in_place_by_zone() rows = by_zone["zones"] total = by_zone["total"] hdr = f'{"zone":6} {"GRV[Mm3]":>10} {"HCPV[m3]":>14} {"STOIIP[MSm3]":>13} {"2-contact":>10}' print(hdr) print("-" * len(hdr)) for row in rows: print( f'{row["zone"]:6} {row["grv_mcm"]:>10.2f} {row["hcpv_m3"]:>14.1f} ' f'{row["stoiip_msm3"]:>13.4f} {str(row["two_contact"]):>10}' ) print("-" * len(hdr)) print(f'{"TOTAL":6} {total["grv_mcm"]:>10.2f} {total["hcpv_m3"]:>14.1f} {total["stoiip_msm3"]:>13.4f}')

zone     GRV[Mm3]       HCPV[m3]  STOIIP[MSm3]  2-contact
---------------------------------------------------------
Z0         247.87            0.0        0.0000      False
Z1         193.40            0.0        0.0000      False
Z2           1.11       149560.1        0.1150      False
Z3         142.45            0.0        0.0000      False
Z4           1.09       143630.9        0.1105       True
Z5         109.01            0.0        0.0000      False
---------------------------------------------------------
TOTAL      694.95       293191.0        0.2255

8. Synthetic Manifest Check¶

Synthetic-only verification lives here. For a real-data run this cell reports that the check is skipped.

In [8]:

if synthetic_manifest is None:
    print("Skipped: synthetic manifest checks apply only to petektools.synth_asset runs.")
else:
    assert OUTLINE == "ModelEdge"
    assert HORIZON_NAMES == synthetic_manifest["horizons"]
    assert ZONE_NAMES == synthetic_manifest["zones"]
    assert synthetic_manifest["surfaces_as_points"] is True

    det = {row["zone"]: row for row in model.in_place_by_zone()["zones"]}
    for zone, plan in synthetic_manifest["contact_plan"].items():
        row = det[zone]
        if plan["type"] == "single":
            assert row["stoiip_msm3"] > 0.0 and not row["two_contact"], (zone, row)
        elif plan["type"] == "two_contact":
            assert row["stoiip_msm3"] > 0.0 and row["two_contact"], (zone, row)
        else:
            assert row["stoiip_msm3"] == 0.0, (zone, row)

    print("Synthetic manifest checks passed.")

if synthetic_manifest is None: print("Skipped: synthetic manifest checks apply only to petektools.synth_asset runs.") else: assert OUTLINE == "ModelEdge" assert HORIZON_NAMES == synthetic_manifest["horizons"] assert ZONE_NAMES == synthetic_manifest["zones"] assert synthetic_manifest["surfaces_as_points"] is True det = {row["zone"]: row for row in model.in_place_by_zone()["zones"]} for zone, plan in synthetic_manifest["contact_plan"].items(): row = det[zone] if plan["type"] == "single": assert row["stoiip_msm3"] > 0.0 and not row["two_contact"], (zone, row) elif plan["type"] == "two_contact": assert row["stoiip_msm3"] > 0.0 and row["two_contact"], (zone, row) else: assert row["stoiip_msm3"] == 0.0, (zone, row) print("Synthetic manifest checks passed.")

Synthetic manifest checks passed.