PyHydroGeophysX package

Subpackages

Module contents

WatershedGeo - A comprehensive package for geophysical modeling and inversion in watershed monitoring.

This package integrates MODFLOW hydrological model outputs with geophysical forward modeling and inversion, specializing in electrical resistivity tomography (ERT) and seismic velocity modeling.

class PyHydroGeophysX.BaseVelocityModel[source]

Bases: object

Base class for seismic velocity models.

calculate_velocity(**kwargs) → ndarray[source]

Calculate seismic velocity from rock properties.

Parameters:

**kwargs – Rock properties specific to each model

Returns:

Seismic velocity values (Vp, Vs, or both)

class PyHydroGeophysX.BrieModel(exponent: float = 3.0)[source]

Bases: object

Brie’s model for calculating the effective bulk modulus of a partially saturated medium.

calculate_fluid_modulus(saturation: float, water_modulus: float = 2.0, gas_modulus: float = 0.01) → float[source]

Calculate effective fluid bulk modulus using Brie’s equation.

Parameters:

  • saturation – Water saturation (0 to 1)

  • water_modulus – Bulk modulus of water (GPa, default: 2.0)

  • gas_modulus – Bulk modulus of gas (GPa, default: 0.01)

Returns:

Effective fluid bulk modulus (GPa)

calculate_saturated_modulus(dry_modulus: float, mineral_modulus: float, porosity: float, saturation: float, water_modulus: float = 2.0, gas_modulus: float = 0.01) → float[source]

Calculate the saturated bulk modulus based on Brie’s equation.

Parameters:

  • dry_modulus – Bulk modulus of the dry rock (GPa)

  • mineral_modulus – Bulk modulus of the mineral matrix (GPa)

  • porosity – Porosity of the rock

  • saturation – Water saturation (0 to 1)

  • water_modulus – Bulk modulus of water (GPa, default: 2.0)

  • gas_modulus – Bulk modulus of gas (GPa, default: 0.01)

Returns:

Saturated bulk modulus (GPa)

class PyHydroGeophysX.CGLSSolver(max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.0, verbose=False)[source]

Bases: LinearSolver

CGLS (Conjugate Gradient Least Squares) solver.

class PyHydroGeophysX.DEMModel[source]

Bases: BaseVelocityModel

Differential Effective Medium (DEM) model for calculating elastic properties and seismic velocities of porous rocks.

calculate_velocity(porosity: ndarray, saturation: ndarray, bulk_modulus: float, shear_modulus: float, mineral_density: float, aspect_ratio: float = 0.1) → Tuple[ndarray, ndarray, ndarray][source]

Calculate P-wave velocity using the DEM model.

Parameters:

  • porosity – Porosity values (array)

  • saturation – Saturation values (array)

  • bulk_modulus – Initial bulk modulus of the solid matrix (GPa)

  • shear_modulus – Initial shear modulus of the solid matrix (GPa)

  • mineral_density – Density of the solid matrix (kg/m³)

  • aspect_ratio – Aspect ratio of pores (default: 0.1)

Returns:

Effective bulk modulus (GPa), effective shear modulus (GPa), and P-wave velocity (m/s)

class PyHydroGeophysX.ERTForwardModeling(mesh: pygimli.Mesh, data: pygimli.DataContainer | None = None)[source]

Bases: object

Class for forward modeling of Electrical Resistivity Tomography (ERT) data.

create_synthetic_data(xpos: ndarray, ypos: ndarray | None = None, mesh: pygimli.Mesh | None = None, res_models: ndarray | None = None, schemeName: str = 'wa', noise_level: float = 0.05, absolute_error: float = 0.0, relative_error: float = 0.05, save_path: str | None = None, show_data: bool = False, seed: int | None = None, xbound: float = 100, ybound: float = 100) → Tuple[pygimli.DataContainer, pygimli.Mesh][source]

Create synthetic ERT data using forward modeling.

This method simulates an ERT survey by placing electrodes, creating a measurement scheme, performing forward modeling to generate synthetic data, and adding noise.

Parameters:

  • xpos – X-coordinates of electrodes

  • ypos – Y-coordinates of electrodes (if None, uses flat surface)

  • mesh – Mesh for forward modeling

  • res_models – Resistivity model values

  • schemeName – Name of measurement scheme (‘wa’, ‘dd’, etc.)

  • noise_level – Level of Gaussian noise to add

  • absolute_error – Absolute error for data estimation

  • relative_error – Relative error for data estimation

  • save_path – Path to save synthetic data (if None, does not save)

  • show_data – Whether to display data after creation

  • seed – Random seed for noise generation

  • xbound – X boundary extension for mesh

  • ybound – Y boundary extension for mesh

Returns:

Tuple of (synthetic ERT data container, simulation mesh)

forward(resistivity_model: ndarray, log_transform: bool = True) → ndarray[source]

Compute forward response for a given resistivity model.

Parameters:

  • resistivity_model – Resistivity model values

  • log_transform – Whether resistivity_model is log-transformed

Returns:

Forward response (apparent resistivity)

forward_and_jacobian(resistivity_model: ndarray, log_transform: bool = True) → Tuple[ndarray, ndarray][source]

Compute forward response and Jacobian matrix.

Parameters:

  • resistivity_model – Resistivity model values

  • log_transform – Whether resistivity_model is log-transformed

Returns:

Tuple of (forward response, Jacobian matrix)

get_coverage(resistivity_model: ndarray, log_transform: bool = True) → ndarray[source]

Compute coverage (resolution) for a given resistivity model.

Parameters:

  • resistivity_model – Resistivity model values

  • log_transform – Whether resistivity_model is log-transformed

Returns:

Coverage values for each cell

set_data(data: pygimli.DataContainer) → None[source]

Set ERT data for forward modeling.

Parameters:

data – ERT data container

set_mesh(mesh: pygimli.Mesh) → None[source]

Set mesh for forward modeling.

Parameters:

mesh – PyGIMLI mesh

class PyHydroGeophysX.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.HertzMindlinModel(critical_porosity: float = 0.4, coordination_number: float = 4.0)[source]

Bases: BaseVelocityModel

Hertz-Mindlin model and Hashin-Shtrikman bounds for seismic velocity in porous rocks.

calculate_velocity(porosity: ndarray, saturation: ndarray, bulk_modulus: float, shear_modulus: float, mineral_density: float, depth: float = 1.0) → Tuple[ndarray, ndarray][source]

Calculate P-wave velocity for porous rocks.

Parameters:

  • porosity – Porosity values (array)

  • saturation – Saturation values (array)

  • bulk_modulus – Bulk modulus of the solid matrix (GPa)

  • shear_modulus – Shear modulus of the solid matrix (GPa)

  • mineral_density – Density of the solid matrix (kg/m³)

  • depth – Depth for pressure estimation (m, default: 1.0)

Returns:

Tuple of high-bound P-wave velocity (m/s) and low-bound P-wave velocity (m/s)

class PyHydroGeophysX.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.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.IterativeRefinement(max_iterations=5, tolerance=1e-10, use_double_precision=True)[source]

Bases: object

Iterative refinement to improve accuracy of a solution to a linear system.

refine(A, b, x0, solver_func)[source]

Perform iterative refinement.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • x0 – Initial solution

  • solver_func – Function that solves A*x = b

Returns:

Improved solution

class PyHydroGeophysX.LSQRSolver(max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.0, verbose=False)[source]

Bases: LinearSolver

LSQR solver for least squares problems.

class PyHydroGeophysX.LinearSolver(method='cgls', max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.0, verbose=False)[source]

Bases: object

Base class for linear system solvers.

solve(A, b, x0=None)[source]

Solve linear system Ax = b.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • x0 – Initial guess (None for zeros)

Returns:

Solution vector

class PyHydroGeophysX.MODFLOWWaterContent(model_directory: str, idomain: ndarray)[source]

Bases: HydroModelOutput

Class for processing water content data from MODFLOW simulations.

get_timestep_info() → List[Tuple[int, int, float, float]][source]

Get information about each timestep in the WaterContent file.

Returns:

List of tuples (kstp, kper, pertim, totim) for each timestep

load_time_range(start_idx: int = 0, end_idx: int | None = None, nlay: int = 3) → ndarray[source]

Load water content for a range of timesteps.

Parameters:

  • start_idx – Starting timestep index (default: 0)

  • end_idx – Ending timestep index (exclusive, default: None loads all)

  • nlay – Number of layers in the model (default: 3)

Returns:

Water content array with shape (timesteps, nlay, nrows, ncols)

load_timestep(timestep_idx: int, nlay: int = 3) → ndarray[source]

Load water content for a specific timestep.

Parameters:

  • timestep_idx – Index of the timestep to load

  • nlay – Number of layers in the model

Returns:

Water content array with shape (nlay, nrows, ncols)

class PyHydroGeophysX.MeshCreator(quality: float = 28, area: float = 40)[source]

Bases: object

Class for creating and managing meshes for geophysical inversion.

create_from_ert_data(data, max_depth: float = 30.0, quality: float = 34)[source]

Create a mesh suitable for ERT inversion from ERT data.

Parameters:

  • data – PyGIMLI ERT data object

  • max_depth – Maximum depth of the mesh

  • quality – Mesh quality parameter

Returns:

PyGIMLI mesh for ERT inversion

create_from_layers(surface: ndarray, layers: List[ndarray], bottom_depth: float = 30.0, markers: List[int] | None = None) → pygimli.Mesh[source]

Create a mesh from surface and layer boundaries.

Parameters:

  • surface – Surface coordinates [[x,z],…]

  • layers – List of layer boundary coordinates

  • bottom_depth – Depth below surface minimum for mesh bottom

  • markers – List of markers for each layer (default: [2, 3, 2, …])

Returns:

PyGIMLI mesh

class PyHydroGeophysX.ProfileInterpolator(point1: List[int], point2: List[int], surface_data: ndarray, origin_x: float = 0.0, origin_y: float = 0.0, pixel_width: float = 1.0, pixel_height: float = -1.0, num_points: int = 200)[source]

Bases: object

Class for handling interpolation of data to/from profiles.

interpolate_3d_data(data: ndarray) → ndarray[source]

Interpolate 3D data (n_layers, ny, nx) to profile.

Parameters:

data – 3D array of values

Returns:

Array of interpolated values (n_layers, n_profile_points)

interpolate_layer_data(layer_data: List[ndarray]) → ndarray[source]

Interpolate multiple layer data to profile.

Parameters:

layer_data – List of 2D arrays for each layer

Returns:

Array of interpolated values (n_layers, n_profile_points)

interpolate_to_mesh(property_values: ndarray, depth_values: ndarray, mesh_x: ndarray, mesh_y: ndarray, mesh_markers: ndarray, ID: ndarray, layer_markers: list = [3, 0, 2]) → ndarray[source]

Interpolate property values from profile to mesh with layer-specific handling.

Parameters:

  • property_values – Property values array (n_points or n_layers, n_points)

  • depth_values – Depth values array (n_layers, n_points)

  • mesh_x – Coordinates of mesh cells

  • mesh_y – Coordinates of mesh cells

  • mesh_markers – Markers indicating different layers in mesh

  • layer_markers – List of marker values for each layer

Returns:

Interpolated values for mesh cells

class PyHydroGeophysX.RRLSQRSolver(max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.1, verbose=False)[source]

Bases: LinearSolver

Regularized LSQR solver.

class PyHydroGeophysX.RRLSSolver(max_iterations=200, tolerance=1e-08, use_gpu=False, parallel=False, n_jobs=-1, damping=0.0, verbose=False)[source]

Bases: LinearSolver

Range-Restricted Least Squares solver.

class PyHydroGeophysX.TikhonvRegularization(regularization_matrix=None, alpha=1.0, regularization_type='identity')[source]

Bases: object

Tikhonov regularization for ill-posed problems.

apply(A, b, solver=None)[source]

Apply Tikhonov regularization to the linear system.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • solver – Solver to use (None for direct solver)

Returns:

Regularized solution

create_regularization_matrix(n)[source]

Create regularization matrix based on the selected type.

Parameters:

n – Size of model vector

Returns:

Regularization matrix

class PyHydroGeophysX.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.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.VRHModel[source]

Bases: BaseVelocityModel

Voigt-Reuss-Hill (VRH) mixing model for effective elastic properties of composites.

calculate_properties(fractions: List[float], bulk_moduli: List[float], shear_moduli: List[float], densities: List[float]) → Tuple[float, float, float][source]

Calculate effective elastic properties using the VRH model.

Parameters:

  • fractions – Volume fractions of each mineral (must sum to 1)

  • bulk_moduli – Bulk moduli of each mineral (GPa)

  • shear_moduli – Shear moduli of each mineral (GPa)

  • densities – Densities of each mineral (kg/m³)

Returns:

Effective bulk modulus (GPa), effective shear modulus (GPa), and effective density (kg/m³)

calculate_velocity(fractions: List[float], bulk_moduli: List[float], shear_moduli: List[float], densities: List[float]) → Tuple[float, float][source]

Calculate P-wave and S-wave velocities using the VRH model.

Parameters:

  • fractions – Volume fractions of each mineral (must sum to 1)

  • bulk_moduli – Bulk moduli of each mineral (GPa)

  • shear_moduli – Shear moduli of each mineral (GPa)

  • densities – Densities of each mineral (kg/m³)

Returns:

P-wave velocity (m/s) and S-wave velocity (m/s)

PyHydroGeophysX.VRH_model(f=[0.35, 0.25, 0.2, 0.125, 0.075], K=[55.4, 36.6, 75.6, 46.7, 50.4], G=[28.1, 45, 25.6, 23.65, 27.4], rho=[2560, 2650, 2630, 2540, 3050])[source]

Implements the Voigt-Reuss-Hill (VRH) mixing model to estimate the effective bulk modulus (Km), shear modulus (Gm), and density (rho_b) of a composite material made from various minerals.

Parameters: f (list): Fraction of each mineral in the composite (must sum to 1). K (list): Bulk modulus of each mineral (GPa). G (list): Shear modulus of each mineral (GPa). rho (list): Density of each mineral (kg/m^3).

Returns: Km (float): Effective bulk modulus of the composite material (GPa). Gm (float): Effective shear modulus of the composite material (GPa). rho_b (float): Effective density of the composite material (kg/m^3).

class PyHydroGeophysX.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.create_mesh_from_layers(surface: ndarray, line1: ndarray, line2: ndarray, bottom_depth: float = 30.0, quality: float = 28, area: float = 40) → Tuple[pygimli.Mesh, ndarray, ndarray][source]

Create mesh from layer boundaries and get cell centers and markers.

Parameters:

  • surface – Surface coordinates [[x,z],…]

  • line1 – First layer boundary coordinates

  • line2 – Second layer boundary coordinates

  • bottom_depth – Depth below surface minimum for mesh bottom

  • quality – Mesh quality parameter

  • area – Maximum cell area

Returns:

PyGIMLI mesh mesh_centers: Array of cell center coordinates markers: Array of cell markers

Return type:

mesh

PyHydroGeophysX.direct_solver(A, b, method='lu', **kwargs)[source]

Solve a linear system using direct methods.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • method – Direct solver method (‘lu’, ‘qr’, ‘svd’, ‘cholesky’)

  • **kwargs – Additional parameters for specific methods

Returns:

Solution vector

PyHydroGeophysX.ertforandjac(fob, rhomodel, xr)[source]

Forward model and Jacobian for ERT.

Parameters:

  • fob (pygimli.ERTModelling) – ERT forward operator.

  • rhomodel (pg.RVector) – Resistivity model.

  • xr (np.ndarray) – Log-transformed model parameter.

Returns:

Log-transformed forward response. J (np.ndarray): Jacobian matrix.

Return type:

dr (np.ndarray)

PyHydroGeophysX.ertforandjac2(fob, xr, mesh)[source]

Alternative ERT forward model and Jacobian using log-resistivity values.

Parameters:

  • fob (pygimli.ERTModelling) – ERT forward operator.

  • xr (np.ndarray) – Log-transformed model parameter.

  • mesh (pg.Mesh) – Mesh for the forward model.

Returns:

Log-transformed forward response. J (np.ndarray): Jacobian matrix.

Return type:

dr (np.ndarray)

PyHydroGeophysX.ertforward(fob, mesh, rhomodel, xr)[source]

Forward model for ERT.

Parameters:

  • fob (pygimli.ERTModelling) – ERT forward operator.

  • mesh (pg.Mesh) – Mesh for the forward model.

  • rhomodel (pg.RVector) – Resistivity model vector.

  • xr (np.ndarray) – Log-transformed model parameter (resistivity).

Returns:

Log-transformed forward response. rhomodel (pg.RVector): Updated resistivity model.

Return type:

dr (np.ndarray)

PyHydroGeophysX.ertforward2(fob, xr, mesh)[source]

Simplified ERT forward model.

Parameters:

  • fob (pygimli.ERTModelling) – ERT forward operator.

  • xr (np.ndarray) – Log-transformed model parameter.

  • mesh (pg.Mesh) – Mesh for the forward model.

Returns:

Log-transformed forward response.

Return type:

dr (np.ndarray)

PyHydroGeophysX.generalized_solver(A, b, method='cgls', x=None, maxiter=200, tol=1e-08, verbose=False, damp=0.0, use_gpu=False, parallel=False, n_jobs=-1)[source]

Generalized solver for Ax = b with optional GPU acceleration and parallelism.

Parameters:

Aarray_like or sparse matrix

The system matrix (Jacobian or forward operator).

barray_like

Right-hand side vector.

methodstr, optional

Solver method: ‘lsqr’, ‘rrlsqr’, ‘cgls’, or ‘rrls’. Default is ‘cgls’.

xarray_like, optional

Initial guess for the solution. If None, zeros are used.

maxiterint, optional

Maximum number of iterations.

tolfloat, optional

Convergence tolerance.

verbosebool, optional

Print progress information every 10 iterations.

dampfloat, optional

Damping factor (Tikhonov regularization).

use_gpubool, optional

Use GPU acceleration with CuPy (if available).

parallelbool, optional

Use parallel CPU computations.

n_jobsint, optional

Number of parallel jobs (if parallel is True).

Returns:

xarray_like

The computed solution vector.

PyHydroGeophysX.get_optimal_solver(A, b, estimate_condition=True, time_limit=None, memory_limit=None)[source]

Automatically select the optimal solver for a given linear system.

Parameters:

  • A – System matrix

  • b – Right-hand side vector

  • estimate_condition – Whether to estimate condition number

  • time_limit – Maximum allowed solution time (seconds)

  • memory_limit – Maximum allowed memory usage (bytes)

Returns:

Tuple of (solver_object, solver_info)

PyHydroGeophysX.interpolate_structure_to_profile(structure_data: List[ndarray], X_grid: ndarray, Y_grid: ndarray, X_pro: ndarray, Y_pro: ndarray) → ndarray[source]

Interpolate multiple structure layers onto profile

Parameters:

  • structure_data – List of 2D arrays for each layer

  • X_grid – X coordinates of original grid

  • Y_grid – Y coordinates of original grid

  • X_pro – X coordinates of profile points

  • Y_pro – Y coordinates of profile points

Returns:

Array of interpolated values with shape (n_layers, n_points)

PyHydroGeophysX.interpolate_to_mesh(property_values: ndarray, profile_distance: ndarray, depth_values: ndarray, mesh_x: ndarray, mesh_y: ndarray, mesh_markers: ndarray, ID, layer_markers: list = [3, 0, 2]) → ndarray[source]

Interpolate property values from profile to mesh with layer-specific handling.

Parameters:

  • property_values – Property values array (n_points)

  • profile_distance – Distance along profile (n_points)

  • depth_values – Depth values array (n_layers, n_points)

  • mesh_x – X coordinates of mesh cells

  • mesh_y – Y coordinates of mesh cells

  • mesh_markers – Markers indicating different layers in mesh

  • layer_markers – List of marker values for each layer

Returns:

Interpolated values for mesh cells

PyHydroGeophysX.interpolate_to_profile(data: ndarray, X_grid: ndarray, Y_grid: ndarray, X_pro: ndarray, Y_pro: ndarray, method: str = 'linear') → ndarray[source]

Interpolate 2D data onto a profile line

Parameters:

  • data – 2D array of values to interpolate

  • X_grid – X coordinates of original grid (meshgrid)

  • Y_grid – Y coordinates of original grid (meshgrid)

  • X_pro – X coordinates of profile points

  • Y_pro – Y coordinates of profile points

  • method – Interpolation method (‘linear’ or ‘nearest’)

Returns:

Interpolated values along profile

PyHydroGeophysX.prepare_2D_profile_data(data: ndarray, XX: ndarray, YY: ndarray, X_pro: ndarray, Y_pro: ndarray) → ndarray[source]

Interpolate multiple 2D gridded data layers onto a profile line.

Parameters:

  • data – 3D array of gridded data (n_layers, ny, nx)

  • XX – Coordinate grids from meshgrid

  • YY – Coordinate grids from meshgrid

  • X_pro – Profile line coordinates

  • Y_pro – Profile line coordinates

Returns:

Interpolated values along profile (n_layers, n_profile_points)

PyHydroGeophysX.satK(Keff, Km, phi, Sat)[source]

Calculate the saturated bulk modulus (K_sat) based on Brie’s equation.

Parameters: Keff (float): Effective bulk modulus of the dry rock (GPa). Km (float): Bulk modulus of the matrix (GPa). phi (float): Porosity of the rock. Sat (float): Saturation level of the fluid in the pores.

Returns: float: Saturated bulk modulus (GPa).

PyHydroGeophysX.setup_profile_coordinates(point1: List[int], point2: List[int], surface_data: ndarray, origin_x: float = 0.0, origin_y: float = 0.0, pixel_width: float = 1.0, pixel_height: float = -1.0, num_points: int = 200) → Tuple[ndarray, ndarray, ndarray, ndarray, ndarray][source]

Set up profile coordinates based on surface elevation data between two points

Parameters:

  • point1 – Starting point indices [col, row]

  • point2 – Ending point indices [col, row]

  • surface_data – 2D array of surface elevation data

  • origin_x – X coordinate of origin

  • origin_y – Y coordinate of origin

  • pixel_width – Width of each pixel

  • pixel_height – Height of each pixel (negative for top-down)

  • num_points – Number of points along profile

Returns:

X coordinates along profile Y_pro: Y coordinates along profile L_profile: Distances along profile XX: X coordinate grid YY: Y coordinate grid

Return type:

X_pro

PyHydroGeophysX.velDEM(phi, Km, Gm, rho_b, Sat, alpha)[source]

Calculate effective bulk modulus (Keff), shear modulus (Geff), and P-wave velocity (Vp) for a rock with varying porosity (phi) based on the DEM model, taking into account the saturation (Sat) and the crack aspect ratio (alpha).

Parameters: phi (np.array): Array of porosities. Km (float): Initial bulk modulus of the material (GPa). Gm (float): Initial shear modulus of the material (GPa). rho_b (float): Density of the solid phase (kg/m^3). Sat (float): Saturation level of the fluid in the cracks (0 to 1, where 1 is fully saturated). alpha (float): Crack aspect ratio.

Returns: Keff1 (np.array): Effective bulk modulus for each porosity value (GPa). Geff1 (np.array): Effective shear modulus for each porosity value (GPa). Vp (np.array): P-wave velocity for each porosity value (m/s).

PyHydroGeophysX.vel_porous(phi, Km, Gm, rho_b, Sat, depth=1)[source]

Calculate P-wave velocity (Vp) for a rock with varying porosity (phi) based on the Hertz-Mindlin model and Hashin-Shtrikman bounds, taking into account the saturation (Sat).

Parameters: phi (np.array): Array of porosities. Km (float): Bulk modulus of the solid phase (GPa). Gm (float): Shear modulus of the solid phase (GPa). rho_b (float): Density of the solid phase (kg/m^3). Sat (float): Saturation level of the fluid in the pores (0 to 1, where 1 is fully saturated). depth (float): depth for pressure estimation (m)

Returns: Vp_h (np.array): P-wave velocity for each porosity value (upper Hashin-Shtrikman bound) (m/s). Vp_l (np.array): P-wave velocity for each porosity value (lower Hashin-Shtrikman bound) (m/s).