2.20.3.7. pycsamt.forward.grid2d#

2-D resistivity grid for MT forward modelling.

The Grid2D dataclass is the central data object shared by the 2-D MT forward solver and any downstream processing that needs the subsurface resistivity structure. It stores the model as a regular (but non-uniform) finite-difference grid together with station positions and padding metadata.

Coordinate convention#

  • x — horizontal distance along the profile [m], increasing to the right.

  • z — depth below surface [m], increasing downward.

  • Cell (i, j) has its top-left corner at (x_nodes[j], z_nodes[i]) where i is the row index (depth) and j is the column index (x).

  • resistivity[i, j] is the resistivity of cell (i, j) [Ω·m].

Grid layout (nx = 4, nz = 3, no padding shown):

x_nodes:   0    dx[0]  dx[0]+dx[1]  …
           ┌────┬────┬────┬────┐
z_nodes[0] │(0,0)│(0,1)│(0,2)│(0,3)│  ← row 0 (shallowest)
           ├────┼────┼────┼────┤
z_nodes[1] │(1,0)│(1,1)│(1,2)│(1,3)│
           ├────┼────┼────┼────┤
z_nodes[2] │(2,0)│(2,1)│(2,2)│(2,3)│  ← row nz-1 (deepest)
           └────┴────┴────┴────┘

Padding#

The FD solver requires padding cells on both horizontal sides and at the bottom to push the Dirichlet boundaries far enough from the model region that edge effects are negligible. Padding cells are constructed with exponentially growing widths/depths. The resistivity in padding cells is copied from the nearest edge column (horizontal) or deepest row (vertical).

The n_pad attribute records how many padding cells were added on each side so that station positions and visualisation can clip to the core region (core_slice).

References

Wannamaker, P.E., Stodt, J.A., Rijo, L. (1987). A stable finite element solution for two-dimensional magnetotelluric modelling. Geophys. J. Int., 88(1), 277-296.

de Groot-Hedlin, C., Constable, S. (1990). Occam’s inversion to generate smooth, two-dimensional models from magnetotelluric data. Geophysics, 55(12), 1613-1624.

Functions

make_padding(core_spacing, n_pad[, factor])

Return cell widths for one padding strip.

Classes

Grid2D(dx, dz, resistivity, x_stations[, ...])

Non-uniform 2-D finite-difference grid for MT forward modelling.
class pycsamt.forward.grid2d.Grid2D(dx, dz, resistivity, x_stations, n_pad=0, name='')[source]

Bases: object

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

Parameters:

  • dx (array-like, shape (nx,)) – Cell widths in the horizontal (x) direction [m].

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

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

  • x_stations (array-like, shape (n_stations,)) – Horizontal positions of MT stations [m]. Must lie within [x_nodes[0], x_nodes[-1]].

  • n_pad (int) – Number of padding cells on each side (left, right) and at the bottom that were added during construction. Used internally to identify the core model region. Set to 0 if no padding is present.

  • name (str) – Optional label for the model.

Variables:

  • nx (int) – Total number of cells in x (includes padding).

  • nz (int) – Total number of cells in z (includes padding).

  • x_nodes (ndarray, shape (nx+1,)) – x-coordinates of vertical node lines [m].

  • z_nodes (ndarray, shape (nz+1,)) – z-coordinates of horizontal node lines (depth) [m].

  • x_centers (ndarray, shape (nx,)) – x-coordinates of cell centres [m].

  • z_centers (ndarray, shape (nz,)) – z-coordinates (depth) of cell centres [m].

  • core_x_slice (slice) – Slice that selects core (non-padding) columns from any [:, j] array.

  • core_z_slice (slice) – Slice that selects core (non-padding) rows from any [i, :] array.

Examples

Uniform halfspace:

>>> g = Grid2D.halfspace(rho=100.0, nx=30, nz=20,
...                       x_max=5_000.0, z_max=3_000.0,
...                       n_stations=10, station_x_max=4_000.0)
>>> g.resistivity.shape
(20, 30)

One conductive anomaly:

>>> g = Grid2D.with_anomaly(
...     bg_rho=500.0,
...     anomaly_rho=5.0,
...     anomaly_bounds=(1_000.0, 3_000.0, 200.0, 800.0),
...     nx=40, nz=25,
...     x_max=6_000.0, z_max=4_000.0,
...     n_stations=12,
... )

1-D layered model extended horizontally:

>>> from pycsamt.forward.synthetic import LayeredModel
>>> m = LayeredModel([100, 10, 500], [300, 800])
>>> g = Grid2D.from_1d_layers(m, nx=40, x_max=6_000.0, n_stations=12)
dx: ndarray
dz: ndarray
resistivity: ndarray
x_stations: ndarray
n_pad: int = 0
name: str = ''
property nx: int[source]

Number of cells in x.

property nz: int[source]

Number of cells in z (depth).

property n_stations: int[source]

Number of MT stations.

property n_nodes: int[source]

Total FD node count (nx+1) × (nz+1).

property x_nodes: ndarray[source]

x-coordinates of vertical node lines [m], shape (nx+1,).

property z_nodes: ndarray[source]

Depth coordinates of horizontal node lines [m], shape (nz+1,).

property x_centers: ndarray[source]

x-coordinates of cell centres [m], shape (nx,).

property z_centers: ndarray[source]

Depth coordinates of cell centres [m], shape (nz,).

property x_extent: float[source]

Total horizontal extent of the grid [m].

property z_extent: float[source]

Total depth extent of the grid [m].

property core_x_slice: slice[source]

Column slice that selects only core (non-padding) cells.

property core_z_slice: slice[source]

Row slice that selects only core (non-padding) cells.

property core_resistivity: ndarray[source]

Resistivity array clipped to the core (non-padding) region.

property core_x_offset: float[source]

x-coordinate where the core (non-padding) region begins [m].

Padding is built identically on both sides of the grid (the same cell-width list, left side reversed), so this is also where along-profile chainage – 0 at the first core column – and the grid’s own absolute x coordinate (0 at the outer edge of the left padding) part ways. See value_at().

value_at(x_m, z_m, *, chainage=False)[source]

Return the resistivity of the cell nearest (x_m, z_m).

Parameters:

  • x_m (float) – Horizontal position, metres. In the grid’s own absolute coordinates (0 at the outer edge of the left padding) by default.

  • z_m (float) – Depth, metres, positive downward.

  • chainage (bool, default False) – When True, x_m is along-profile chainage instead – 0 at the first core (station-carrying) column, matching pycsamt.geology.Borehole’s x convention. Set this when sampling ground truth for a borehole or any other profile-relative position.

Return type:

float

Examples

>>> g = Grid2D.halfspace(rho=100.0, nx=20, x_max=2000.0, n_stations=5)
>>> g.value_at(0.0, 50.0, chainage=True)
100.0
column_profile(col)[source]

Return the 1-D resistivity/thickness profile of column col.

Used by the FD solver to compute Dirichlet boundary conditions at the left and right edges from the 1-D MT response of the edge columns.

Parameters:

col (int) – Column index (0 = leftmost, nx-1 = rightmost).

Returns:

  • rho (ndarray, shape (nz,)) – Layer resistivities top-to-bottom [Ω·m].

  • thick (ndarray, shape (nz-1,)) – Layer thicknesses [m]. (Bottom layer is a halfspace.)

Return type:

tuple[ndarray, ndarray]

station_cell_indices()[source]

Return the column index of the cell that contains each station.

Returns:

indices

Return type:

ndarray of int, shape (n_stations,)

station_node_indices()[source]

Return the index of the nearest surface node for each station.

Surface nodes are at z = 0 and numbered 0, 1, …, nx.

Returns:

indices

Return type:

ndarray of int, shape (n_stations,)

property conductivity: ndarray[source]

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

harmonic_mean_x()[source]

Harmonic-mean conductivity at x-directed cell interfaces.

Used by the TM-mode FD assembler.

Returns:

sigma_hx[i, j] is the harmonic mean of sigma[i, j] and sigma[i, j+1].

Return type:

ndarray, shape (nz, nx-1)

harmonic_mean_z()[source]

Harmonic-mean conductivity at z-directed cell interfaces.

Returns:

sigma_hz[i, j] is the harmonic mean of sigma[i, j] and sigma[i+1, j].

Return type:

ndarray, shape (nz-1, nx)

plot(ax=None, *, log_scale=True, cmap='jet_r', clip_core=True, show_stations=True, figsize=(10, 5), vmin=None, vmax=None)[source]

Plot the 2-D resistivity model.

Parameters:

  • ax (Axes or None) – Target axes. Created when not given.

  • log_scale (bool) – Plot log₁₀(ρ) rather than ρ.

  • cmap (str) – Matplotlib colormap name.

  • clip_core (bool) – Clip to the non-padding core region.

  • show_stations (bool) – Mark station positions on the surface.

  • figsize ((float, float)) – Figure size when ax is not provided.

  • vmin (float or None) – Colormap limits in log₁₀(Ω·m) space when log_scale is True.

  • vmax (float or None) – Colormap limits in log₁₀(Ω·m) space when log_scale is True.

Returns:

ax

Return type:

Axes

classmethod halfspace(rho=100.0, *, nx=40, nz=30, x_max=10000.0, z_max=5000.0, n_pad=8, pad_factor=1.3, n_stations=10, station_x_max=None, name='halfspace')[source]

Create a uniform resistivity halfspace grid.

Parameters:

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

  • nx (int) – Number of core cells in x (padding added on both sides).

  • nz (int) – Number of core cells in z (padding added at the bottom).

  • x_max (float) – Total horizontal extent of the core region [m].

  • z_max (float) – Total depth of the core region [m].

  • n_pad (int) – Number of padding cells on each side (left, right, bottom).

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

  • n_stations (int) – Number of equally spaced stations across the core region.

  • station_x_max (float or None) – Right-most station position [m]. Defaults to x_max.

  • name (str) – Label for the model.

Return type:

Grid2D

classmethod from_1d_layers(model, *, nx=40, x_max=10000.0, n_pad=8, pad_factor=1.3, n_stations=10, station_x_max=None, dz_min=10.0, name='')[source]

Build a 2-D grid by extending a 1-D LayeredModel.

Every column has the same 1-D layer stack. This is useful as a background model for synthetic anomaly tests and inversion starting models.

Parameters:

  • model (LayeredModel) – Source 1-D model. model.thickness defines the cell thicknesses; the deepest cell becomes the halfspace.

  • nx (int) – Number of core cells in x.

  • x_max (float) – Horizontal extent of the core region [m].

  • n_pad (int) – Padding cells on each side.

  • pad_factor (float) – Padding growth factor.

  • n_stations (int) – Number of equally spaced surface stations.

  • station_x_max (float or None) – Right-most station x [m].

  • dz_min (float) – Minimum cell height [m]; avoids excessively thin cells at depth.

  • name (str) – Model label.

Return type:

Grid2D

classmethod with_anomaly(bg_rho=100.0, anomaly_rho=5.0, anomaly_bounds=(2000.0, 4000.0, 200.0, 800.0), *, nx=50, nz=35, x_max=8000.0, z_max=4000.0, n_pad=8, pad_factor=1.3, n_stations=12, station_x_max=None, name='')[source]

Create a background halfspace with one rectangular anomaly.

Parameters:

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

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

  • anomaly_bounds ((x_lo, x_hi, z_lo, z_hi)) – Anomaly extent in core coordinates [m]: horizontal from x_lo to x_hi, depth from z_lo to z_hi.

  • nx (int) – Core cell counts.

  • nz (int) – Core cell counts.

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

  • z_max (float) – Core region extents [m].

  • n_pad (int) – Padding cells per side.

  • pad_factor (float) – Padding growth factor.

  • n_stations (int) – Number of equally spaced surface stations.

  • station_x_max (float or None) – Right-most station x [m].

  • name (str) – Model label.

Return type:

Grid2D

Examples

Resistive body at mid-depth:

>>> g = Grid2D.with_anomaly(
...     bg_rho=50.0, anomaly_rho=2000.0,
...     anomaly_bounds=(2000.0, 4000.0, 500.0, 1500.0),
... )
classmethod layered_with_fault(model, *, fault_x_m, apparent_dip_deg, throw_m, downthrown_side='right', nx=60, nz=50, x_max=10000.0, z_max=5000.0, n_pad=10, pad_factor=1.35, n_stations=10, station_x_max=None, name='')[source]

Extend a 1-D LayeredModel across x on a uniform fine grid (as halfspace() does, not from_1d_layers()’s coarse one-cell-per-layer discretisation – accurate forward modelling near a thin near-surface layer needs cells much thinner than the layer itself), then offset it across one dipping fault plane.

The fault plane migrates laterally with depth (z / tan(apparent_dip_deg)), so a column’s resistivity comes from evaluating the same 1-D layer stack at a shifted depth (z - throw_m on the downthrown block) rather than a fixed per-column swap – the whole stack drops by throw_m, it is not replaced by a different one.

Parameters:

  • model (LayeredModel) – Background 1-D layer stack.

  • fault_x_m (float) – Along-profile position where the fault reaches the surface, in the same core coordinates as station_x_max (0 at the left edge of the core region, not the padded grid).

  • apparent_dip_deg (float) – Dip of the fault as it appears in this 2-D section, degrees below horizontal, (0, 90]. This is not necessarily the fault’s true 3-D dip unless the profile runs perpendicular to strike – see pycsamt.geology.FaultTrace for the general apparent-dip relationship.

  • throw_m (float) – Vertical displacement of the layer stack on the downthrown side, metres (magnitude only; direction is carried by downthrown_side).

  • downthrown_side ({'right', 'left'}) – Which side of fault_x_m – toward increasing or decreasing x – has the layer stack shifted deeper.

  • nx (int) – Core cell counts, horizontal and vertical.

  • nz (int) – Core cell counts, horizontal and vertical.

  • x_max (float) – Core region extent, metres.

  • z_max (float) – Core region extent, metres.

  • n_pad (int) – Padding cells on each side (left, right, bottom).

  • pad_factor (float) – Padding growth factor.

  • n_stations (int) – Number of equally spaced surface stations.

  • station_x_max (float or None) – Right-most station x [m]. Defaults to x_max.

  • name (str) – Model label.

Return type:

Grid2D

Examples

>>> from pycsamt.forward.synthetic import LayeredModel
>>> from pycsamt.forward.grid2d import Grid2D
>>> m = LayeredModel(resistivity=[80.0, 300.0, 3000.0], thickness=[60.0, 260.0])
>>> g = Grid2D.layered_with_fault(
...     m, fault_x_m=1200.0, apparent_dip_deg=65.0,
...     throw_m=700.0, downthrown_side="right",
...     nx=60, nz=50, x_max=2400.0, z_max=1600.0, n_stations=25,
... )
>>> g.resistivity.shape
(60, 80)
classmethod random(*, nx=40, nz=30, x_max=8000.0, z_max=4000.0, n_pad=8, pad_factor=1.3, n_stations=10, rho_min=1.0, rho_max=10000.0, n_layers=4, lateral_variation=True, seed=None, name='random')[source]

Generate a random layered 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 perturbed by ±30 % in log₁₀ space.

Parameters:

  • nx (int) – Core cell counts.

  • nz (int) – Core cell counts.

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

  • z_max (float) – Core extents [m].

  • n_pad (int) – Padding cells per side.

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

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

  • n_layers (int) – Number of horizontal layers.

  • lateral_variation (bool) – Add column-wise log-normal perturbation.

  • seed (int or Generator or None)

  • name (str)

  • pad_factor (float)

  • n_stations (int)

Return type:

Grid2D

pycsamt.forward.grid2d.make_padding(core_spacing, n_pad, factor=1.3)[source]

Return cell widths for one padding strip.

The first padding cell has width core_spacing * factor and each subsequent cell is factor times wider.

Parameters:

  • core_spacing (float) – Width/depth of the last core cell adjacent to the padding strip [m].

  • n_pad (int) – Number of padding cells.

  • factor (float) – Growth factor between successive padding cells. Values in [1.2, 1.5] give good results for MT; 1.3 is a safe default.

Returns:

Cell widths ordered from the core edge outward.

Return type:

numpy.ndarray, shape (n_pad,)