pycsamt.forward.grid3d#

3-D resistivity grid for MT forward modelling.

Grid3D stores the subsurface resistivity as a regular (non-uniform) finite-difference grid in x–y–z together with the surface station positions needed by MT3DForward.

Coordinate convention#

  • x — easting [m], increasing to the right.

  • y — northing [m], increasing into the page.

  • z — depth [m], increasing downward.

  • resistivity[iz, iy, ix] is the resistivity of the cell whose top-left-front corner is at (x_nodes[ix], y_nodes[iy], z_nodes[iz]).

Padding#

The quasi-3D solver extracts 2-D XZ and YZ slices and runs MT2DForward on each. Those 2-D solvers need padding to push Dirichlet BCs away from the model region. n_pad records the padding count added on each side in x and y and at the bottom in z, consistent with Grid2D.

Classes

Grid3D(dx, dy, dz, resistivity, stations_xy)

Non-uniform 3-D finite-difference grid for MT forward modelling.
class pycsamt.forward.grid3d.Grid3D(dx, dy, dz, resistivity, stations_xy, n_pad=0, name='')[source]#

Bases: object

Non-uniform 3-D finite-difference grid for MT forward modelling.

Parameters:

  • dx (array-like, shape (nx,)) – Cell widths in x (easting) [m], including padding.

  • dy (array-like, shape (ny,)) – Cell widths in y (northing) [m], including padding.

  • dz (array-like, shape (nz,)) – Cell heights in z (depth) [m], including padding.

  • resistivity (array-like, shape (nz, ny, nx)) – Cell resistivities [Ω·m], top→bottom, front→back, left→right.

  • stations_xy (array-like, shape (n_stations, 2)) – Surface station (x, y) coordinates [m]. Must lie within the grid extent.

  • n_pad (int) – Padding cells added on each side in x and y and at the bottom in z during construction.

  • name (str) – Optional label.

Examples

Uniform halfspace:

>>> g = Grid3D.halfspace(rho=100.0, nx=20, ny=20, nz=15,
...     x_max=8_000.0, y_max=8_000.0, z_max=4_000.0,
...     nx_stations=5, ny_stations=5)
>>> g.resistivity.shape   # includes padding
(23, 28, 28)

3-D conductive block:

>>> g = Grid3D.block_anomaly(
...     bg_rho=500.0, anomaly_rho=5.0,
...     bounds=(2_000.0, 6_000.0, 2_000.0, 6_000.0, 300.0, 1_500.0),
...     nx=25, ny=25, nz=18,
...     x_max=8_000.0, y_max=8_000.0, z_max=5_000.0,
...     nx_stations=5, ny_stations=5)
dx: ndarray#
dy: ndarray#
dz: ndarray#
resistivity: ndarray#
stations_xy: ndarray#
n_pad: int = 0#
name: str = ''#
property nx: int[source]#
property ny: int[source]#
property nz: int[source]#
property n_stations: int[source]#
property x_nodes: ndarray[source]#
property y_nodes: ndarray[source]#
property z_nodes: ndarray[source]#
property x_centers: ndarray[source]#
property y_centers: ndarray[source]#
property z_centers: ndarray[source]#
property x_extent: float[source]#
property y_extent: float[source]#
property z_extent: float[source]#
xz_slice(yi)[source]#

Extract an XZ (east-west) 2-D slice at y-cell yi.

Returns a Grid2D whose horizontal axis is x (easting) and vertical axis is z (depth). The station x-positions are those of all stations whose y-coordinate falls in cell yi.

Returns:

  • grid2d (Grid2D)

  • station_indices (ndarray of int) – Global indices (into stations_xy) of stations in this y-row. grid2d.x_stations[k] is the x-coordinate of stations_xy[station_indices[k]].

Parameters:

yi (int)

Return type:

tuple[Grid2D, ndarray]

yz_slice(xi)[source]#

Extract a YZ (north-south) 2-D slice at x-cell xi.

The Grid2D horizontal axis is mapped to y (northing); its dx array contains dy and its resistivity is resistivity[:, :, xi].

Returns:

  • grid2d (Grid2D)

  • station_indices (ndarray of int) – Global indices of stations in this x-column. grid2d.x_stations[k] is the y-coordinate of the station (used as the horizontal axis in the yz-plane solver).

Parameters:

xi (int)

Return type:

tuple[Grid2D, ndarray]

column_profile_3d(xi, yi)[source]#

Return the 1-D resistivity/thickness profile at cell (xi, yi).

Parameters:

  • xi (int Column index (x direction).)

  • yi (int Row index (y direction).)

Returns:

  • rho (ndarray, shape (nz,))

  • thick (ndarray, shape (nz-1,))

Return type:

tuple[ndarray, ndarray]

property conductivity: ndarray[source]#

Cell conductivities [S/m], shape (nz, ny, nx).

plot(*, cmap='jet_r', log_scale=True, clip_core=True, show_stations=True, figsize=(13, 4.5), vmin=None, vmax=None)[source]#

Plot orthogonal slices: XZ (mid-y), YZ (mid-x), XY (mid-z).

Parameters:
Returns:

  • fig (matplotlib.figure.Figure)

  • axes (ndarray of Axes, shape (3,))

classmethod halfspace(rho=100.0, *, nx=20, ny=20, nz=15, x_max=8000.0, y_max=8000.0, z_max=4000.0, n_pad=8, pad_factor=1.3, nx_stations=5, ny_stations=5, name='halfspace')[source]#

Create a uniform resistivity 3-D halfspace with a regular station grid.

Parameters:

  • rho (float) – Background resistivity [Ω·m].

  • nx (int) – Core cell counts (padding added automatically).

  • ny (int) – Core cell counts (padding added automatically).

  • nz (int) – Core cell counts (padding added automatically).

  • x_max (float) – Core horizontal extents [m].

  • y_max (float) – Core horizontal extents [m].

  • z_max (float) – Core depth extent [m].

  • n_pad (int) – Padding cells per side (x, y) and at the bottom (z).

  • pad_factor (float) – Exponential growth factor for padding cells.

  • nx_stations (int) – Regular station grid dimensions.

  • ny_stations (int) – Regular station grid dimensions.

  • name (str)

Return type:

Grid3D

classmethod block_anomaly(bg_rho=100.0, anomaly_rho=5.0, bounds=(2000.0, 6000.0, 2000.0, 6000.0, 300.0, 1500.0), *, nx=25, ny=25, nz=18, x_max=8000.0, y_max=8000.0, z_max=5000.0, n_pad=8, pad_factor=1.3, nx_stations=5, ny_stations=5, name='')[source]#

Create a background halfspace with one rectangular 3-D anomaly.

Parameters:

  • bg_rho (float) – Background resistivity [Ω·m].

  • anomaly_rho (float) – Anomaly resistivity [Ω·m].

  • bounds ((x_lo, x_hi, y_lo, y_hi, z_lo, z_hi)) – Anomaly extents in core coordinates [m].

  • nx (int Core cell counts.)

  • ny (int Core cell counts.)

  • nz (int Core cell counts.)

  • x_max (float Core extents [m].)

  • y_max (float Core extents [m].)

  • z_max (float Core extents [m].)

  • n_pad (int, float)

  • pad_factor (int, float)

  • nx_stations (int)

  • ny_stations (int)

  • name (str)

Return type:

Grid3D

Examples

3-D conductive fault zone:

>>> g = Grid3D.block_anomaly(
...     bg_rho=500.0, anomaly_rho=3.0,
...     bounds=(2_500., 5_500., 2_500., 5_500., 200., 1_800.))
classmethod random_layered(*, n_layers=4, nx=20, ny=20, nz=15, x_max=8000.0, y_max=8000.0, z_max=4000.0, n_pad=8, pad_factor=1.3, nx_stations=5, ny_stations=5, rho_min=1.0, rho_max=10000.0, lateral_variation=True, corr_length=2000.0, seed=None, name='random_3d')[source]#

Generate a random layered 3-D model with optional lateral variation.

Resistivities for n_layers horizontal layers are drawn from a log-uniform prior. When lateral_variation is True, each column’s resistivity is smoothly perturbed using a Gaussian random field with correlation length corr_length.

Parameters:

  • n_layers (int Number of horizontal layers.)

  • nx (as above.)

  • ny (as above.)

  • nz (as above.)

  • x_max (as above.)

  • y_max (as above.)

  • z_max (as above.)

  • n_pad (as above.)

  • pad_factor (as above.)

  • rho_min (float Resistivity bounds [Ω·m].)

  • rho_max (float Resistivity bounds [Ω·m].)

  • lateral_variation (bool)

  • corr_length (float GRF correlation length [m].)

  • seed (int or Generator or None.)

  • name (str)

  • nx_stations (int)

  • ny_stations (int)

Return type:

Grid3D