PyHydroGeophysX.inversion package#

Submodules#

PyHydroGeophysX.inversion.base module#

Base classes for geophysical inversion frameworks.

This module defines:

  • InversionResult: A base class to store and manage common results from an inversion process, including final model, predicted data, convergence history, and plotting utilities.

  • TimeLapseInversionResult: A specialized version of InversionResult for time-lapse studies, handling multiple models over time and providing time-slice plotting and animation.

  • InversionBase: An abstract base class outlining the common structure and interface for various geophysical inversion methods (e.g., ERT, SRT). It handles data, mesh, and basic parameter management.

class PyHydroGeophysX.inversion.base.InversionBase(data: pygimli.DataContainer, mesh: pygimli.Mesh | None = None, **kwargs: Any)[source]#

Bases: object

Abstract base class for geophysical inversion methods.

This class provides a foundational structure for specific inversion techniques (e.g., ERT, SRT). It manages observed data, mesh, common inversion parameters, and an InversionResult object. Subclasses must implement methods for setting up the inversion, running the inversion loop, computing Jacobians, and evaluating the objective function.

compute_jacobian(model: ndarray) ndarray[source]#

Abstract method to compute the Jacobian matrix (sensitivity matrix).

The Jacobian J_ij = ∂d_i / ∂m_j relates changes in model parameters (m_j) to changes in observed data (d_i).

Parameters:

model (np.ndarray) – The current model parameter vector for which to compute the Jacobian.

Returns:,

np.ndarray: The computed Jacobian matrix, typically of shape (n_data, n_model_params).

objective_function(model: ndarray, data_to_fit: ndarray | None = None) float[source]#

Abstract method to compute the value of the objective function.

The objective function typically includes a data misfit term and one or more regularization terms: Φ(m) = Φ_d(m) + λ * Φ_m(m).

Parameters:
  • model (np.ndarray) – The current model parameter vector.

  • data_to_fit (Optional[np.ndarray], optional) – The observed data to fit. If None, self.data (from DataContainer) is typically used. Defaults to None.

Returns:

The calculated value of the objective function.

Return type:

float

run() InversionResult[source]#

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() None[source]#

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

class PyHydroGeophysX.inversion.base.InversionResult[source]#

Bases: object

Base class to store, save, load, and plot results from a geophysical inversion.

final_model#

The final inverted model parameters (e.g., resistivity, velocity). Typically a 1D array corresponding to mesh cells.

Type:

Optional[np.ndarray]

predicted_data#

The data predicted by the forward model using the final_model.

Type:

Optional[np.ndarray]

coverage#

Coverage or sensitivity values for the model parameters, often derived from the Jacobian or resolution matrix.

Type:

Optional[np.ndarray]

mesh#

The PyGIMLi mesh object used in the inversion.

Type:

Optional[pg.Mesh]

iteration_models#

A list storing the model parameters at each iteration of the inversion.

Type:

List[np.ndarray]

iteration_data_errors#

A list storing the data misfit (e.g., residuals) at each iteration.

Type:

List[np.ndarray]

iteration_chi2#

A list storing the chi-squared (χ²) value or a similar misfit metric at each iteration.

Type:

List[float]

meta#

A dictionary to store any additional metadata about the inversion run (e.g., inversion parameters, timings, comments).

Type:

Dict[str, Any]

classmethod load(filename: str) InversionResult[source]#

Load inversion results from a file previously saved by the save method.

Parameters:

filename (str) – The base path to the saved results file. If saved as name.pkl and name.bms, provide name.pkl or name.

Returns:

An instance of InversionResult (or a subclass if called on one)

populated with the loaded data.

Return type:

InversionResult

Raises:
  • FileNotFoundError – If the main data file or associated mesh file (if referenced) is not found.

  • IOError – If there’s an error during file reading.

  • pickle.UnpicklingError – If the file cannot be unpickled.

plot_convergence(ax: Axes | None = None, **kwargs: Any) Tuple[Figure, Axes][source]#

Plot the convergence curve (chi-squared misfit vs. iteration number).

Parameters:
  • ax (Optional[plt.Axes], optional) – A matplotlib Axes object to plot on. If None, a new figure and axes are created. Defaults to None.

  • **kwargs (Any) – Additional keyword arguments passed directly to ax.plot (e.g., color, marker, linestyle).

Returns:

The matplotlib Figure and Axes objects of the plot.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If self.iteration_chi2 is empty (no convergence data).

plot_model(ax: Axes | None = None, cmap: str = 'viridis', coverage_threshold: float | None = None, **kwargs: Any) Tuple[Figure, Axes][source]#

Plot the final inverted model on its associated mesh.

Parameters:
  • ax (Optional[plt.Axes], optional) – A matplotlib Axes object to plot on. If None, a new figure and axes are created. Defaults to None.

  • cmap (str, optional) – The colormap to use for visualizing model values. Defaults to ‘viridis’.

  • coverage_threshold (Optional[float], optional) – If provided, cells with coverage values below this threshold will be masked (made transparent or semi-transparent) in the plot. Requires self.coverage to be populated. Defaults to None (no masking).

  • **kwargs (Any) – Additional keyword arguments passed directly to pygimli.show (e.g., cMin, cMax, orientation, logScale).

Returns:

The matplotlib Figure and Axes objects of the plot.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If self.final_model or self.mesh is None.

save(filename: str) None[source]#

Save the inversion results to a file using Python’s pickle format. If a mesh is present, it is saved separately as a PyGIMLi binary mesh file (.bms).

Parameters:

filename (str) – The base path (including filename without extension) to save the results. The main data will be saved as filename.pkl (or just filename if user includes .pkl). The mesh will be saved as filename.bms or filename.pkl.bms. It’s recommended to provide filename without .pkl.

Raises:
  • IOError – If there’s an error during file writing.

  • pickle.PicklingError – If an object cannot be pickled.

class PyHydroGeophysX.inversion.base.TimeLapseInversionResult[source]#

Bases: InversionResult

Specialized class to store and manage results from time-lapse inversions. Inherits from InversionResult and adds attributes specific to time-lapse data.

final_models#

A 2D NumPy array where each column represents the inverted model for a specific timestep (shape: num_cells x num_timesteps).

Type:

Optional[np.ndarray]

timesteps#

An array or list of time values (e.g., hours, days) corresponding to each model slice in final_models.

Type:

Optional[np.ndarray]

all_coverage#

A list where each element is the coverage array for the corresponding timestep’s model.

Type:

List[np.ndarray]

all_chi2#

A list where each element is a list of chi-squared values per iteration for the inversion of that specific timestep or window. (Note: Original was List[float], if chi2 is per window, might need adjustment)

Type:

List[Any]

create_time_lapse_animation(output_filename: str, cmap: str = 'viridis', coverage_threshold: float | None = None, dpi: int = 100, fps: int = 2, **kwargs: Any) None[source]#

Create and save an animation (e.g., MP4 video) of the time-lapse inversion results.

Requires ffmpeg or another Matplotlib-supported animation writer to be installed.

Parameters:
  • output_filename (str) – The filename for the output animation (e.g., ‘timelapse_animation.mp4’).

  • cmap (str, optional) – Colormap for model values. Defaults to ‘viridis’.

  • coverage_threshold (Optional[float], optional) – Threshold for coverage masking. Defaults to None.

  • dpi (int, optional) – Dots Per Inch for the output animation. Defaults to 100.

  • fps (int, optional) – Frames Per Second for the animation. Defaults to 2.

  • **kwargs (Any) – Additional keyword arguments passed to plot_time_slice for each frame.

Raises:
  • ValueError – If final_models or mesh is missing.

  • ImportError – If matplotlib.animation cannot be imported.

classmethod load(filename: str) TimeLapseInversionResult[source]#

Load time-lapse inversion results. Overrides base class load.

plot_time_slice(timestep_idx: int, ax: Axes | None = None, cmap: str = 'viridis', coverage_threshold: float | None = None, **kwargs: Any) Tuple[Figure, Axes][source]#

Plot a single time slice (inverted model at a specific timestep) from the results.

Parameters:
  • timestep_idx (int) – The zero-based index of the timestep to plot.

  • ax (Optional[plt.Axes], optional) – Matplotlib Axes to plot on. If None, creates new.

  • cmap (str, optional) – Colormap for model values. Defaults to ‘viridis’.

  • coverage_threshold (Optional[float], optional) – Threshold for coverage masking. Defaults to None.

  • **kwargs (Any) – Additional arguments passed to pg.show.

Returns:

The matplotlib Figure and Axes objects.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If models or mesh are missing, or if timestep_idx is out of range.

save(filename: str) None[source]#

Save time-lapse inversion results. Overrides base class save.

PyHydroGeophysX.inversion.cross_constraints module#

Cross-method constraint utilities for joint/cooperative inversion.

These utilities enable information sharing between geophysical methods and provide hydrology-to-geophysics coupling helpers.

class PyHydroGeophysX.inversion.cross_constraints.PetrophysicalCoupling[source]#

Bases: object

Coupling helpers from hydrological state to multi-method geophysics.

static compare_inversions_to_hydro(hydro_wc, ert_result, srt_result=None, em_result=None, petro_params=None)[source]#

Compare inversion products back to hydrological water content.

Returns misfit statistics for available methods.

static water_content_to_all_geophysics(water_content, porosity, **params)[source]#

Convert water content to resistivity, velocity, and conductivity.

Uses existing petrophysical functions in the package.

class PyHydroGeophysX.inversion.cross_constraints.StructuralConstraint[source]#

Bases: object

Build structural constraint terms from one method for another.

static apply_structural_weights_to_Wm(Wm, boundary_weights)[source]#

Modify an existing smoothness matrix Wm to respect boundaries.

static build_cross_gradient_operator(mesh, model_a, model_b)[source]#

Build a cross-gradient operator surrogate.

Returns a diagonal sparse matrix weighted by local cross-gradient magnitude, where small values indicate better structural agreement.

static build_linearized_cross_gradient_blocks(RCM: ndarray, X: ndarray, model_a: ndarray, model_b: ndarray, mode: str = 'direct') Tuple[ndarray, ndarray][source]#

Build linearized cross-gradient blocks B1 and B2.

B1 multiplies model_a and B2 multiplies model_b. The resulting penalty terms are ||B1 m_a||^2 and ||B2 m_b||^2.

static build_local_design_matrix(mesh) ndarray[source]#

Build the local design matrix X=[x, z, 1] used by cross-gradient.

This mirrors the legacy xyzpos[:, 2] = 1 construction.

static build_neighborhood_matrix(mesh, Wm: Any | None = None, source: str = 'smoothness', correlation_lengths: Sequence[float] = (4.0, 4.0), threshold: float = 0.01, binarize: bool = True) ndarray[source]#

Build the neighborhood/correlation matrix used by cross-gradient.

Parameters:
  • mesh – Mesh object.

  • Wm – Smoothness matrix. Required when source='smoothness'.

  • source – One of 'smoothness' (legacy direct mode) or 'covariance' (legacy spatial mode).

  • correlation_lengths – Correlation lengths passed to pg.utils.covarianceMatrix when source='covariance'.

  • threshold – Entries with absolute value below this threshold are zeroed.

  • binarize – If True, convert nonzero entries to 1.

static from_conductivity_model(conductivity, mesh, gradient_threshold: float = 0.3)[source]#

Extract structural boundaries from an EM conductivity model.

static from_velocity_model(velocity_model, mesh, gradient_threshold: float = 0.3)[source]#

Extract structural boundaries from an SRT velocity model.

Returns a pairwise-boundary weight map that can reduce smoothing across detected structural interfaces.

PyHydroGeophysX.inversion.ert_inversion module#

Single-time ERT inversion functionality.

class PyHydroGeophysX.inversion.ert_inversion.ERTInversion(data_file: str, mesh: pygimli.Mesh | None = None, **kwargs)[source]#

Bases: InversionBase

Single-time ERT inversion class.

run(initial_model: ndarray | None = None) InversionResult[source]#

Run ERT inversion.

Parameters:

initial_model – Initial model parameters (if None, a homogeneous model is used)

Returns:

InversionResult with inversion results

setup()[source]#

Set up ERT inversion (create operators, matrices, etc.)

PyHydroGeophysX.inversion.fdem_inversion module#

Inversion utilities for Frequency-Domain Electromagnetic (FDEM) data.

Provides 1D layered-earth FDEM inversion using SimPEG.

class PyHydroGeophysX.inversion.fdem_inversion.FDEMInversion(frequencies: ndarray, dobs: ndarray, uncertainties: ndarray, thicknesses: ndarray | None = None, source_location: ndarray | None = None, source_radius: float = 10.0, receiver_location: ndarray | None = None, receiver_orientation: str = 'z', receiver_component: str = 'secondary', waveform_type: str = 'dipole', n_layers: int = 25, min_thickness: float = 1.0, max_thickness: float = 30.0, **kwargs)[source]#

Bases: object

1D FDEM inversion using SimPEG.

Follows the same pattern used in TDEMInversion.

run(starting_model: ndarray | None = None) FDEMInversionResult[source]#
setup() None[source]#
class PyHydroGeophysX.inversion.fdem_inversion.FDEMInversionResult(recovered_model: ndarray | None = None, recovered_conductivity: ndarray | None = None, l2_model: ndarray | None = None, l2_conductivity: ndarray | None = None, predicted_data: ndarray | None = None, mesh: Any | None = None, thicknesses: ndarray | None = None, chi2: float | None = None, frequencies: ndarray | None = None)[source]#

Bases: object

Container for FDEM inversion results.

chi2: float = None#
frequencies: ndarray = None#
l2_conductivity: ndarray = None#
l2_model: ndarray = None#
mesh: Any = None#
predicted_data: ndarray = None#
recovered_conductivity: ndarray = None#
recovered_model: ndarray = None#
thicknesses: ndarray = None#

PyHydroGeophysX.inversion.joint_ert_srt module#

Joint ERT-SRT inversion with structural and geostatistical constraints.

This module integrates the legacy alternating inversion strategy from cg_joint_inversion_hang into the package architecture.

class PyHydroGeophysX.inversion.joint_ert_srt.JointERTSRTInversion(ert_data: str | PathLike | pygimli.DataContainer, srt_data: str | PathLike | pygimli.DataContainer, mesh: pygimli.Mesh | None = None, **kwargs: Any)[source]#

Bases: InversionBase

Joint ERT-SRT inversion using alternating Gauss-Newton updates.

  • ERT model parameter: log-resistivity

  • SRT model parameter: log-slowness

  • Coupling: linearized cross-gradient terms (B1, B2)

  • Regularization: smoothness or geostatistical

run() JointERTSRTResult[source]#

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() None[source]#

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

class PyHydroGeophysX.inversion.joint_ert_srt.JointERTSRTResult(ert_log_resistivity: ~numpy.ndarray | None = None, srt_log_slowness: ~numpy.ndarray | None = None, ert_resistivity: ~numpy.ndarray | None = None, srt_velocity: ~numpy.ndarray | None = None, ert_predicted: ~numpy.ndarray | None = None, srt_predicted: ~numpy.ndarray | None = None, ert_coverage: ~numpy.ndarray | None = None, srt_coverage: ~numpy.ndarray | None = None, chi2_ert: float | None = None, chi2_srt: float | None = None, iteration_history: list = <factory>, mesh: ~typing.Any | None = None, meta: ~typing.Dict[str, ~typing.Any] = <factory>)[source]#

Bases: object

Container for joint ERT-SRT inversion outputs.

chi2_ert: float | None = None#
chi2_srt: float | None = None#
ert_coverage: ndarray | None = None#
ert_log_resistivity: ndarray | None = None#
ert_predicted: ndarray | None = None#
ert_resistivity: ndarray | None = None#
iteration_history: list#
mesh: Any = None#
meta: Dict[str, Any]#
srt_coverage: ndarray | None = None#
srt_log_slowness: ndarray | None = None#
srt_predicted: ndarray | None = None#
srt_velocity: ndarray | None = None#

PyHydroGeophysX.inversion.multi_method module#

Unified multi-method geophysical inversion interface.

class PyHydroGeophysX.inversion.multi_method.GeophysicalInversion(method: str, **kwargs)[source]#

Bases: object

Unified factory for multi-method geophysical inversion.

Dispatches to the correct inversion engine.

SUPPORTED = {'ert', 'fdem', 'joint', 'joint_ert_srt', 'srt', 'tdem'}#
property engine#
run(**kwargs)[source]#
setup(**kwargs)[source]#

PyHydroGeophysX.inversion.srt_inversion module#

Seismic Refraction Tomography (SRT) inversion functionality.

Uses PyGIMLi’s TravelTimeManager for forward modeling and Jacobian computation. Provides custom Gauss-Newton inversion with the same architecture as ERTInversion but for travel-time data.

class PyHydroGeophysX.inversion.srt_inversion.SRTInversion(data_file: str, mesh: pygimli.Mesh | None = None, **kwargs: Any)[source]#

Bases: InversionBase

Seismic Refraction Tomography inversion class.

Inverts travel-time data for subsurface velocity via log-slowness, and uses a Gauss-Newton optimization loop.

run(initial_model: ndarray | None = None) InversionResult[source]#

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() None[source]#

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

PyHydroGeophysX.inversion.srt_time_lapse module#

Time-lapse Seismic Refraction Tomography (SRT) inversion functionality.

Jointly inverts multiple travel-time datasets with spatial and temporal regularization using log-slowness parameterization.

class PyHydroGeophysX.inversion.srt_time_lapse.TimeLapseSRTInversion(data_files: List[str], measurement_times: List[float], mesh: pygimli.Mesh | None = None, **kwargs: Any)[source]#

Bases: InversionBase

Time-lapse Seismic Refraction Tomography inversion.

Architecture mirrors TimeLapseERTInversion while using TravelTimeManager and log-slowness model parameters.

run(initial_model: ndarray | None = None) TimeLapseInversionResult[source]#

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() None[source]#

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

PyHydroGeophysX.inversion.tdem_inversion module#

Inversion utilities for Time-Domain Electromagnetic (TDEM) data.

This module provides classes for 1D TDEM inversion using SimPEG, with support for both L2 and sparse (IRLS) regularization.

class PyHydroGeophysX.inversion.tdem_inversion.TDEMInversion(times: ndarray, dobs: ndarray, uncertainties: ndarray, source_radius: float = 10.0, source_location: ndarray | None = None, receiver_location: ndarray | None = None, n_layers: int = 25, min_thickness: float = 1.0, max_thickness: float = 30.0, **kwargs)[source]#

Bases: object

Class for 1D TDEM inversion using SimPEG.

This class provides functionality for inverting TDEM sounding data to recover 1D layered Earth conductivity models.

Example

>>> # Load or create data
>>> times = np.logspace(-5, -2, 31)
>>> dobs = ...  # observed data
>>> uncertainties = ...  # data uncertainties
>>>
>>> # Create inversion
>>> inv = TDEMInversion(
...     times=times,
...     dobs=dobs,
...     uncertainties=uncertainties,
...     source_radius=10.0
... )
>>>
>>> # Run inversion
>>> result = inv.run()
plot_result(result: TDEMInversionResult, true_model: Tuple[ndarray, ndarray] | None = None, save_path: str | None = None) None[source]#

Plot inversion results.

Parameters:
  • result – TDEMInversionResult from run()

  • true_model – Tuple of (thicknesses, conductivity) for true model

  • save_path – Path to save figure (optional)

run(starting_model: ndarray | None = None) TDEMInversionResult[source]#

Run TDEM inversion.

Parameters:

starting_model – Initial log-conductivity model (optional)

Returns:

TDEMInversionResult with inversion results

setup() None[source]#

Set up inversion components (survey, mesh, simulation).

class PyHydroGeophysX.inversion.tdem_inversion.TDEMInversionResult(recovered_model: ~numpy.ndarray | None = None, recovered_conductivity: ~numpy.ndarray | None = None, l2_model: ~numpy.ndarray | None = None, l2_conductivity: ~numpy.ndarray | None = None, predicted_data: ~numpy.ndarray | None = None, mesh: discretize.TensorMesh | None = None, thicknesses: ~numpy.ndarray | None = None, chi2: float | None = None, iterations: int = 0, convergence_history: ~typing.List[float] = <factory>)[source]#

Bases: object

Container for TDEM inversion results.

recovered_model#

Final recovered model (log-conductivity)

Type:

numpy.ndarray

recovered_conductivity#

Recovered conductivity (S/m)

Type:

numpy.ndarray

l2_model#

L2 model before IRLS (log-conductivity)

Type:

numpy.ndarray

l2_conductivity#

L2 conductivity (S/m)

Type:

numpy.ndarray

predicted_data#

Predicted data from recovered model

Type:

numpy.ndarray

mesh#

TensorMesh used for inversion

Type:

discretize.TensorMesh

thicknesses#

Layer thicknesses used in inversion

Type:

numpy.ndarray

chi2#

Final chi-squared misfit

Type:

float

iterations#

Number of iterations

Type:

int

convergence_history#

History of data misfit per iteration

Type:

List[float]

chi2: float = None#
convergence_history: List[float]#
iterations: int = 0#
l2_conductivity: ndarray = None#
l2_model: ndarray = None#
mesh: discretize.TensorMesh = None#
predicted_data: ndarray = None#
recovered_conductivity: ndarray = None#
recovered_model: ndarray = None#
thicknesses: ndarray = None#
PyHydroGeophysX.inversion.tdem_inversion.run_tdem_inversion(times: ndarray, dobs: ndarray, uncertainties: ndarray, source_radius: float = 10.0, n_layers: int = 25, use_irls: bool = True, verbose: bool = True, **kwargs) TDEMInversionResult[source]#

Convenience function to run TDEM inversion.

Parameters:
  • times – Time channels (s)

  • dobs – Observed data (T)

  • uncertainties – Data uncertainties (T)

  • source_radius – Loop radius (m)

  • n_layers – Number of inversion layers

  • use_irls – Use IRLS for sparse inversion

  • verbose – Print progress

  • **kwargs – Additional parameters for TDEMInversion

Returns:

TDEMInversionResult with inversion results

PyHydroGeophysX.inversion.time_lapse module#

Time-lapse ERT inversion functionality.

class PyHydroGeophysX.inversion.time_lapse.TimeLapseERTInversion(data_files: List[str], measurement_times: List[float], mesh: pygimli.Mesh | None = None, **kwargs)[source]#

Bases: InversionBase

Time-lapse ERT inversion class.

run(initial_model: ndarray | None = None) TimeLapseInversionResult[source]#

Run time-lapse ERT inversion.

Parameters:

initial_model – Initial model parameters (if None, a homogeneous model is used)

Returns:

TimeLapseInversionResult with inversion results

setup()[source]#

Set up time-lapse ERT inversion (load data, create operators, matrices, etc.)

PyHydroGeophysX.inversion.windowed module#

Windowed time-lapse ERT inversion for handling large temporal datasets.

class PyHydroGeophysX.inversion.windowed.WindowedTimeLapseERTInversion(data_dir: str, ert_files: List[str], measurement_times: List[float], window_size: int = 3, mesh: pygimli.Mesh | str | None = None, **kwargs)[source]#

Bases: object

Class for windowed time-lapse ERT inversion to handle large temporal datasets.

run(window_parallel: bool = False, max_window_workers: int | None = None) TimeLapseInversionResult[source]#

Run windowed time-lapse ERT inversion.

Parameters:
  • window_parallel – Whether to process windows in parallel

  • max_window_workers – Maximum number of parallel workers (None for auto)

Returns:

TimeLapseInversionResult with stitched results

Module contents#

Inversion framework for geophysical applications.

class PyHydroGeophysX.inversion.ERTInversion(data_file: str, mesh: pygimli.Mesh | None = None, **kwargs)[source]#

Bases: InversionBase

Single-time ERT inversion class.

run(initial_model: ndarray | None = None) InversionResult[source]#

Run ERT inversion.

Parameters:

initial_model – Initial model parameters (if None, a homogeneous model is used)

Returns:

InversionResult with inversion results

setup()[source]#

Set up ERT inversion (create operators, matrices, etc.)

class PyHydroGeophysX.inversion.FDEMInversion(frequencies: ndarray, dobs: ndarray, uncertainties: ndarray, thicknesses: ndarray | None = None, source_location: ndarray | None = None, source_radius: float = 10.0, receiver_location: ndarray | None = None, receiver_orientation: str = 'z', receiver_component: str = 'secondary', waveform_type: str = 'dipole', n_layers: int = 25, min_thickness: float = 1.0, max_thickness: float = 30.0, **kwargs)[source]#

Bases: object

1D FDEM inversion using SimPEG.

Follows the same pattern used in TDEMInversion.

run(starting_model: ndarray | None = None) FDEMInversionResult[source]#
setup() None[source]#
class PyHydroGeophysX.inversion.FDEMInversionResult(recovered_model: ndarray | None = None, recovered_conductivity: ndarray | None = None, l2_model: ndarray | None = None, l2_conductivity: ndarray | None = None, predicted_data: ndarray | None = None, mesh: Any | None = None, thicknesses: ndarray | None = None, chi2: float | None = None, frequencies: ndarray | None = None)[source]#

Bases: object

Container for FDEM inversion results.

chi2: float = None#
frequencies: ndarray = None#
l2_conductivity: ndarray = None#
l2_model: ndarray = None#
mesh: Any = None#
predicted_data: ndarray = None#
recovered_conductivity: ndarray = None#
recovered_model: ndarray = None#
thicknesses: ndarray = None#
class PyHydroGeophysX.inversion.GeophysicalInversion(method: str, **kwargs)[source]#

Bases: object

Unified factory for multi-method geophysical inversion.

Dispatches to the correct inversion engine.

SUPPORTED = {'ert', 'fdem', 'joint', 'joint_ert_srt', 'srt', 'tdem'}#
property engine#
run(**kwargs)[source]#
setup(**kwargs)[source]#
class PyHydroGeophysX.inversion.InversionBase(data: pygimli.DataContainer, mesh: pygimli.Mesh | None = None, **kwargs: Any)[source]#

Bases: object

Abstract base class for geophysical inversion methods.

This class provides a foundational structure for specific inversion techniques (e.g., ERT, SRT). It manages observed data, mesh, common inversion parameters, and an InversionResult object. Subclasses must implement methods for setting up the inversion, running the inversion loop, computing Jacobians, and evaluating the objective function.

compute_jacobian(model: ndarray) ndarray[source]#

Abstract method to compute the Jacobian matrix (sensitivity matrix).

The Jacobian J_ij = ∂d_i / ∂m_j relates changes in model parameters (m_j) to changes in observed data (d_i).

Parameters:

model (np.ndarray) – The current model parameter vector for which to compute the Jacobian.

Returns:,

np.ndarray: The computed Jacobian matrix, typically of shape (n_data, n_model_params).

objective_function(model: ndarray, data_to_fit: ndarray | None = None) float[source]#

Abstract method to compute the value of the objective function.

The objective function typically includes a data misfit term and one or more regularization terms: Φ(m) = Φ_d(m) + λ * Φ_m(m).

Parameters:
  • model (np.ndarray) – The current model parameter vector.

  • data_to_fit (Optional[np.ndarray], optional) – The observed data to fit. If None, self.data (from DataContainer) is typically used. Defaults to None.

Returns:

The calculated value of the objective function.

Return type:

float

run() InversionResult[source]#

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() None[source]#

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

class PyHydroGeophysX.inversion.InversionResult[source]#

Bases: object

Base class to store, save, load, and plot results from a geophysical inversion.

final_model#

The final inverted model parameters (e.g., resistivity, velocity). Typically a 1D array corresponding to mesh cells.

Type:

Optional[np.ndarray]

predicted_data#

The data predicted by the forward model using the final_model.

Type:

Optional[np.ndarray]

coverage#

Coverage or sensitivity values for the model parameters, often derived from the Jacobian or resolution matrix.

Type:

Optional[np.ndarray]

mesh#

The PyGIMLi mesh object used in the inversion.

Type:

Optional[pg.Mesh]

iteration_models#

A list storing the model parameters at each iteration of the inversion.

Type:

List[np.ndarray]

iteration_data_errors#

A list storing the data misfit (e.g., residuals) at each iteration.

Type:

List[np.ndarray]

iteration_chi2#

A list storing the chi-squared (χ²) value or a similar misfit metric at each iteration.

Type:

List[float]

meta#

A dictionary to store any additional metadata about the inversion run (e.g., inversion parameters, timings, comments).

Type:

Dict[str, Any]

classmethod load(filename: str) InversionResult[source]#

Load inversion results from a file previously saved by the save method.

Parameters:

filename (str) – The base path to the saved results file. If saved as name.pkl and name.bms, provide name.pkl or name.

Returns:

An instance of InversionResult (or a subclass if called on one)

populated with the loaded data.

Return type:

InversionResult

Raises:
  • FileNotFoundError – If the main data file or associated mesh file (if referenced) is not found.

  • IOError – If there’s an error during file reading.

  • pickle.UnpicklingError – If the file cannot be unpickled.

plot_convergence(ax: Axes | None = None, **kwargs: Any) Tuple[Figure, Axes][source]#

Plot the convergence curve (chi-squared misfit vs. iteration number).

Parameters:
  • ax (Optional[plt.Axes], optional) – A matplotlib Axes object to plot on. If None, a new figure and axes are created. Defaults to None.

  • **kwargs (Any) – Additional keyword arguments passed directly to ax.plot (e.g., color, marker, linestyle).

Returns:

The matplotlib Figure and Axes objects of the plot.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If self.iteration_chi2 is empty (no convergence data).

plot_model(ax: Axes | None = None, cmap: str = 'viridis', coverage_threshold: float | None = None, **kwargs: Any) Tuple[Figure, Axes][source]#

Plot the final inverted model on its associated mesh.

Parameters:
  • ax (Optional[plt.Axes], optional) – A matplotlib Axes object to plot on. If None, a new figure and axes are created. Defaults to None.

  • cmap (str, optional) – The colormap to use for visualizing model values. Defaults to ‘viridis’.

  • coverage_threshold (Optional[float], optional) – If provided, cells with coverage values below this threshold will be masked (made transparent or semi-transparent) in the plot. Requires self.coverage to be populated. Defaults to None (no masking).

  • **kwargs (Any) – Additional keyword arguments passed directly to pygimli.show (e.g., cMin, cMax, orientation, logScale).

Returns:

The matplotlib Figure and Axes objects of the plot.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If self.final_model or self.mesh is None.

save(filename: str) None[source]#

Save the inversion results to a file using Python’s pickle format. If a mesh is present, it is saved separately as a PyGIMLi binary mesh file (.bms).

Parameters:

filename (str) – The base path (including filename without extension) to save the results. The main data will be saved as filename.pkl (or just filename if user includes .pkl). The mesh will be saved as filename.bms or filename.pkl.bms. It’s recommended to provide filename without .pkl.

Raises:
  • IOError – If there’s an error during file writing.

  • pickle.PicklingError – If an object cannot be pickled.

class PyHydroGeophysX.inversion.JointERTSRTInversion(ert_data: str | PathLike | pygimli.DataContainer, srt_data: str | PathLike | pygimli.DataContainer, mesh: pygimli.Mesh | None = None, **kwargs: Any)[source]#

Bases: InversionBase

Joint ERT-SRT inversion using alternating Gauss-Newton updates.

  • ERT model parameter: log-resistivity

  • SRT model parameter: log-slowness

  • Coupling: linearized cross-gradient terms (B1, B2)

  • Regularization: smoothness or geostatistical

run() JointERTSRTResult[source]#

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() None[source]#

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

class PyHydroGeophysX.inversion.JointERTSRTResult(ert_log_resistivity: ~numpy.ndarray | None = None, srt_log_slowness: ~numpy.ndarray | None = None, ert_resistivity: ~numpy.ndarray | None = None, srt_velocity: ~numpy.ndarray | None = None, ert_predicted: ~numpy.ndarray | None = None, srt_predicted: ~numpy.ndarray | None = None, ert_coverage: ~numpy.ndarray | None = None, srt_coverage: ~numpy.ndarray | None = None, chi2_ert: float | None = None, chi2_srt: float | None = None, iteration_history: list = <factory>, mesh: ~typing.Any | None = None, meta: ~typing.Dict[str, ~typing.Any] = <factory>)[source]#

Bases: object

Container for joint ERT-SRT inversion outputs.

chi2_ert: float | None = None#
chi2_srt: float | None = None#
ert_coverage: ndarray | None = None#
ert_log_resistivity: ndarray | None = None#
ert_predicted: ndarray | None = None#
ert_resistivity: ndarray | None = None#
iteration_history: list#
mesh: Any = None#
meta: Dict[str, Any]#
srt_coverage: ndarray | None = None#
srt_log_slowness: ndarray | None = None#
srt_predicted: ndarray | None = None#
srt_velocity: ndarray | None = None#
class PyHydroGeophysX.inversion.PetrophysicalCoupling[source]#

Bases: object

Coupling helpers from hydrological state to multi-method geophysics.

static compare_inversions_to_hydro(hydro_wc, ert_result, srt_result=None, em_result=None, petro_params=None)[source]#

Compare inversion products back to hydrological water content.

Returns misfit statistics for available methods.

static water_content_to_all_geophysics(water_content, porosity, **params)[source]#

Convert water content to resistivity, velocity, and conductivity.

Uses existing petrophysical functions in the package.

class PyHydroGeophysX.inversion.SRTInversion(data_file: str, mesh: pygimli.Mesh | None = None, **kwargs: Any)[source]#

Bases: InversionBase

Seismic Refraction Tomography inversion class.

Inverts travel-time data for subsurface velocity via log-slowness, and uses a Gauss-Newton optimization loop.

run(initial_model: ndarray | None = None) InversionResult[source]#

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() None[source]#

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

class PyHydroGeophysX.inversion.StructuralConstraint[source]#

Bases: object

Build structural constraint terms from one method for another.

static apply_structural_weights_to_Wm(Wm, boundary_weights)[source]#

Modify an existing smoothness matrix Wm to respect boundaries.

static build_cross_gradient_operator(mesh, model_a, model_b)[source]#

Build a cross-gradient operator surrogate.

Returns a diagonal sparse matrix weighted by local cross-gradient magnitude, where small values indicate better structural agreement.

static build_linearized_cross_gradient_blocks(RCM: ndarray, X: ndarray, model_a: ndarray, model_b: ndarray, mode: str = 'direct') Tuple[ndarray, ndarray][source]#

Build linearized cross-gradient blocks B1 and B2.

B1 multiplies model_a and B2 multiplies model_b. The resulting penalty terms are ||B1 m_a||^2 and ||B2 m_b||^2.

static build_local_design_matrix(mesh) ndarray[source]#

Build the local design matrix X=[x, z, 1] used by cross-gradient.

This mirrors the legacy xyzpos[:, 2] = 1 construction.

static build_neighborhood_matrix(mesh, Wm: Any | None = None, source: str = 'smoothness', correlation_lengths: Sequence[float] = (4.0, 4.0), threshold: float = 0.01, binarize: bool = True) ndarray[source]#

Build the neighborhood/correlation matrix used by cross-gradient.

Parameters:
  • mesh – Mesh object.

  • Wm – Smoothness matrix. Required when source='smoothness'.

  • source – One of 'smoothness' (legacy direct mode) or 'covariance' (legacy spatial mode).

  • correlation_lengths – Correlation lengths passed to pg.utils.covarianceMatrix when source='covariance'.

  • threshold – Entries with absolute value below this threshold are zeroed.

  • binarize – If True, convert nonzero entries to 1.

static from_conductivity_model(conductivity, mesh, gradient_threshold: float = 0.3)[source]#

Extract structural boundaries from an EM conductivity model.

static from_velocity_model(velocity_model, mesh, gradient_threshold: float = 0.3)[source]#

Extract structural boundaries from an SRT velocity model.

Returns a pairwise-boundary weight map that can reduce smoothing across detected structural interfaces.

class PyHydroGeophysX.inversion.TDEMInversion(times: ndarray, dobs: ndarray, uncertainties: ndarray, source_radius: float = 10.0, source_location: ndarray | None = None, receiver_location: ndarray | None = None, n_layers: int = 25, min_thickness: float = 1.0, max_thickness: float = 30.0, **kwargs)[source]#

Bases: object

Class for 1D TDEM inversion using SimPEG.

This class provides functionality for inverting TDEM sounding data to recover 1D layered Earth conductivity models.

Example

>>> # Load or create data
>>> times = np.logspace(-5, -2, 31)
>>> dobs = ...  # observed data
>>> uncertainties = ...  # data uncertainties
>>>
>>> # Create inversion
>>> inv = TDEMInversion(
...     times=times,
...     dobs=dobs,
...     uncertainties=uncertainties,
...     source_radius=10.0
... )
>>>
>>> # Run inversion
>>> result = inv.run()
plot_result(result: TDEMInversionResult, true_model: Tuple[ndarray, ndarray] | None = None, save_path: str | None = None) None[source]#

Plot inversion results.

Parameters:
  • result – TDEMInversionResult from run()

  • true_model – Tuple of (thicknesses, conductivity) for true model

  • save_path – Path to save figure (optional)

run(starting_model: ndarray | None = None) TDEMInversionResult[source]#

Run TDEM inversion.

Parameters:

starting_model – Initial log-conductivity model (optional)

Returns:

TDEMInversionResult with inversion results

setup() None[source]#

Set up inversion components (survey, mesh, simulation).

class PyHydroGeophysX.inversion.TDEMInversionResult(recovered_model: ~numpy.ndarray | None = None, recovered_conductivity: ~numpy.ndarray | None = None, l2_model: ~numpy.ndarray | None = None, l2_conductivity: ~numpy.ndarray | None = None, predicted_data: ~numpy.ndarray | None = None, mesh: discretize.TensorMesh | None = None, thicknesses: ~numpy.ndarray | None = None, chi2: float | None = None, iterations: int = 0, convergence_history: ~typing.List[float] = <factory>)[source]#

Bases: object

Container for TDEM inversion results.

recovered_model#

Final recovered model (log-conductivity)

Type:

numpy.ndarray

recovered_conductivity#

Recovered conductivity (S/m)

Type:

numpy.ndarray

l2_model#

L2 model before IRLS (log-conductivity)

Type:

numpy.ndarray

l2_conductivity#

L2 conductivity (S/m)

Type:

numpy.ndarray

predicted_data#

Predicted data from recovered model

Type:

numpy.ndarray

mesh#

TensorMesh used for inversion

Type:

discretize.TensorMesh

thicknesses#

Layer thicknesses used in inversion

Type:

numpy.ndarray

chi2#

Final chi-squared misfit

Type:

float

iterations#

Number of iterations

Type:

int

convergence_history#

History of data misfit per iteration

Type:

List[float]

chi2: float = None#
convergence_history: List[float]#
iterations: int = 0#
l2_conductivity: ndarray = None#
l2_model: ndarray = None#
mesh: discretize.TensorMesh = None#
predicted_data: ndarray = None#
recovered_conductivity: ndarray = None#
recovered_model: ndarray = None#
thicknesses: ndarray = None#
class PyHydroGeophysX.inversion.TimeLapseERTInversion(data_files: List[str], measurement_times: List[float], mesh: pygimli.Mesh | None = None, **kwargs)[source]#

Bases: InversionBase

Time-lapse ERT inversion class.

run(initial_model: ndarray | None = None) TimeLapseInversionResult[source]#

Run time-lapse ERT inversion.

Parameters:

initial_model – Initial model parameters (if None, a homogeneous model is used)

Returns:

TimeLapseInversionResult with inversion results

setup()[source]#

Set up time-lapse ERT inversion (load data, create operators, matrices, etc.)

class PyHydroGeophysX.inversion.TimeLapseInversionResult[source]#

Bases: InversionResult

Specialized class to store and manage results from time-lapse inversions. Inherits from InversionResult and adds attributes specific to time-lapse data.

final_models#

A 2D NumPy array where each column represents the inverted model for a specific timestep (shape: num_cells x num_timesteps).

Type:

Optional[np.ndarray]

timesteps#

An array or list of time values (e.g., hours, days) corresponding to each model slice in final_models.

Type:

Optional[np.ndarray]

all_coverage#

A list where each element is the coverage array for the corresponding timestep’s model.

Type:

List[np.ndarray]

all_chi2#

A list where each element is a list of chi-squared values per iteration for the inversion of that specific timestep or window. (Note: Original was List[float], if chi2 is per window, might need adjustment)

Type:

List[Any]

create_time_lapse_animation(output_filename: str, cmap: str = 'viridis', coverage_threshold: float | None = None, dpi: int = 100, fps: int = 2, **kwargs: Any) None[source]#

Create and save an animation (e.g., MP4 video) of the time-lapse inversion results.

Requires ffmpeg or another Matplotlib-supported animation writer to be installed.

Parameters:
  • output_filename (str) – The filename for the output animation (e.g., ‘timelapse_animation.mp4’).

  • cmap (str, optional) – Colormap for model values. Defaults to ‘viridis’.

  • coverage_threshold (Optional[float], optional) – Threshold for coverage masking. Defaults to None.

  • dpi (int, optional) – Dots Per Inch for the output animation. Defaults to 100.

  • fps (int, optional) – Frames Per Second for the animation. Defaults to 2.

  • **kwargs (Any) – Additional keyword arguments passed to plot_time_slice for each frame.

Raises:
  • ValueError – If final_models or mesh is missing.

  • ImportError – If matplotlib.animation cannot be imported.

classmethod load(filename: str) TimeLapseInversionResult[source]#

Load time-lapse inversion results. Overrides base class load.

plot_time_slice(timestep_idx: int, ax: Axes | None = None, cmap: str = 'viridis', coverage_threshold: float | None = None, **kwargs: Any) Tuple[Figure, Axes][source]#

Plot a single time slice (inverted model at a specific timestep) from the results.

Parameters:
  • timestep_idx (int) – The zero-based index of the timestep to plot.

  • ax (Optional[plt.Axes], optional) – Matplotlib Axes to plot on. If None, creates new.

  • cmap (str, optional) – Colormap for model values. Defaults to ‘viridis’.

  • coverage_threshold (Optional[float], optional) – Threshold for coverage masking. Defaults to None.

  • **kwargs (Any) – Additional arguments passed to pg.show.

Returns:

The matplotlib Figure and Axes objects.

Return type:

Tuple[plt.Figure, plt.Axes]

Raises:

ValueError – If models or mesh are missing, or if timestep_idx is out of range.

save(filename: str) None[source]#

Save time-lapse inversion results. Overrides base class save.

class PyHydroGeophysX.inversion.TimeLapseSRTInversion(data_files: List[str], measurement_times: List[float], mesh: pygimli.Mesh | None = None, **kwargs: Any)[source]#

Bases: InversionBase

Time-lapse Seismic Refraction Tomography inversion.

Architecture mirrors TimeLapseERTInversion while using TravelTimeManager and log-slowness model parameters.

run(initial_model: ndarray | None = None) TimeLapseInversionResult[source]#

Abstract method to run the main inversion loop.

Derived classes must implement this to perform the iterative optimization process. The method should populate and return an InversionResult (or subclass) object.

Returns:

An object containing the results of the inversion.

Return type:

InversionResult

setup() None[source]#

Abstract method to set up the inversion specifics.

This should be implemented by derived classes to prepare everything needed for the inversion, such as: - Creating or validating the mesh if not already done. - Initializing forward modeling operators. - Preparing data weighting and model regularization matrices. - Setting up initial models.

class PyHydroGeophysX.inversion.WindowedTimeLapseERTInversion(data_dir: str, ert_files: List[str], measurement_times: List[float], window_size: int = 3, mesh: pygimli.Mesh | str | None = None, **kwargs)[source]#

Bases: object

Class for windowed time-lapse ERT inversion to handle large temporal datasets.

run(window_parallel: bool = False, max_window_workers: int | None = None) TimeLapseInversionResult[source]#

Run windowed time-lapse ERT inversion.

Parameters:
  • window_parallel – Whether to process windows in parallel

  • max_window_workers – Maximum number of parallel workers (None for auto)

Returns:

TimeLapseInversionResult with stitched results

PyHydroGeophysX.inversion.run_tdem_inversion(times: ndarray, dobs: ndarray, uncertainties: ndarray, source_radius: float = 10.0, n_layers: int = 25, use_irls: bool = True, verbose: bool = True, **kwargs) TDEMInversionResult[source]#

Convenience function to run TDEM inversion.

Parameters:
  • times – Time channels (s)

  • dobs – Observed data (T)

  • uncertainties – Data uncertainties (T)

  • source_radius – Loop radius (m)

  • n_layers – Number of inversion layers

  • use_irls – Use IRLS for sparse inversion

  • verbose – Print progress

  • **kwargs – Additional parameters for TDEMInversion

Returns:

TDEMInversionResult with inversion results