2.17. pycsamt.interp#

Hydrogeological interpretation helpers, borehole calibration, uncertainty, fusion, and report/export routines, built on top of the geological domain knowledge in pycsamt.geology.

pycsamt.interp — geophysical-geological interpretation for EM methods.

This package converts 2-D EM resistivity models (from any inversion code — Occam2D, ModEM, AI-based, or generic arrays) into geological deliverables: pseudo-stratigraphic logs, borehole-calibrated models, and industry-standard exports.

It supersedes the geodrill module from pyCSAMT v1, extending its scope beyond OCCAM2D and CSAMT to cover MT, AMT, CSAMT, and CSEM.

RockDatabase, RockEntry, Layer, StratigraphicLog, Borehole, and Interval moved to pycsamt.geology — general geological domain knowledge with no EM dependency. They remain importable from here (pycsamt.interp.RockDatabase and so on) for backward compatibility; see pycsamt.geology for their canonical location and pycsamt.geology.lithology/.borehole for the submodules that used to live at pycsamt.interp.lithology/.borehole (removed; import from pycsamt.geology directly if you relied on that submodule path).

2.17.1. Quickstart#

>>> from pycsamt.interp import ResistivityModel, ModelCalibrator
>>> from pycsamt.interp import export, plot as iplot
>>>
>>> # Adapt any inversion result
>>> model = ResistivityModel.from_occam2d(result)
>>>
>>> # Optional: calibrate with borehole data
>>> from pycsamt.interp import Borehole
>>> bh = Borehole.from_csv("boreholes/Bo.csv", x=1050.0)
>>> cal = ModelCalibrator(ptol=0.10).fit(model, [bh])
>>>
>>> # Build stratigraphic logs
>>> logs = cal.stratigraphic_logs()
>>>
>>> # Export
>>> export.to_oasis_montaj_xyz(logs, "profile.xyz")
>>> export.to_las(logs[0], "S17.las")
>>> export.to_csv(logs, "all_stations.csv")
>>>
>>> # Plot
>>> fig = iplot.PlotStratigraphicLog(logs[0]).plot()
>>> fig = iplot.PlotFenceDiagram(logs).plot()
>>> fig = iplot.PlotCalibratedModel(
...     model, cal.calibrated_model(), cal.misfit_map()
... ).plot()

2.17.2. Package layout#

_base: ResistivityModel — unified model container; adapters for Occam2D, ModEM, AI outputs, and raw numpy arrays.
calibrate: ModelCalibrator — generalised NM construction (Kouadio et al. 2022), borehole-constrained or database-only; imports RockDatabase and Borehole from pycsamt.geology.
export: to_oasis_montaj_xyz(), to_las(), to_csv(), to_vtk() — industry-format writers.
plot: PlotStratigraphicLog, PlotFenceDiagram, PlotCalibratedModel — matplotlib-based visuals.

References

class pycsamt.interp.ResistivityModel(x_centers, z_centers, rho_2d, station_x=<factory>, station_names=<factory>, method='generic', rms=nan, metadata=<factory>)#

Bases: PyCSAMTObject, MetadataMixin

Unified 2-D resistivity model container.

All spatial arrays use metres; depth is positive downward.

Parameters:

x_centers (ndarray (n_x,)) – Horizontal cell-centre positions along the profile, metres.
z_centers (ndarray (n_z,)) – Depth cell-centre positions, metres (positive downward).
rho_2d (ndarray (n_z, n_x)) – $\log_{10}(\rho / \Omega\mathrm{m})$ on the cell grid.
station_x (ndarray (n_sta,)) – Station positions along the profile, metres.
station_names (list of str) – Station labels, same length as station_x.
method (str) – Source inversion method tag ('occam2d', 'modem', 'ai', 'generic').
rms (float) – Final RMS misfit from the inversion; nan if unknown.
metadata (dict) – Free-form provenance (coordinate reference, backend version, original file paths, and similar). Empty unless populated explicitly or via update_metadata().

x_centers: ndarray#

z_centers: ndarray#

rho_2d: ndarray#

station_x: ndarray#

station_names: list[str]#

method: str = 'generic'#

rms: float = nan#

metadata: dict[str, Any]#

property n_x: int#

property n_z: int#

property depth_max: float#

property profile_length: float#

column_nearest(x)#

Return the log10-rho column (n_z,) nearest to position x.

Parameters:: x (float)
Return type:: ndarray

station_column(station_name)#

Return the column nearest to the named station.

Parameters:: station_name (str)
Return type:: ndarray

clip_to_stations(*, pad_m=150.0)#

Return a copy restricted to the station-carrying core.

A model built through from_occam2d() spans Occam2D’s own wide (often multi-kilometre) boundary-padding columns on both sides of the stations – real, but rarely useful for display or for anything that treats x_centers as “the profile”. This drops everything outside [station_x.min() - pad_m, station_x.max() + pad_m].

Parameters:: pad_m (float, default 150.0) – Margin kept on each side of the outermost stations, metres.
Raises:: ValueError – If station_x is empty (nothing to clip around).
Return type:: ResistivityModel

Examples

>>> import numpy as np
>>> from pycsamt.interp import ResistivityModel
>>> x = np.linspace(-5000.0, 5000.0, 50)
>>> z = np.array([10.0, 20.0])
>>> rho = np.zeros((2, 50))
>>> m = ResistivityModel.from_array(rho, x, z, station_x=np.array([0.0, 100.0]))
>>> clipped = m.clip_to_stations(pad_m=150.0)
>>> clipped.x_centers.min() >= -150.0, clipped.x_centers.max() <= 250.0
(True, True)

clip_to_depth(max_depth_m)#

Return a copy restricted to z_centers <= max_depth_m.

Solvers that need a numerically distant lower boundary (Occam2D, ModEM, MARE2DEM) commonly extend their mesh to tens of kilometres depth even when the geological zone of interest is a couple of kilometres at most. Plotting the full mesh then squeezes that zone into a sliver of the figure. This keeps only the shallow cells actually worth displaying.

Parameters:: max_depth_m (float) – Deepest cell centre to keep, metres.
Raises:: ValueError – If no cell centre is at or above max_depth_m.
Return type:: ResistivityModel

Examples

>>> import numpy as np
>>> from pycsamt.interp import ResistivityModel
>>> x = np.array([0.0, 100.0])
>>> z = np.array([10.0, 500.0, 9000.0])
>>> rho = np.zeros((3, 2))
>>> m = ResistivityModel.from_array(rho, x, z)
>>> m.clip_to_depth(1000.0).z_centers
array([ 10., 500.])

classmethod from_occam2d(result)#

Build from an InversionResult.

Parameters:: result (InversionResult) – A loaded Occam2D post-inversion result.
Raises:: ValueError – If result does not contain a valid rho_2d grid.
Return type:: ResistivityModel

Notes

mesh.x_nodes is in the mesh’s own coordinate frame, which starts at the outer edge of from_data()’s geometrically-expanding horizontal padding – not at the first station. result.data.offsets (station_x below), by contrast, is real station chainage. Left uncorrected, every station’s nearest x_centers column collapses onto the same outermost padding cell (column_nearest/station_column silently return the same, meaningless column for every station). Padding is built identically on both sides (from_data’s left_pad/right_pad share one pad list), so the station-carrying core is exactly centred in the mesh; that symmetry is what lets the shift below be recovered from x_nodes and data.offsets alone, without assuming a specific padding cell count.

classmethod from_array(rho_2d, x_centers, z_centers, *, station_x=None, station_names=None, method='generic', rms=nan)#

Build directly from numpy arrays.

Parameters:

rho_2d (ndarray (n_z, n_x)) – $\log_{10}(\rho)$.
x_centers (ndarray) – Cell-centre coordinates, metres.
z_centers (ndarray) – Cell-centre coordinates, metres.
station_x (ndarray, optional) – Station positions. Defaults to x_centers.
station_names (list of str, optional) – Station labels.
method (str) – Metadata.
rms (float) – Metadata.

Return type:

ResistivityModel

classmethod from_any(model, *, station_x=None, station_names=None, unit='m')#

Build from any supported 2-D inversion result.

A single entry point on top of pycsamt.topo.section.extract_grid() (the same “accept anything” adapter plot_topo_section() uses for topography draping), so callers that just want a ResistivityModel — for pycsamt.interp.export’s Surfer writers, calibration, classification, or anything else in this package — do not need to know which of Occam2D, ModEM, the backend-neutral InversionResult, or an AI agent result actually produced it.

Parameters:

model (object) – Already a ResistivityModel (returned as-is), a raw (x_centers, z_centers, rho_2d) triple, a backend-neutral InversionResult, a native Occam2D or 2-D ModEM InversionResult, or an AI agent result exposing pred_rho.
station_x (ndarray, optional) – Overrides the station positions carried by model.
station_names (list of str, optional) – Overrides the station labels carried by model.
unit ({"m", "km"}) – Unit of model’s own coordinates when model is a raw array triple. Ignored for input forms that carry their own unit; the returned model is always in metres.

Raises:

TypeError – If model is not a recognised input form.
NotImplementedError – For a native MARE2DEM result (unstructured triangular mesh — regrid onto a regular grid first).
ValueError – For a native 3-D ModEM result (extract a 2-D cut first).

Return type:

ResistivityModel

Examples

>>> import numpy as np
>>> from pycsamt.interp import ResistivityModel
>>> x = np.array([0.0, 100.0, 200.0])
>>> z = np.array([10.0, 50.0])
>>> rho_log10 = np.array([[2.0, 2.1, 2.2], [2.5, 2.6, 2.7]])
>>> model = ResistivityModel.from_any((x, z, rho_log10))
>>> model.method, model.n_x, model.n_z
('array', 3, 2)

class pycsamt.interp.ModelCalibrator(ptol=0.1, *, max_borehole_distance=500.0, db=None, verbose=True)#

Bases: PyCSAMTObject

Constrain a 2-D EM resistivity model with borehole TRES values.

Parameters:

ptol (float, default 0.10) – Soft-error threshold — fractional tolerance between a CRM resistivity value and a TRES value for accepting a replace. The paper uses 10 % (ptol = 0.10).
max_borehole_distance (float, default 500.0) – Maximum profile-distance (metres) from a station to the nearest borehole. Stations farther than this threshold are calibrated using the rock database only (no TRES constraint).
db (RockDatabase, optional) – Rock physics database for autolayer assignment. Defaults to default().
verbose (bool, default True)

fit(model, boreholes=None)#

Calibrate model against boreholes.

Parameters:

model (ResistivityModel) – The Calculated Resistivity Model (CRM) from any EM inversion.
boreholes (list of Borehole, optional) – Ground-truth borehole logs. If None or empty, the calibration falls back to rock-database autolayer assignment for every cell.

Return type:

self

calibrated_model()#

Return a ResistivityModel containing the NM.

Raises:: RuntimeError – If fit() has not been called.
Return type:: ResistivityModel

misfit_map()#

Return the per-column misfit G (%), shape (n_z, n_x).

Values range 0 – 100 %. Low values indicate good CRM–TRES agreement; high values flag where the inversion needed the most correction.

Return type:: ndarray

stratigraphic_logs(*, db=None, model='nm', merge_tolerance=0.2)#

Build a StratigraphicLog for every station.

Parameters:

db (RockDatabase, optional) – Override the classifier database.
model ({'nm', 'crm'}) – Which model to classify: the calibrated NM or the original CRM.
merge_tolerance (float) – Log-space merging tolerance passed to from_column().

Return type:

list of StratigraphicLog

class pycsamt.interp.AquiferZone(station, x, top, bottom, mean_rho_ohm_m, confidence, zone_type='aquifer', metadata=<factory>)#

Bases: PyCSAMTObject, MetadataMixin

Contiguous aquifer-favourable interval at one station/profile column.

Parameters:

station (str)
x (float)
top (float)
bottom (float)
mean_rho_ohm_m (float)
confidence (float)
zone_type (str)
metadata (dict[str, Any])

station: str#

x: float#

top: float#

bottom: float#

mean_rho_ohm_m: float#

confidence: float#

zone_type: str = 'aquifer'#

metadata: dict[str, Any]#

property thickness: float#: Interval thickness in metres.

class pycsamt.interp.HydroGeophysicalModel(resistivity_model, unit_map, confidence, zones=<factory>, logs=<factory>, metadata=<factory>)#

Bases: PyCSAMTObject, MetadataMixin

Hydrogeological interpretation product for one resistivity model.

Parameters:

resistivity_model (ResistivityModel)
unit_map (ndarray)
confidence (ndarray)
zones (list[AquiferZone])
logs (list[StratigraphicLog])
metadata (dict[str, Any])

resistivity_model: ResistivityModel#

unit_map: ndarray#

confidence: ndarray#

zones: list[AquiferZone]#

logs: list[StratigraphicLog]#

metadata: dict[str, Any]#

aquifer_zones(*, min_confidence=0.0)#

Return aquifer-favourable zones above a confidence threshold.

Parameters:: min_confidence (float)
Return type:: list[AquiferZone]

station_summary()#

Return compact per-station hydrogeological summaries.

Return type:: list[dict[str, Any]]

to_csv(path)#

Write cell-level hydro units and aquifer zones to CSV.

Parameters:: path (str | Path)
Return type:: Path

zones_to_csv(path)#

Write interpreted aquifer/fracture target intervals to CSV.

Parameters:: path (str | Path)
Return type:: Path

class pycsamt.interp.HydroInterpreter(*, context='', db=None, water_table_depth=None, aquifer_range=(5.0, 300.0), fracture_range=(50.0, 800.0), clay_max=20.0, saline_max=3.0, basement_min=1000.0, min_zone_thickness=1.0, calibration_ptol=0.1, max_borehole_distance=500.0)#

Bases: PyCSAMTObject

Rule-based hydrogeophysical interpreter for EM resistivity sections.

The default thresholds are conservative and intended as a first-pass screening model. Local calibration can be introduced by passing boreholes, a custom RockDatabase, or explicit threshold overrides.

Parameters:

context (str)
db (RockDatabase | None)
water_table_depth (float | None)
aquifer_range (tuple[float, float])
fracture_range (tuple[float, float])
clay_max (float)
saline_max (float)
basement_min (float)
min_zone_thickness (float)
calibration_ptol (float)
max_borehole_distance (float)

fit(model, *, boreholes=None)#

Interpret a resistivity model or inversion result.

Parameters:

model (Any)
boreholes (Sequence[Borehole] | None)

Return type:

HydroGeophysicalModel

aquifer_zones(*, min_confidence=0.0)#

Return aquifer zones from the last fitted model.

Parameters:: min_confidence (float)
Return type:: list[AquiferZone]

class pycsamt.interp.HydroUnit(x, z, rho_ohm_m, rho_log10, unit, lithology, confidence=1.0, metadata=<factory>)#

Bases: PyCSAMTObject, MetadataMixin

One classified hydrogeophysical cell.

Parameters:

x (float)
z (float)
rho_ohm_m (float)
rho_log10 (float)
unit (str)
lithology (str)
confidence (float)
metadata (dict[str, Any])

x: float#

z: float#

rho_ohm_m: float#

rho_log10: float#

unit: str#

lithology: str#

confidence: float = 1.0#

metadata: dict[str, Any]#

class pycsamt.interp.PetrophysicalConfig(petro=<factory>, rho_w=20.0, porosity_prior=0.25, Sw_water_table_threshold=0.85, d50_m=0.00025, kozeny_C=180.0, kozeny_tortuosity=0.5, fracture_depth_m=None, fracture_rho_matrix=5000.0, specific_storage=0.0001, min_wt_search_depth=0.5)#

Bases: PyCSAMTObject

All petrophysical and hydraulic parameters needed by EMHydroModel.

Parameters:

petro (ArchieModel or WaxmanSmitsModel) – Petrophysical model for ρ ↔ (φ, Sw) transforms. Default is ArchieModel(m=1.8, n=2.0) — appropriate for clean sands targeted by shallow TDEM and AMT surveys.
rho_w (float) – Pore-water resistivity (Ω·m). Default 20 Ω·m (EC ≈ 0.5 mS/cm, potable fresh water at 25 °C; range: 0.2 Ω·m seawater – 100+ Ω·m pristine groundwater). Constrained by water-quality measurements or EC logs via pycsamt.interp.constraints.
porosity_prior (float) – Reference porosity used in the vadose zone and as an upper-clip for unrealistic Archie-inferred porosities (default 0.25).
Sw_water_table_threshold (float) – Saturation value that marks the water table in the Archie inverse (default 0.85; lower it for coarser-grained aquifers).
d50_m (float) – Median grain size (m) for Kozeny-Carman K (default 2.5 × 10⁻⁴ m, fine sand). Increase to ~1 × 10⁻³ m for coarse sand/gravel.
kozeny_C (float) – Kozeny constant × shape factor (default 180 for spheres).
kozeny_tortuosity (float) – Tortuosity correction in Kozeny-Carman (default 0.5).
fracture_depth_m (float or None) – Depth (m) below which fractured_zone_K() replaces Kozeny-Carman. Use for AMT/MT basement targets. None disables fracture K everywhere (default: None).
fracture_rho_matrix (float) – Background resistivity of intact rock for the cubic-law fracture K (Ω·m; default 5 000).
specific_storage (float) – Specific storage Ss (m⁻¹) for confined storativity (default 10⁻⁴).
min_wt_search_depth (float) – Minimum depth (m) for water-table detection; skips near-surface noise (default 0.5 m).

petro: ArchieModel | WaxmanSmitsModel#

rho_w: float = 20.0#

porosity_prior: float = 0.25#

Sw_water_table_threshold: float = 0.85#

d50_m: float = 0.00025#

kozeny_C: float = 180.0#

kozeny_tortuosity: float = 0.5#

fracture_depth_m: float | None = None#

fracture_rho_matrix: float = 5000.0#

specific_storage: float = 0.0001#

min_wt_search_depth: float = 0.5#

class pycsamt.interp.EMHydroModel(resistivity_model, config=None, *, petro=None, rho_w=None, porosity_prior=None, method_tag='')#

Bases: PyCSAMTObject

Quantitative hydrogeological model from an EM resistivity section.

Wires the petrophysical transforms in pycsamt.interp.petrophysics onto a full 2-D ResistivityModel to produce spatially continuous porosity, saturation, hydraulic-conductivity, and water-table maps.

The computation follows a two-pass strategy to break the Archie chicken-and-egg (φ needs Sw, Sw needs φ):

Water-table pass — scan each column with the Archie-inverse Sw estimator to find the saturation front (Sw ≥ threshold).
Cell-property pass — for cells below the water table assume full saturation (Sw = 1) and solve for φ; for cells above use porosity_prior and solve for Sw.

Parameters:

resistivity_model (ResistivityModel) – Source inversion model.
config (PetrophysicalConfig, optional) – Full parameter bundle. If None, default parameters are used. Individual keyword arguments below override specific fields.
petro (ArchieModel or WaxmanSmitsModel, optional) – Petrophysical model.
rho_w (float, optional) – Pore-water resistivity (Ω·m).
porosity_prior (float, optional) – Prior porosity.
method_tag (str, optional) – EM method label ('TDEM', 'AMT', 'MT', 'EMAP').

Examples

>>> cfg = PetrophysicalConfig(
...     rho_w=0.03, porosity_prior=0.20, fracture_depth_m=300.0
... )
>>> result = EMHydroModel(rm, cfg, method_tag="AMT").fit()
>>> result.water_table  # array of depths, one per x-column
>>> result.transmissivity  # T profile (m²/s)

fit()#

Run all petrophysical transforms and return EMHydroResult.

Return type:: EMHydroResult

class pycsamt.interp.EMHydroResult(resistivity_model, config, porosity, saturation, hydraulic_K, water_table, transmissivity, storativity_confined, storativity_unconfined, dar_zarrouk_TR, dar_zarrouk_S, tds, method_tag='', metadata=<factory>)#

Bases: PyCSAMTObject

Quantitative hydrogeological output of one EM resistivity section.

All 2-D arrays have shape (n_z, n_x) matching the source ResistivityModel; 1-D station arrays have shape (n_x,).

Variables:

resistivity_model (ResistivityModel)
config (PetrophysicalConfig)
porosity (ndarray (n_z, n_x)) – Effective porosity per cell (dimensionless, 0–1). Below the water table: Archie inverse at Sw = 1. Above the water table: config.porosity_prior.
saturation (ndarray (n_z, n_x)) – Water saturation Sw per cell (0–1). Below the water table: 1.0. Above the water table: Archie inverse at φ = porosity_prior.
hydraulic_K (ndarray (n_z, n_x)) – Hydraulic conductivity K (m/s). Porous media (z < fracture_depth_m): Kozeny-Carman from φ. Fractured basement (z ≥ fracture_depth_m): cubic-law from ρ contrast.
water_table (ndarray (n_x,)) – Estimated water-table depth (m, positive downward). nan where detection failed (resistivity column shows no saturation transition above the Sw threshold).
transmissivity (ndarray (n_x,)) – Aquifer transmissivity T = ∫K dz over the saturated zone (m²/s).
storativity_confined (ndarray (n_x,)) – Confined storativity S = Ss × b_sat (dimensionless).
storativity_unconfined (ndarray (n_x,)) – Unconfined storativity ≈ mean porosity in saturated zone (≈ specific yield).
dar_zarrouk_TR (ndarray (n_x,)) – Transverse resistance TR = Σ ρᵢ hᵢ (Ω·m²) over all cells.
dar_zarrouk_S (ndarray (n_x,)) – Longitudinal conductance S = Σ hᵢ/ρᵢ (siemens) over all cells.
tds (float) – Estimated total dissolved solids (mg/L) from config.rho_w.
method_tag (str) – Source EM method label (e.g. 'TDEM', 'AMT').
metadata (dict) – Provenance and parameter snapshot.

Parameters:

resistivity_model (ResistivityModel)
config (PetrophysicalConfig)
porosity (ndarray)
saturation (ndarray)
hydraulic_K (ndarray)
water_table (ndarray)
transmissivity (ndarray)
storativity_confined (ndarray)
storativity_unconfined (ndarray)
dar_zarrouk_TR (ndarray)
dar_zarrouk_S (ndarray)
tds (float)
method_tag (str)
metadata (dict)

resistivity_model: ResistivityModel#

config: PetrophysicalConfig#

porosity: ndarray#

saturation: ndarray#

hydraulic_K: ndarray#

water_table: ndarray#

transmissivity: ndarray#

storativity_confined: ndarray#

storativity_unconfined: ndarray#

dar_zarrouk_TR: ndarray#

dar_zarrouk_S: ndarray#

tds: float#

method_tag: str = ''#

metadata: dict#

station_report()#

Return one dict per model column with key hydro indicators.

Return type:: list[dict]

to_csv(path)#

Write cell-level hydro parameters to a flat CSV file.

Parameters:: path (str | Path)
Return type:: Path

station_report_csv(path)#

Write per-station hydro summary to CSV.

Parameters:: path (str | Path)
Return type:: Path

to_dataframe()#

Return station summary as a pandas.DataFrame (requires pandas).

Return type:: Any

class pycsamt.interp.ConstrainedCalibrator(constraints, *, calibrate_rho_w=True, calibrate_m=False, calibrate_phi_prior=False, rho_w_bounds=(0.003, 2.0), m_bounds=(1.2, 2.8), phi_bounds=(0.03, 0.6), n_restarts=1, verbose=False)#

Bases: PyCSAMTObject

Calibrate petrophysical parameters to match hydrogeological field data.

Uses scipy.optimize.minimize (L-BFGS-B) to find the parameter set that minimises the weighted sum of squared residuals between the EM-derived hydrogeological quantities and the supplied observations.

The optimisation is deterministic (single run from the current parameter values). If convergence is poor, set n_restarts > 1 to try multiple random starting points within the parameter bounds.

Parameters:

constraints (sequence of constraint objects) – Any mix of WaterLevelConstraint, PumpingTestConstraint, SlugTestConstraint, ECConstraint.
calibrate_rho_w (bool) – Fit pore-water resistivity ρ_w (default True). Almost always needed.
calibrate_m (bool) – Fit Archie cementation exponent m (default False).
calibrate_phi_prior (bool) – Fit prior porosity (default False).
rho_w_bounds ((float, float)) – Search bounds for ρ_w in Ω·m (default 0.003 – 2.0).
m_bounds ((float, float)) – Search bounds for Archie m (default 1.2 – 2.8).
phi_bounds ((float, float)) – Search bounds for porosity prior (default 0.03 – 0.60).
n_restarts (int) – Number of optimisation restarts with random initialisations (default 1).
verbose (bool) – Print iteration count and final misfit (default False).

:param Attributes (set after fit()): :param ———————————–: :param calibrated_config_: The fitted configuration. :type calibrated_config_: PetrophysicalConfig :param misfit_history_: Final misfit value per restart. :type misfit_history_: list of float :param opt_result_: Raw scipy result from the best restart. :type opt_result_: OptimizeResult or None

calibrated_config_: PetrophysicalConfig | None#

misfit_history_: list[float]#

opt_result_: Any#

fit(em_hydro_model)#

Calibrate and return the best EMHydroResult.

Parameters:: em_hydro_model (EMHydroModel) – The model to calibrate. Its config provides the starting point for the optimisation.
Returns:: Model result evaluated at the calibrated parameters.
Return type:: EMHydroResult

constraint_residuals(result)#

Return per-constraint residuals for the given result (for diagnostics).

Parameters:: result (EMHydroResult)
Return type:: list[dict]

class pycsamt.interp.WaterLevelConstraint(x, depth_m, uncertainty_m=1.0, station='')#

Bases: PyCSAMTObject

Piezometer or well water-level measurement.

The calibrator matches the EM-derived water-table depth at the nearest model column to the observed depth.

Parameters:

x (float) – Profile position of the measurement (m, same reference as the model).
depth_m (float) – Observed water-table depth below surface (m, positive downward).
uncertainty_m (float) – Measurement + representativeness uncertainty (m; default 1.0 m).
station (str) – Optional label for reporting.

x: float#

depth_m: float#

uncertainty_m: float = 1.0#

station: str = ''#

class pycsamt.interp.PumpingTestConstraint(x, T_m2s, S=None, uncertainty_factor=3.0, station='')#

Bases: PyCSAMTObject

Transmissivity (and optionally storativity) from a pumping test.

Misfit is computed in log₁₀ space because T spans orders of magnitude.

Parameters:

x (float) – Profile position of the pumping well (m).
T_m2s (float) – Observed transmissivity (m²/s).
S (float, optional) – Observed storativity (dimensionless). Currently used for reporting only; storativity calibration will be added in a future release.
uncertainty_factor (float) – Multiplicative uncertainty on T (default 3 → ±0.5 in log₁₀).
station (str) – Optional label.

x: float#

T_m2s: float#

S: float | None = None#

uncertainty_factor: float = 3.0#

station: str = ''#

class pycsamt.interp.SlugTestConstraint(x, K_ms, depth_m, uncertainty_factor=5.0, station='')#

Bases: PyCSAMTObject

Hydraulic conductivity K from a slug test or bail test.

Matches the model K at the nearest cell to (x, depth_m).

Parameters:

x (float) – Profile position (m).
K_ms (float) – Observed K (m/s).
depth_m (float) – Depth of the tested interval mid-point (m).
uncertainty_factor (float) – Multiplicative uncertainty (default 5 → ~0.7 in log₁₀).
station (str) – Optional label.

x: float#

K_ms: float#

depth_m: float#

uncertainty_factor: float = 5.0#

station: str = ''#

class pycsamt.interp.ECConstraint(x, ec_mscm, uncertainty_mscm=2.0, station='')#

Bases: PyCSAMTObject

Electrical conductivity of formation water (e.g., from an EC log or sample).

Directly constrains rho_w via the EC ↔ ρ_w conversion.

Parameters:

x (float) – Profile position (m).
ec_mscm (float) – Observed pore-water EC (mS/cm; ≈ 1 mS/cm for fresh water at 25 °C).
uncertainty_mscm (float) – Measurement uncertainty (mS/cm; default 2.0).
station (str) – Optional label.

x: float#

ec_mscm: float#

uncertainty_mscm: float = 2.0#

station: str = ''#

class pycsamt.interp.TimeLapseEM(surveys, *, times=None, labels=None)#

Bases: PyCSAMTObject

Time-lapse EM analysis for hydrogeological change detection.

Parameters:

surveys (sequence of ResistivityModel) – Time-ordered EM inversion results. All must share the same grid (checked on construction via assert_compatible_grids()).
times (sequence of float, optional) – Time stamps for each survey (any consistent unit: days, months, etc.). Used for labelling only — no numerical time-derivative is computed.
labels (sequence of str, optional) – Human-readable survey labels (e.g. ['dry2022', 'wet2023']).

Variables:

n_surveys (int)
n_z (n_x,)

property n_surveys: int#

property n_x: int#

property n_z: int#

resistivity_change(baseline_idx=0)#

Resistivity change relative to the baseline survey.

\[\Delta\log_{10}\rho_i = \log_{10}\rho_i - \log_{10}\rho_\text{baseline}\]

Parameters:: baseline_idx (int) – Index of the baseline survey (default 0).
Returns:: One array per non-baseline survey, in time order. Positive values indicate resistivity increase (drying/desaturation). Negative values indicate resistivity decrease (wetting/salinisation).
Return type:: list of ndarray (n_z, n_x)

saturation_change(petro, *, phi=0.25, rho_w=0.025, baseline_idx=0)#

Water-saturation change Δ*Sw* via Archie/WS inverse.

\[\Delta S_{w,i} = S_w(\rho_i) - S_w(\rho_\text{baseline})\]

Parameters:

petro (ArchieModel or WaxmanSmitsModel) – Petrophysical model for ρ → Sw inversion.
phi (float or ndarray) – Porosity used in the inversion. Scalar for a uniform profile; 2-D array (n_z, n_x) for a spatially varying prior.
rho_w (float) – Pore-water resistivity (Ω·m).
baseline_idx (int) – Index of the baseline survey (default 0).

Returns:

ΔSw arrays, one per non-baseline survey. Positive = wetting (saturation increase). Negative = drying (saturation decrease).

Return type:

list of ndarray (n_z, n_x)

water_content_change(petro, *, phi=0.25, rho_w=0.025, baseline_idx=0)#

Volumetric water-content change Δθ = φ · ΔSw.

Parameters:

petro (ArchieModel | WaxmanSmitsModel) – Same as saturation_change().
phi (float | ndarray) – Same as saturation_change().
rho_w (float) – Same as saturation_change().
baseline_idx (int) – Same as saturation_change().

Returns:

Δθ arrays (dimensionless, range ≈ −φ to +φ).

Return type:

list of ndarray (n_z, n_x)

water_table_displacement(petro, *, rho_w=0.025, Sw_threshold=0.85, min_depth=0.5, baseline_idx=0)#

Water-table depth change relative to the baseline survey.

\[\Delta z_{wt} = z_{wt}(t_i) - z_{wt}(t_\text{baseline})\]

Positive values → water table dropped (deeper, e.g. dry season). Negative values → water table rose (shallower, e.g. recharge).

nan is returned where the water table cannot be detected in either the baseline or the comparison survey.

Parameters:

petro (ArchieModel or WaxmanSmitsModel)
rho_w (float)
Sw_threshold (float) – Saturation level that defines the water table (default 0.85).
min_depth (float) – Minimum search depth in metres (default 0.5).
baseline_idx (int)

Returns:

Water-table displacement in metres, one row per non-baseline survey.

Return type:

ndarray (n_surveys−1, n_x)

water_table_map(petro, *, rho_w=0.025, Sw_threshold=0.85, min_depth=0.5)#

Water-table depth (m) for every survey and every column.

Returns:

nan where the water table could not be detected.

Return type:

ndarray (n_surveys, n_x)

Parameters:

petro (ArchieModel | WaxmanSmitsModel)
rho_w (float)
Sw_threshold (float)
min_depth (float)

resistivity_stats(baseline_idx=0)#

Summary statistics of resistivity change across all surveys.

Returns:: mean_delta, std_delta, max_increase, max_decrease — each an ndarray (n_z, n_x).
Return type:: dict with keys
Parameters:: baseline_idx (int)

pycsamt.interp.assert_compatible_grids(surveys, *, rtol=0.0001)#

Raise ValueError if surveys have incompatible grids.

Parameters:

surveys (sequence of ResistivityModel)
rtol (float) – Relative tolerance for coordinate comparison (default 1e-4).

Return type:

None

class pycsamt.interp.MultiMethodEMModel(primary, secondary, *, primary_max_depth=None, secondary_min_depth=None, blend='linear', blend_overlap=None, z_grid=None, sigmoid_k=0.02)#

Bases: PyCSAMTObject

Fuse two EM resistivity models onto a single depth grid.

Parameters:

primary (ResistivityModel) – The model trusted at shallow depths (e.g. TDEM, EMAP).
secondary (ResistivityModel) – The model trusted at deeper depths (e.g. AMT, MT).
primary_max_depth (float, optional) – Override the primary model’s maximum contributing depth (m). Defaults to primary.z_centers[-1].
secondary_min_depth (float, optional) – Override the secondary model’s minimum contributing depth (m). Defaults to secondary.z_centers[0].
blend (str) –
Blend strategy in the overlap zone:

'linear' (default)
Linear weight ramp from primary (top of overlap) to secondary (bottom of overlap).

'sigmoid'
Smooth S-curve parameterised by sigmoid_k. Avoids the slope discontinuity at the ends of the transition zone.

'rms_weighted'
Constant weights throughout: primary weight = rms_secondary / (rms_primary + rms_secondary). Requires both models to carry a valid rms value. Falls back to 'linear' if either RMS is nan.
blend_overlap (float, optional) – Restrict the blend transition to a window of this width (m) centred on the mid-point of the natural overlap zone. None uses the full overlap.
z_grid (ndarray, optional) – Explicit output depth-cell centres (m). Overrides the automatic union grid. Both models are interpolated onto this grid.
sigmoid_k (float) – Shape parameter for the sigmoid blend (m⁻¹; default 0.02, giving a smooth ~100 m transition for a 500 m overlap zone).

Variables:

diagnostics (FusionDiagnostics or None) – Set after merge() is called.

Examples

>>> fused_model = MultiMethodEMModel(
...     tdem_model, amt_model, blend="sigmoid", sigmoid_k=0.03
... ).merge()
>>> fused_model.method
'TDEM+AMT'

diagnostics_: FusionDiagnostics | None#

merge()#

Produce the fused ResistivityModel.

Returns:: Unified model on the output depth grid. method is set to '<primary_method>+<secondary_method>' (e.g. 'TDEM+AMT').
Return type:: ResistivityModel

class pycsamt.interp.FusionDiagnostics(z_overlap_start, z_overlap_end, has_overlap, blend_mode, primary_method, secondary_method, primary_rms, secondary_rms, n_z_fused, blend_weights)#

Bases: PyCSAMTObject

Metadata about a completed fusion operation.

Variables:

z_overlap_start (float) – Top of the depth zone where both methods contribute (m).
z_overlap_end (float) – Bottom of the overlap zone (m).
has_overlap (bool) – False if the two depth ranges are disjoint (simple concatenation).
blend_mode (str)
primary_method (str)
secondary_method (str)
primary_rms (float)
secondary_rms (float)
n_z_fused (int) – Total number of depth cells in the fused model.
blend_weights (ndarray (n_z,)) – Primary-model weight at each depth cell (1 = all primary, 0 = all secondary).

Parameters:

z_overlap_start (float)
z_overlap_end (float)
has_overlap (bool)
blend_mode (str)
primary_method (str)
secondary_method (str)
primary_rms (float)
secondary_rms (float)
n_z_fused (int)
blend_weights (ndarray)

z_overlap_start: float#

z_overlap_end: float#

has_overlap: bool#

blend_mode: str#

primary_method: str#

secondary_method: str#

primary_rms: float#

secondary_rms: float#

n_z_fused: int#

blend_weights: ndarray#

class pycsamt.interp.UncertaintyBounds(rho_w_range=None, m_range=None, n_range=None, phi_prior_range=None, dist='uniform')#

Bases: PyCSAMTObject

Prior distribution specification for Monte Carlo sampling.

Each *_range parameter activates that parameter as a free variable. Un-specified parameters are held fixed at the central PetrophysicalConfig values.

Parameters:

rho_w_range ((float, float), optional) – Pore-water resistivity range (Ω·m). uniform: (low, high); normal: (mean, std). Typical fresh-water range: (5, 100); brackish: (0.5, 20).
m_range ((float, float), optional) – Archie cementation exponent. Typical: (1.3, 2.5) for mixed lithology.
n_range ((float, float), optional) – Archie saturation exponent. Typical: (1.8, 2.5) for clean sand.
phi_prior_range ((float, float), optional) – Porosity prior. Typical: (0.10, 0.45).
dist (str) – Sampling distribution: 'uniform' (default, non-informative) or 'normal' (Gaussian — range is interpreted as (mean, std)).

Notes

At least one *_range must be set. rho_w_range alone covers the dominant source of uncertainty in most shallow EM surveys.

rho_w_range: tuple[float, float] | None = None#

m_range: tuple[float, float] | None = None#

n_range: tuple[float, float] | None = None#

phi_prior_range: tuple[float, float] | None = None#

dist: str = 'uniform'#

property n_free: int#: Number of free (uncertain) parameters.

property free_names: list[str]#: Names of the free parameters.

sample(cfg, n, rng)#

Draw n parameter samples and return a list of configs.

Parameters:

cfg (PetrophysicalConfig) – Central (best-estimate) configuration — used as fallback for fixed parameters.
n (int) – Number of Monte Carlo samples.
rng (np.random.Generator)

Return type:

list of PetrophysicalConfig, length n

class pycsamt.interp.UncertaintyResult(resistivity_model, config, bounds, n_samples, method_tag='', sampled_params=<factory>, mean_K=<factory>, std_K=<factory>, p10_K=<factory>, p50_K=<factory>, p90_K=<factory>, cv_K=<factory>, mean_Sw=<factory>, std_Sw=<factory>, p10_Sw=<factory>, p90_Sw=<factory>, mean_phi=<factory>, std_phi=<factory>, mean_wt=<factory>, std_wt=<factory>, p10_wt=<factory>, p50_wt=<factory>, p90_wt=<factory>, wt_detection_rate=<factory>, mean_T=<factory>, std_T=<factory>, p10_T=<factory>, p50_T=<factory>, p90_T=<factory>)#

Bases: PyCSAMTObject

Ensemble statistics from a MonteCarloHydro run.

All 2-D arrays have shape (n_z, n_x); 1-D station arrays have shape (n_x,). K values are stored on the linear scale (m/s).

Variables:

resistivity_model (ResistivityModel)
config (PetrophysicalConfig) – Central configuration used to seed the Monte Carlo.
bounds (UncertaintyBounds)
n_samples (int)
method_tag (str)
sampled_params (ndarray (n_samples, n_free)) – Matrix of sampled parameter values, columns in bounds.free_names order.
(m/s) (Hydraulic conductivity K)
------------------------------
p90_K (mean_K, std_K, p10_K, p50_K,)
cv_K (ndarray (n_z, n_x)) – Coefficient of variation std_K / mean_K. Dimensionless uncertainty map.
Sw (Water saturation)
-------------------
p90_Sw (mean_Sw, std_Sw, p10_Sw,)
φ (Porosity)
----------
std_phi (mean_phi,)
(m) (Water-table depth)
---------------------
p90_wt (mean_wt, std_wt, p10_wt, p50_wt,)
wt_detection_rate (ndarray (n_x,)) – Fraction of samples in which the water table was detected (0–1).
(m²/s) (Transmissivity T)
-----------------------
p90_T (mean_T, std_T, p10_T, p50_T,)

Parameters:

resistivity_model (ResistivityModel)
config (PetrophysicalConfig)
bounds (UncertaintyBounds)
n_samples (int)
method_tag (str)
sampled_params (ndarray)
mean_K (ndarray)
std_K (ndarray)
p10_K (ndarray)
p50_K (ndarray)
p90_K (ndarray)
cv_K (ndarray)
mean_Sw (ndarray)
std_Sw (ndarray)
p10_Sw (ndarray)
p90_Sw (ndarray)
mean_phi (ndarray)
std_phi (ndarray)
mean_wt (ndarray)
std_wt (ndarray)
p10_wt (ndarray)
p50_wt (ndarray)
p90_wt (ndarray)
wt_detection_rate (ndarray)
mean_T (ndarray)
std_T (ndarray)
p10_T (ndarray)
p50_T (ndarray)
p90_T (ndarray)

resistivity_model: ResistivityModel#

config: PetrophysicalConfig#

bounds: UncertaintyBounds#

n_samples: int#

method_tag: str = ''#

sampled_params: ndarray#

mean_K: ndarray#

std_K: ndarray#

p10_K: ndarray#

p50_K: ndarray#

p90_K: ndarray#

cv_K: ndarray#

mean_Sw: ndarray#

std_Sw: ndarray#

p10_Sw: ndarray#

p90_Sw: ndarray#

mean_phi: ndarray#

std_phi: ndarray#

mean_wt: ndarray#

std_wt: ndarray#

p10_wt: ndarray#

p50_wt: ndarray#

p90_wt: ndarray#

wt_detection_rate: ndarray#

mean_T: ndarray#

std_T: ndarray#

p10_T: ndarray#

p50_T: ndarray#

p90_T: ndarray#

prob_wt_shallower_than(depth_m, *, wt_ensemble=None)#

P(water-table depth < depth_m) per station column.

Parameters:

depth_m (float) – Reference depth (m, positive downward).
wt_ensemble (ndarray (n_samples, n_x), optional) – Full water-table ensemble from MonteCarloHydro.run_ensemble(). If None, approximated from the stored P10/P50/P90 assuming a Gaussian distribution with mean=mean_wt, std=std_wt.

Return type:

ndarray (n_x,) — probability in [0, 1]

station_report()#

Per-station uncertainty summary.

Return type:: list[dict]

to_csv(path)#

Write per-station uncertainty summary to CSV.

Parameters:: path (str | Path)
Return type:: Path

class pycsamt.interp.MonteCarloHydro(resistivity_model, config, bounds, *, n_samples=200, seed=42, method_tag='', verbose=False)#

Bases: PyCSAMTObject

Monte Carlo uncertainty propagation for a quantitative EM hydro model.

Parameters:

resistivity_model (ResistivityModel) – Source inversion model.
config (PetrophysicalConfig) – Central (best-estimate) petrophysical configuration.
bounds (UncertaintyBounds) – Prior distributions for the free parameters.
n_samples (int) – Number of Monte Carlo realisations (default 200). Typical: 100 for quick runs, 500–1000 for publication-quality.
seed (int) – Random seed for reproducibility (default 42).
method_tag (str, optional) – EM method label inherited by the output.
verbose (bool) – Print progress every 50 iterations (default False).

Examples

>>> bounds = UncertaintyBounds(
...     rho_w_range=(5.0, 100.0), m_range=(1.4, 2.4)
... )
>>> mc = MonteCarloHydro(rm, cfg, bounds, n_samples=300)
>>> unc = mc.run()
>>> unc.cv_K.max()  # worst-case K uncertainty (fraction)
>>> unc.p90_wt - unc.p10_wt  # WT depth spread per station (m)

run()#

Execute the Monte Carlo ensemble and return UncertaintyResult.

Return type:: UncertaintyResult

run_ensemble()#

Like run() but also returns the raw WT and T ensembles.

Returns:

unc (UncertaintyResult)
wt_ensemble (ndarray (n_samples, n_x))
T_ensemble (ndarray (n_samples, n_x))

Return type:

tuple[UncertaintyResult, ndarray, ndarray]

class pycsamt.interp.RockDatabase(entries, *, metadata=None)#

Bases: PyCSAMTObject, MetadataMixin

Extensible rock physics database for EM resistivity interpretation.

Parameters:

entries (list of RockEntry) – The database entries. default() returns the built-in set.
metadata (dict)

Example

>>> db = RockDatabase.default()
>>> db.classify(180.0).name
'Granite (weathered)'

metadata: dict#

property entries: tuple[RockEntry, ...]#: Read-only view of the database entries, in insertion order.

classmethod default()#

Return a database pre-loaded with the built-in rock entries.

Entries come from pycsamt.geology.rock_library.BUILTIN_ROCKS; see that module for the literature the ranges are drawn from.

Return type:: RockDatabase

classmethod from_csv(path)#

Load from a CSV file.

Required columns: name, rho_min, rho_max Optional columns: color, description, code, source

Parameters:: path (str | Path)
Return type:: RockDatabase

classmethod from_provider(provider)#

Build a database from any RockPropertyProvider.

This is the generic entry point behind from_url(); call it directly when provider is something other than a plain URL fetch, such as a provider pre-configured with authentication or a provider composing several sources.

Parameters:: provider (RockPropertyProvider)
Return type:: RockDatabase

classmethod from_url(url, *, cache_dir=None, ttl_seconds=86400.0, timeout=10.0, force=False, fallback=True)#

Fetch a database from url, with local caching and fallback.

url must serve a JSON array of objects using the same fields as RockEntry (name, rho_min, rho_max, and optionally color, description, code, source). No public rock-resistivity service is bundled or assumed by default – point this at a project- or organisation-controlled endpoint (an internal API, a shared JSON file on a file server, and so on).

Parameters:

url (str) – Location to fetch the JSON entry list from.
cache_dir (path-like, optional) – Override the local cache location. Defaults to $PYCSAMT_ROCKDB_CACHE or ~/.pycsamt/rock_db.
ttl_seconds (float) – Reuse a cached response younger than this many seconds instead of re-fetching. Default one day.
timeout (float) – Network timeout in seconds for the fetch itself.
force (bool) – Re-fetch even if a fresh cache entry exists.
fallback (bool) – If the fetch fails (network, timeout, malformed response) and no usable cache entry exists, fall back to default() instead of raising. Set to False to surface the failure instead.

Return type:

RockDatabase

2.17.3. Interpretation Modules#

`pycsamt.interp.calibrate`	ModelCalibrator — constrain a 2-D EM model with borehole ground truth.
`pycsamt.interp.constraints`	Field-measurement constraints for petrophysical calibration.
`pycsamt.interp.export`	export — write EM interpretation results to industry formats.
`pycsamt.interp.fusion`	Multi-method EM fusion — merge shallow and deep resistivity models.
`pycsamt.interp.hydro`	Hydrogeophysical interpretation from EM resistivity models.
`pycsamt.interp.hydromodel`	Quantitative hydro-geophysical model from EM resistivity sections.
`pycsamt.interp.petrophysics`	Petrophysical transforms for EM hydrogeophysics.
`pycsamt.interp.plot`	plot — visualization for pycsamt.interp results.
`pycsamt.interp.timelapse`	Time-lapse EM monitoring for hydrogeological change detection.
`pycsamt.interp.uncertainty`	Monte Carlo uncertainty propagation for EM hydro-geophysical models.