02 - Volumes, Bundles, And Segment Aggregation¶
This notebook keeps the same simple data setup as notebook 01, then focuses on volumetrics, view bundles, and the current segment pattern.
To use your own export, edit only the data-source cell and the role-binding literals. The final synthetic checks are isolated from the workflow.
1. Imports¶
Only the library imports live here.
In [1]:
Copied!
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]:
Copied!
# 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_vbhlw62u
3. Load Data And Inspect Inventory¶
Loading is separate from synthetic generation. Inspect the inventory before changing the role-binding cell.
In [3]:
Copied!
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]:
Copied!
# 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 Base Model¶
The remaining cells use this populated StaticModel to showcase volumes, bundle payloads, and segment aggregation.
In [5]:
Copied!
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=['SW', 'PORO', 'NTG']) properties: ['SW', 'PORO', 'NTG']
6. Per-Zone GRV / In-Place Table¶
The full breakdown includes gross rock volume, hydrocarbon pore volume, STOIIP, and the model rollup.
In [6]:
Copied!
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}'
print(hdr)
print("-" * len(hdr))
for row in rows:
print(f'{row["zone"]:6} {row["grv_mcm"]:>10.2f} {row["hcpv_m3"]:>14.1f} {row["stoiip_msm3"]:>13.4f}')
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}'
print(hdr)
print("-" * len(hdr))
for row in rows:
print(f'{row["zone"]:6} {row["grv_mcm"]:>10.2f} {row["hcpv_m3"]:>14.1f} {row["stoiip_msm3"]:>13.4f}')
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] ---------------------------------------------- Z0 247.87 0.0 0.0000 Z1 193.40 0.0 0.0000 Z2 1.11 149560.1 0.1150 Z3 142.45 0.0 0.0000 Z4 1.09 143630.9 0.1105 Z5 109.01 0.0 0.0000 ---------------------------------------------- TOTAL 694.95 293191.0 0.2255
7. Contact Scenario¶
Contacts.replace derives a new contact spec. The geometry and role bindings stay fixed; only the contact role value changes.
In [7]:
Copied!
oil_zone = next(zone for zone in ZONE_NAMES if (contacts.for_zone(zone) or {}).get("owc") is not None)
base_owc = contacts.for_zone(oil_zone)["owc"]
deep_owc = base_owc + 40.0
contacts_deep = contacts.replace(oil_zone, owc=deep_owc)
model_deep = grid.model(props, contacts_deep, fluid="oil", fvf=1.30, gas_fvf=0.005)
def zone_stoiip_msm3(static_model, zone_name):
return next(
row["stoiip_msm3"]
for row in static_model.in_place_by_zone()["zones"]
if row["zone"] == zone_name
)
base_v = zone_stoiip_msm3(model, oil_zone)
deep_v = zone_stoiip_msm3(model_deep, oil_zone)
print(f"zone {oil_zone}")
print(f" base OWC {base_owc:7.1f} m -> STOIIP {base_v:.4f} MSm3")
print(f" deep OWC {deep_owc:7.1f} m -> STOIIP {deep_v:.4f} MSm3")
print(f" uplift: +{deep_v - base_v:.4f} MSm3")
oil_zone = next(zone for zone in ZONE_NAMES if (contacts.for_zone(zone) or {}).get("owc") is not None)
base_owc = contacts.for_zone(oil_zone)["owc"]
deep_owc = base_owc + 40.0
contacts_deep = contacts.replace(oil_zone, owc=deep_owc)
model_deep = grid.model(props, contacts_deep, fluid="oil", fvf=1.30, gas_fvf=0.005)
def zone_stoiip_msm3(static_model, zone_name):
return next(
row["stoiip_msm3"]
for row in static_model.in_place_by_zone()["zones"]
if row["zone"] == zone_name
)
base_v = zone_stoiip_msm3(model, oil_zone)
deep_v = zone_stoiip_msm3(model_deep, oil_zone)
print(f"zone {oil_zone}")
print(f" base OWC {base_owc:7.1f} m -> STOIIP {base_v:.4f} MSm3")
print(f" deep OWC {deep_owc:7.1f} m -> STOIIP {deep_v:.4f} MSm3")
print(f" uplift: +{deep_v - base_v:.4f} MSm3")
zone Z2 base OWC 1951.4 m -> STOIIP 0.1150 MSm3 deep OWC 1991.4 m -> STOIIP 7.4402 MSm3 uplift: +7.3252 MSm3
8. View Bundles¶
The model emits self-contained JSON bundles for downstream viewers. This cell prints bundle shape only.
In [8]:
Copied!
volume_bundle = model.volume_bundle(property="PORO")
map_bundle = model.map_bundle(property="PORO")
print("volume_bundle keys:", sorted(volume_bundle.keys()))
print(" property :", volume_bundle["property"])
print(" cells :", volume_bundle["cell_count"])
print(" triangles :", volume_bundle["triangle_count"])
print(" encoding :", volume_bundle["encoding"])
print()
frame = map_bundle["frame"]
print("map_bundle keys:", sorted(map_bundle.keys()))
print(f' frame : {frame["ncol"]} x {frame["nrow"]}, {frame["spacing_x"]:.0f} m spacing')
print(" zone maps :", len(map_bundle["zone_averages"]))
volume_bundle = model.volume_bundle(property="PORO")
map_bundle = model.map_bundle(property="PORO")
print("volume_bundle keys:", sorted(volume_bundle.keys()))
print(" property :", volume_bundle["property"])
print(" cells :", volume_bundle["cell_count"])
print(" triangles :", volume_bundle["triangle_count"])
print(" encoding :", volume_bundle["encoding"])
print()
frame = map_bundle["frame"]
print("map_bundle keys:", sorted(map_bundle.keys()))
print(f' frame : {frame["ncol"]} x {frame["nrow"]}, {frame["spacing_x"]:.0f} m spacing')
print(" zone maps :", len(map_bundle["zone_averages"]))
volume_bundle keys: ['blocks', 'cell_count', 'encoding', 'inputs_ref', 'kind', 'property', 'schema_version', 'shell_cell_count', 'triangle_count', 'value_range', 'vertex_count', 'zone_names'] property : PORO cells : 76800 triangles : 35776 encoding : base64 map_bundle keys: ['contacts', 'frame', 'horizons', 'inputs_ref', 'k_slices', 'outline', 'schema_version', 'wells', 'zone_averages'] frame : 80 x 80, 50 m spacing zone maps : 6
9. Current Segment Pattern¶
For segmented fields, keep each segment as its own run. Build/run segment A, build/run segment B, and so on, then combine segment-level Uncertainty results with ps.aggregate. petekStatic should not create a zone-by-segment single-model rollup.
In [9]:
Copied!
# Synthetic has one segment, so this builds one segment-level model and
# aggregates one segment. A real multi-segment workflow repeats this block per
# segment, usually with a segment-specific outline, then passes all runs to
# ps.aggregate.
segment_horizons = ps.Horizons(
ps.hz(HORIZON_NAMES[0]),
ps.hz(HORIZON_NAMES[-1]),
ties=ps.TieSettings(method="convergent"),
gridding=ps.Gridding(collapse=True),
)
segment_grid = geom.build(
segment_horizons,
layering=ps.Layering(nk=6),
collapse_negative=True,
outline=OUTLINE,
min_thickness_m=0.0,
)
segment_contact = (
synthetic_manifest["contacts"]["fwl_z4"] + 60.0
if synthetic_manifest is not None
else 2455.0
)
segment_model = segment_grid.model(
props,
{"fwl": segment_contact},
fluid="oil",
fvf=1.30,
gas_fvf=0.005,
)
segment_runs = [
segment_model.uncertainty(
ps.Mc(contacts=4.0, porosity=0.01, n=64, seed=42)
)
]
field_total = ps.aggregate(segment_runs, correlation="independent")
print("segments:", len(segment_runs))
print(
"field STOIIP P90/P50/P10:",
f'{field_total["p90_msm3"]:.4f} / {field_total["p50_msm3"]:.4f} / {field_total["p10_msm3"]:.4f} MSm3',
)
# Synthetic has one segment, so this builds one segment-level model and
# aggregates one segment. A real multi-segment workflow repeats this block per
# segment, usually with a segment-specific outline, then passes all runs to
# ps.aggregate.
segment_horizons = ps.Horizons(
ps.hz(HORIZON_NAMES[0]),
ps.hz(HORIZON_NAMES[-1]),
ties=ps.TieSettings(method="convergent"),
gridding=ps.Gridding(collapse=True),
)
segment_grid = geom.build(
segment_horizons,
layering=ps.Layering(nk=6),
collapse_negative=True,
outline=OUTLINE,
min_thickness_m=0.0,
)
segment_contact = (
synthetic_manifest["contacts"]["fwl_z4"] + 60.0
if synthetic_manifest is not None
else 2455.0
)
segment_model = segment_grid.model(
props,
{"fwl": segment_contact},
fluid="oil",
fvf=1.30,
gas_fvf=0.005,
)
segment_runs = [
segment_model.uncertainty(
ps.Mc(contacts=4.0, porosity=0.01, n=64, seed=42)
)
]
field_total = ps.aggregate(segment_runs, correlation="independent")
print("segments:", len(segment_runs))
print(
"field STOIIP P90/P50/P10:",
f'{field_total["p90_msm3"]:.4f} / {field_total["p50_msm3"]:.4f} / {field_total["p10_msm3"]:.4f} MSm3',
)
segments: 1 field STOIIP P90/P50/P10: 64.9202 / 71.0416 / 77.1838 MSm3
10. Synthetic Manifest Check¶
Synthetic-only verification lives here. For a real-data run this cell reports that the check is skipped.
In [10]:
Copied!
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 set(CONTACTS_BY_ZONE) == {
zone for zone, plan in synthetic_manifest["contact_plan"].items()
if plan["contacts"] is not None
}
assert model.volume_bundle(property="PORO")["cell_count"] > 0
assert model.map_bundle(property="PORO")["zone_averages"]
print("Synthetic manifest and bundle 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 set(CONTACTS_BY_ZONE) == {
zone for zone, plan in synthetic_manifest["contact_plan"].items()
if plan["contacts"] is not None
}
assert model.volume_bundle(property="PORO")["cell_count"] > 0
assert model.map_bundle(property="PORO")["zone_averages"]
print("Synthetic manifest and bundle checks passed.")
Synthetic manifest and bundle checks passed.