pycsamt.ai.geology.layers#

Correlated stratigraphic interfaces and layered resistivity priors.

Layer interfaces are stored as depth below the grid datum and have shape (n_interface, nx) in 2-D or (n_interface, ny, nx) in 3-D. Resulting resistivity arrays follow the geology convention (nz, nx) or (nz, ny, nx).

Functions

generate_layered_geology(grid, layers, ...)

Generate a correlated 2-D section or 3-D layered electrical volume.

Classes

`ElectricalLayer`(name, resistivity_ohm_m[, ...])	Define the electrical distribution of one stratigraphic unit.
`LayeredGeology`(grid, layers, ...[, ...])	Immutable discretized layered electrical geology.

class pycsamt.ai.geology.layers.ElectricalLayer(name, resistivity_ohm_m, log10_std=0.0, heterogeneity=None, resistivity_bounds_ohm_m=None)[source]

Bases: object

Define the electrical distribution of one stratigraphic unit.

Parameters:

name (str) – Unique non-empty unit name.
resistivity_ohm_m (float) – Positive median resistivity in ohm metres.
log10_std (float, default=0.0) – Non-negative standard deviation of within-unit log10 resistivity.
heterogeneity (GaussianCorrelation or None, optional) – Spatial correlation model required when log10_std is positive.
resistivity_bounds_ohm_m (tuple of float or None, optional) – Positive inclusive lower/upper clipping bounds.

Examples

A homogeneous conductive cover:

>>> cover = ElectricalLayer("conductive cover", 10.0)
>>> cover.log10_resistivity
1.0

A heterogeneous basement:

>>> basement = ElectricalLayer(
...     "basement",
...     1000.0,
...     log10_std=0.2,
...     heterogeneity=GaussianCorrelation(1000, 200),
...     resistivity_bounds_ohm_m=(100, 5000),
... )

name: str

resistivity_ohm_m: float

log10_std: float = 0.0

heterogeneity: GaussianCorrelation | None = None

resistivity_bounds_ohm_m: tuple[float, float] | None = None

property log10_resistivity: float[source]

Return the median resistivity in log10 ohm metres.

Returns:: log10(resistivity_ohm_m).
Return type:: float

Examples

>>> ElectricalLayer("unit", 100).log10_resistivity
2.0

to_dict()[source]

Return a JSON-serializable unit definition.

Returns:: Electrical distribution and optional correlation state.
Return type:: dict

Examples

>>> ElectricalLayer("unit", 100).to_dict()["resistivity_ohm_m"]
100.0

classmethod from_dict(data)[source]

Restore a validated electrical layer.

Parameters:: data (mapping) – State returned by to_dict().
Returns:: Immutable unit definition.
Return type:: ElectricalLayer

Examples

>>> unit = ElectricalLayer("unit", 50)
>>> ElectricalLayer.from_dict(unit.to_dict()) == unit
True

class pycsamt.ai.geology.layers.LayeredGeology(grid, layers, interface_depth_m, layer_index, resistivity_ohm_m, seed, boundary, interface_policy, minimum_thickness_m, adjusted_interface_fraction=0.0, generation_config=<factory>)[source]

Bases: object

Immutable discretized layered electrical geology.

Parameters:

grid (GeologyGrid) – Target model grid.
layers (sequence of ElectricalLayer) – Units ordered shallowest to deepest.
interface_depth_m (ndarray) – Interface surfaces shaped (n_layer - 1, *horizontal_shape).
layer_index (ndarray of int) – Zero-based unit index for every model cell, shaped like grid.
resistivity_ohm_m (ndarray) – Positive cell resistivities shaped like grid.
seed (int) – Root generation seed.
boundary (str) – Recorded generation policies.
interface_policy (str) – Recorded generation policies.
minimum_thickness_m (float) – Enforced minimum vertical separation.
adjusted_interface_fraction (float) – Fraction of interface samples changed by projection.
generation_config (mapping, optional) – Finite JSON-compatible requested interface means, relief amplitudes, and correlation models needed to regenerate the surfaces.

Examples

>>> grid = GeologyGrid.regular_2d(nx=8, nz=6, dx_m=100, dz_m=50)
>>> layers = [
...     ElectricalLayer("cover", 10),
...     ElectricalLayer("basement", 1000),
... ]
>>> model = generate_layered_geology(grid, layers, [150], seed=1)
>>> model.resistivity_ohm_m.shape
(6, 8)
>>> model.n_layers
2

grid: GeologyGrid

layers: tuple[ElectricalLayer, ...]

interface_depth_m: ndarray

layer_index: ndarray

resistivity_ohm_m: ndarray

seed: int

boundary: str

interface_policy: str

minimum_thickness_m: float

adjusted_interface_fraction: float = 0.0

generation_config: Mapping[str, Any]

property n_layers: int[source]

Return the number of stratigraphic units.

Returns:: Number of electrical layers.
Return type:: int

Examples

>>> grid = GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1)
>>> model = generate_layered_geology(
...     grid,
...     [ElectricalLayer("a", 1), ElectricalLayer("b", 2)],
...     [1],
...     seed=0,
... )
>>> model.n_layers
2

property model_hash: str[source]

Return a platform-stable digest of model values and provenance.

Returns:: SHA-256 digest.
Return type:: str

Examples

>>> grid = GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1)
>>> model = generate_layered_geology(
...     grid,
...     [ElectricalLayer("a", 1), ElectricalLayer("b", 2)],
...     [1],
...     seed=0,
... )
>>> len(model.model_hash)
64

layer_mask(layer)[source]

Return a read-only Boolean mask for one unit.

Parameters:

layer (int or str) – Zero-based index or exact layer name.

Returns:

Mask shaped like the geological grid.

Return type:

ndarray of bool

Raises:

KeyError – If a named unit is unavailable.
IndexError – If a numeric index is outside the layer range.

Examples

>>> grid = GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1)
>>> model = generate_layered_geology(
...     grid,
...     [ElectricalLayer("a", 1), ElectricalLayer("b", 2)],
...     [1],
...     seed=0,
... )
>>> model.layer_mask("a").shape
(2, 2)

interface(index)[source]

Return one read-only interface-depth surface.

Parameters:: index (int) – Zero-based interface between layers index and index + 1.
Returns:: Shape (nx,) in 2-D or (ny, nx) in 3-D.
Return type:: ndarray

Examples

>>> grid = GeologyGrid.regular_2d(nx=3, nz=2, dx_m=1, dz_m=1)
>>> model = generate_layered_geology(
...     grid,
...     [ElectricalLayer("a", 1), ElectricalLayer("b", 2)],
...     [1],
...     seed=0,
... )
>>> model.interface(0).shape
(3,)

summary()[source]

Return compact JSON-compatible geological diagnostics.

Returns:: Shape, unit fractions, interface ranges, resistivity range, and interface-adjustment fraction.
Return type:: dict

Examples

>>> grid = GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1)
>>> model = generate_layered_geology(
...     grid,
...     [ElectricalLayer("a", 1), ElectricalLayer("b", 2)],
...     [1],
...     seed=0,
... )
>>> model.summary()["n_layers"]
2

provenance()[source]

Return generation provenance without cell arrays.

Returns:: Grid, layer definitions, seed, and interface policies.
Return type:: dict

Examples

>>> grid = GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1)
>>> model = generate_layered_geology(
...     grid,
...     [ElectricalLayer("a", 1), ElectricalLayer("b", 2)],
...     [1],
...     seed=0,
... )
>>> model.provenance()["seed"]
0

to_npz(path)[source]

Persist the layered model in a pickle-free compressed archive.

Parameters:: path (str or pathlib.Path) – Destination NPZ path.
Returns:: Requested destination.
Return type:: pathlib.Path

Examples

>>> from tempfile import TemporaryDirectory
>>> grid = GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1)
>>> model = generate_layered_geology(
...     grid,
...     [ElectricalLayer("a", 1), ElectricalLayer("b", 2)],
...     [1],
...     seed=0,
... )
>>> with TemporaryDirectory() as directory:
...     path = model.to_npz(Path(directory) / "layers.npz")
...     restored = LayeredGeology.from_npz(path)
>>> restored.model_hash == model.model_hash
True

classmethod from_npz(path)[source]

Load and validate a layered model without enabling pickle.

Parameters:: path (str or pathlib.Path) – Archive written by to_npz().
Returns:: Immutable restored model.
Return type:: LayeredGeology

Examples

>>> from tempfile import TemporaryDirectory
>>> model = generate_layered_geology(
...     GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1),
...     [ElectricalLayer("a", 1), ElectricalLayer("b", 2)],
...     [1],
...     seed=0,
... )
>>> with TemporaryDirectory() as directory:
...     path = model.to_npz(Path(directory) / "m.npz")
...     restored = LayeredGeology.from_npz(path)
>>> np.array_equal(restored.layer_index, model.layer_index)
True

pycsamt.ai.geology.layers.generate_layered_geology(grid, layers, mean_interface_depth_m, *, seed, interface_relief_std_m=0.0, interface_correlation=None, minimum_thickness_m=None, interface_policy='project', boundary='reflect')[source]

Generate a correlated 2-D section or 3-D layered electrical volume.

Parameters:

grid (GeologyGrid) – Regular target grid.
layers (sequence of ElectricalLayer) – Units ordered shallowest to deepest. At least one is required.
mean_interface_depth_m (sequence of float) – One mean depth for every adjacent layer pair.
seed (int) – Explicit root seed. Labeled child seeds independently drive interfaces and within-layer heterogeneity.
interface_relief_std_m (float or sequence, default=0.0) – Interface depth standard deviations in metres.
interface_correlation (GaussianCorrelation, sequence, or None, optional) – Horizontal correlation model(s). Required for every interface with positive relief. A single model is shared across interfaces.
minimum_thickness_m (float or None, optional) – Minimum separation between top boundary, interfaces, and bottom boundary. Default is one vertical cell.
interface_policy ({"project", "raise"}, default="project") – Project invalid/crossing surfaces into the feasible domain or reject the realization. Projection fraction is recorded on the result.
boundary ({"reflect", "periodic"}, default="reflect") – Boundary policy for all correlated fields.

Returns:

Immutable interfaces, unit indices, resistivity, and provenance.

Return type:

LayeredGeology

Raises:

ValueError – If layer/interface counts disagree, thickness is infeasible, relief lacks correlation, a 3-D correlation lacks length_y_m, or strict interface validation fails.

Examples

Laterally varying three-layer 2-D model:

>>> grid = GeologyGrid.regular_2d(nx=24, nz=16, dx_m=100, dz_m=50)
>>> units = [
...     ElectricalLayer("cover", 20),
...     ElectricalLayer("sediments", 100),
...     ElectricalLayer("basement", 1000),
... ]
>>> model = generate_layered_geology(
...     grid,
...     units,
...     [200, 500],
...     seed=5,
...     interface_relief_std_m=[30, 50],
...     interface_correlation=GaussianCorrelation(600, 100),
... )
>>> model.interface_depth_m.shape
(2, 24)
>>> set(np.unique(model.layer_index)) <= {0, 1, 2}
True