2.29.3.2. pycsamt.topo.drape#

Terrain-following coordinate transform for 2-D section plots.

A standard 2-D resistivity section stores:

  • x — along-profile distance (km)

  • z — depth below a flat datum at z = 0 (km, positive downward)

When real surface topography is available, the datum is no longer flat. This module transforms the flat depth-grid into a terrain-following coordinate frame in which z = 0 at every x position corresponds to the local surface elevation:

z_real(x, z) = elev(x) - z [both in km]

The result can be fed directly to matplotlib.pyplot.pcolormesh() using a 2-D Z argument (supported since Matplotlib 3.3) so that the mesh cells drape over the terrain.

Typical usage:

from pycsamt.topo.drape import interp_elev, drape_section

elev_at_x = interp_elev(chain_km, elev_km, x_centres_km)
x_nodes, z_draped, data = drape_section(
    x_nodes, z_nodes, rho_2d, elev_at_x
)
ax.pcolormesh(x_nodes, z_draped, data, ...)

Functions

drape_section(x_nodes, z_nodes, data, ...[, ...])

Transform a flat depth section into terrain-following coordinates.

interp_elev(chainage_km, elev_km, x_query_km)

Interpolate station elevations to arbitrary profile positions.

mask_above_topo(x_nodes, z_nodes, data, ...)

Set data cells that lie above the terrain surface to NaN.

station_surface_z(chainage_km, elev_km, ...)

Return the terrain-draped z-coordinate for station marker positions.
pycsamt.topo.drape.interp_elev(chainage_km, elev_km, x_query_km, method='linear')[source]

Interpolate station elevations to arbitrary profile positions.

Clamps extrapolated values to the boundary station elevations so the terrain surface never shoots up unexpectedly at the section edges.

Parameters:

  • chainage_km (array_like, shape (n_stations,)) – Along-profile distances of the stations (km).

  • elev_km (array_like, shape (n_stations,)) – Terrain elevation at each station (km a.s.l.).

  • x_query_km (array_like, shape (m,)) – Positions at which to evaluate the interpolated elevation (km).

  • method ({"linear", "cubic", "nearest"}) – Interpolation method. "cubic" requires scipy.

Returns:

Interpolated elevation in km a.s.l.

Return type:

numpy.ndarray, shape (m,)

pycsamt.topo.drape.drape_section(x_nodes, z_nodes, data, elev_at_centres, exaggeration=1.0, clip_above_surface=False)[source]

Transform a flat depth section into terrain-following coordinates.

Builds a 2-D z_draped array (shape (nz+1, nx+1)) where each column is shifted so that z_nodes[0] (the surface) aligns with the local terrain elevation. The result can be passed directly to pcolormesh() as the Y argument.

Parameters:

  • x_nodes (array_like, shape (nx+1,)) – Horizontal pcolormesh node positions (km).

  • z_nodes (array_like, shape (nz+1,)) – Depth node positions (km, positive downward from flat datum). z_nodes[0] should be 0 (surface) or the shallowest depth.

  • data (array_like, shape (nz, nx)) – 2-D data values (e.g. log10(rho)).

  • elev_at_centres (array_like, shape (nx,)) – Terrain elevation at each cell-centre x position (km a.s.l.).

  • exaggeration (float) – Vertical exaggeration applied to both the elevation offset and the depth axis. Values > 1 amplify relief for display purposes.

  • clip_above_surface (bool) – If True, set cells that are above the terrain surface (unreachable subsurface) to NaN.

Returns:

  • x_nodes (numpy.ndarray, shape (nx+1,)) – Unchanged horizontal node positions.

  • z_draped (numpy.ndarray, shape (nz+1, nx+1)) – 2-D elevation node array. Column j equals elev_at_nodes[j] - z_nodes * exaggeration.

  • data_out (numpy.ndarray, shape (nz, nx)) – Input data, optionally NaN-masked above terrain.

Return type:

tuple[ndarray, ndarray, ndarray]

Notes

The terrain-following coordinate for node (k, j) is:

z_draped[k, j] = elev_nodes[j] - z_nodes[k] * exaggeration

where elev_nodes is the terrain interpolated to the x node positions (nx+1 values) from the cell-centre values.

pycsamt.topo.drape.mask_above_topo(x_nodes, z_nodes, data, elev_at_centres, exaggeration=1.0)[source]

Set data cells that lie above the terrain surface to NaN.

A cell at column j and depth-row k has absolute elevation:

cell_elev = elev_at_centres[j] - z_centre[k] * exaggeration

where z_centre[k] = (z_nodes[k] + z_nodes[k+1]) / 2.

A cell is above the surface when cell_elev > elev_at_centres[j], which simplifies to z_centre[k] < 0. For standard meshes (all z ≥ 0) this never occurs; the masking is only relevant for meshes whose z-origin is below the deepest station (uncommon).

The more practically useful masking is for varying terrain: station A sits at 800 m, station B sits at 200 m. A cell at depth 500 m below A has absolute elevation 300 m — it is above the surface at B. This function masks such cells so they do not appear in the plot.

Parameters:

  • x_nodes ((nx+1,) node positions)

  • z_nodes ((nz+1,) depth nodes (km, positive down))

  • data ((nz, nx) data array)

  • elev_at_centres ((nx,) elevation at cell centres (km))

  • exaggeration (float)

Returns:

Data with cells above the terrain set to NaN.

Return type:

numpy.ndarray, shape (nz, nx)

pycsamt.topo.drape.station_surface_z(chainage_km, elev_km, station_x_km, exaggeration=1.0)[source]

Return the terrain-draped z-coordinate for station marker positions.

In the draped coordinate frame the station sits at the surface elevation, not at z = 0.

Parameters:

  • chainage_km ((n,) station chainage values)

  • elev_km ((n,) elevation at each station (km))

  • station_x_km ((m,) x positions where markers should be drawn)

  • exaggeration (float)

Returns:

z-coordinates (km) at which to place station markers.

Return type:

numpy.ndarray, shape (m,)