pycsamt.interp#

Geological and hydrogeological interpretation helpers, lithology catalogues, borehole calibration, uncertainty, fusion, and report/export routines.

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.

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()

Package layout#

_base: ResistivityModel — unified model container; adapters for Occam2D, ModEM, AI outputs, and raw numpy arrays.
borehole: Borehole, Interval — ground-truth depth-interval data; readers for CSV and LAS 2.0.
calibrate: ModelCalibrator — generalised NM construction (Kouadio et al. 2022), borehole-constrained or database-only.
lithology: RockDatabase, StratigraphicLog, Layer — resistivity-to-lithology classification; built-in 25-rock EM database.
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)#

Bases: object

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.

x_centers: ndarray#

z_centers: ndarray#

rho_2d: ndarray#

station_x: ndarray#

station_names: list[str]#

method: str = 'generic'#

rms: float = nan#

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

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

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

class pycsamt.interp.Borehole(name, x, intervals, *, collar_elevation=0.0)#

Bases: object

Borehole / well log with depth-interval data.

Parameters:

name (str) – Well / borehole identifier.
x (float) – Position along the survey profile, metres.
collar_elevation (float) – Surface elevation at the borehole collar, metres a.s.l. Used only when elevation-corrected exports are requested.
intervals (list of Interval) – Depth-interval log, sorted ascending by top.

interval_at_depth(z)#

Return the interval that contains depth z, or None.

Parameters:: z (float)
Return type:: Interval | None

tres_at_depth(z)#

Return TRES (Ω·m) at depth z, or None if unknown.

Parameters:: z (float)
Return type:: float | None

lithology_at_depth(z)#

Return the lithology name at depth z, or None.

Parameters:: z (float)
Return type:: str | None

tres_column(z_centers)#

Return TRES values at z_centers as a float array.

Depths not covered by any interval are nan.

Parameters:: z_centers (ndarray)
Return type:: ndarray

property max_depth: float#

property min_depth: float#

classmethod from_csv(path, *, name=None, x=0.0, collar_elevation=0.0, delimiter=',', top_col='top', bottom_col='bottom', lithology_col='lithology', resistivity_col='resistivity')#

Load from a CSV file.

Expected columns (case-insensitive header): top, bottom, lithology[, resistivity]

Parameters:

path (path-like)
name (str, optional) – Defaults to the file stem.
x (float) – Profile position.
collar_elevation (float)
delimiter (str)
top_col (str) – Column header names.
bottom_col (str) – Column header names.
lithology_col (str) – Column header names.
resistivity_col (str) – Column header names.

Return type:

Borehole

classmethod from_las(path, *, x=0.0, collar_elevation=0.0, depth_curve='DEPT', resistivity_curve='RESD', lithology_curve='LITH', null_value=-9999.25, step=None)#

Load from a LAS 2.0 well-log file.

Converts the continuous depth log into discrete intervals by grouping consecutive samples with the same lithology code.

Parameters:

path (path-like)
x (float) – Profile position.
depth_curve (str) – Curve mnemonics. lithology_curve may be None to assign a generic label.
resistivity_curve (str) – Curve mnemonics. lithology_curve may be None to assign a generic label.
lithology_curve (str | None) – Curve mnemonics. lithology_curve may be None to assign a generic label.
null_value (float) – LAS null sentinel replaced with nan.
step (float, optional) – If provided, override the step value from the LAS header.
collar_elevation (float)

Return type:

Borehole

to_dataframe()#: Return intervals as a pandas.DataFrame.

to_dict()#

Return type:: dict

class pycsamt.interp.Interval(top, bottom, lithology, resistivity=None)#

Bases: object

A single depth interval in a borehole log.

Parameters:

top (float) – Depth to the top of the interval, metres.
bottom (float) – Depth to the bottom of the interval, metres.
lithology (str) – Geological formation / lithology name.
resistivity (float or None) – True resistivity (TRES) in Ω·m (linear, not log₁₀). None when no electrical measurement is available.

top: float#

bottom: float#

lithology: str#

resistivity: float | None = None#

property thickness: float#

contains(z)#

Parameters:: z (float)
Return type:: bool

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

Bases: object

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

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'

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)#

Bases: object

Extensible rock physics database for EM resistivity interpretation.

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

Example

>>> db = RockDatabase.default()
>>> db.classify(180.0).name
'Fractured zone'

classmethod default()#

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

Return type:: RockDatabase

classmethod from_csv(path)#

Load from a CSV file.

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

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

classify(rho_ohm_m, method='nearest')#

Return the best-matching rock entry for rho_ohm_m.

Parameters:

rho_ohm_m (float) – Resistivity in Ω·m (linear).
method ({'nearest', 'overlap'}) – 'nearest': log-distance to midpoint. 'overlap': first entry whose range brackets rho_ohm_m.

Return type:

RockEntry

classify_column(rho_log10)#

Classify every cell in a log10-rho depth column.

Parameters:: rho_log10 (ndarray)
Return type:: list[RockEntry]

class pycsamt.interp.RockEntry(name, rho_min, rho_max, color='#AAAAAA', description='', code=0)#

Bases: object

A single entry in the rock physics database.

Parameters:

name (str) – Geological unit / lithology name.
rho_min (float) – Resistivity range in Ω·m (linear scale).
rho_max (float) – Resistivity range in Ω·m (linear scale).
color (str) – Hex colour code for plotting (e.g. '#28B463').
description (str) – Optional free-text note.
code (int) – Integer code used in LAS exports.

name: str#

rho_min: float#

rho_max: float#

color: str = '#AAAAAA'#

description: str = ''#

code: int = 0#

property rho_mid: float#

property log_rho_mid: float#

contains(rho_ohm_m)#

Parameters:: rho_ohm_m (float)
Return type:: bool

class pycsamt.interp.Layer(top, bottom, rho_log10, lithology, color='#AAAAAA', confidence=1.0)#

Bases: object

One geological unit in a pseudo-stratigraphic log.

Parameters:

top (float) – Depth in metres (positive downward).
bottom (float) – Depth in metres (positive downward).
rho_log10 (float) – Representative \(\log_{10}(\rho)\) of the layer.
lithology (str) – Rock name from RockDatabase.
color (str) – Hex colour for plotting.
confidence (float) – Fraction of depth cells whose DB classification matches the reported lithology (0 – 1).

top: float#

bottom: float#

rho_log10: float#

lithology: str#

color: str = '#AAAAAA'#

confidence: float = 1.0#

property thickness: float#

property rho_ohm_m: float#

class pycsamt.interp.StratigraphicLog(station_name, station_x, z_centers, rho_log10, layers)#

Bases: object

Per-station pseudo-stratigraphic depth profile.

Constructed from a 1-D log10-rho column and a RockDatabase, it merges adjacent cells that share the same lithology into discrete Layer objects.

Parameters:

station_name (str)
station_x (float)
z_centers (ndarray (n_z,)) – Depth cell centres, metres.
rho_log10 (ndarray (n_z,)) – \(\log_{10}(\rho)\) for each depth cell.
layers (list of Layer) – Merged geological units (assembled by from_column()).

classmethod from_column(station_name, x, z_centers, rho_log10, db=None, *, merge_tolerance=0.2)#

Build a log from a 1-D resistivity column.

Parameters:

station_name (str)
x (float) – Station position, metres.
z_centers (ndarray (n_z,))
rho_log10 (ndarray (n_z,))
db (RockDatabase, optional) – Defaults to RockDatabase.default().
merge_tolerance (float) – Log10-rho difference threshold for merging adjacent cells into one layer (default 0.2 decade).

Return type:

StratigraphicLog

to_dataframe()#: Return layers as a pandas.DataFrame.

to_dict()#

Return type:: dict

Interpretation Modules#

`pycsamt.interp.borehole`	Borehole — ground-truth data model for EM geological interpretation.
`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.lithology`	Lithology — resistivity-to-geology classification for EM methods.
`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.