pycsamt.ai.geology.topography#

Physics-facing topographic surfaces on geological model grids.

Survey extraction remains in pycsamt.topo. This module converts those station samples, or explicit projected samples, into immutable elevation rasters and derives air/earth masks for geological and Maxwell meshes.

Functions

`interpolate_topography`(grid, ...[, ...])	Interpolate projected elevation samples onto a geological grid.
`topography_from_sites`(sites, grid, *[, ...])	Extract station elevations with `pycsamt.topo` and rasterize them.

Classes

TopographicSurface(grid, elevation_m, ...[, ...])

Immutable terrain elevation raster aligned to a geological grid.

class pycsamt.ai.geology.topography.TopographicSurface(grid, elevation_m, reference_elevation_m, vertical_datum='metres above sea level', source='array', interpolation_method='linear', sample_coordinates_m=None, sample_elevation_m=None, station_names=())[source]

Bases: object

Immutable terrain elevation raster aligned to a geological grid.

Parameters:

grid (GeologyGrid) – Grid whose horizontal cell centres are sampled by the surface.
elevation_m (ndarray) – Elevation above vertical_datum. Shape is (nx,) in 2-D or (ny, nx) in 3-D.
reference_elevation_m (float) – Elevation corresponding to geological depth zero. Using the maximum terrain elevation keeps all surface depths non-negative.
vertical_datum (str, default="metres above sea level") – Human-readable vertical datum and unit description.
source (str, default="array") – Provenance label such as "sites", "dem", or "array".
interpolation_method ({"linear", "nearest", "cubic"}, default="linear") – Method used to rasterize samples.
sample_coordinates_m (ndarray or None, optional) – Projected sample x positions in 2-D or x/y positions in 3-D.
sample_elevation_m (ndarray or None, optional) – Elevations corresponding to sample coordinates.
station_names (sequence of str, optional) – Station identifiers corresponding to samples.

Examples

>>> grid = GeologyGrid.regular_2d(nx=4, nz=3, dx_m=100, dz_m=50)
>>> surface = TopographicSurface(grid, [100, 110, 105, 95], 110)
>>> surface.relief_m
15.0
>>> surface.earth_mask().shape
(3, 4)

grid: GeologyGrid

elevation_m: ndarray

reference_elevation_m: float

vertical_datum: str = 'metres above sea level'

source: str = 'array'

interpolation_method: str = 'linear'

sample_coordinates_m: ndarray | None = None

sample_elevation_m: ndarray | None = None

station_names: tuple[str, ...] = ()

property relief_m: float[source]

Return maximum minus minimum terrain elevation.

Returns:: Relief in metres.
Return type:: float

Examples

>>> grid = GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1)
>>> TopographicSurface(grid, [10, 13], 13).relief_m
3.0

property surface_depth_m: ndarray[source]

Return terrain depth below the configured reference elevation.

Returns:: Horizontal surface shaped like elevation_m. Positive values lie below reference depth zero; negative values lie above it.
Return type:: ndarray

Examples

>>> grid = GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1)
>>> TopographicSurface(grid, [10, 8], 10).surface_depth_m.tolist()
[0.0, 2.0]

property surface_hash: str[source]

Return a platform-stable terrain and provenance digest.

Returns:: SHA-256 digest.
Return type:: str

Examples

>>> grid = GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1)
>>> len(TopographicSurface(grid, [10, 9], 10).surface_hash)
64

local_depth_m()[source]

Return each cell centre’s signed depth below local terrain.

Returns:: Array shaped like grid. Non-negative values are in the earth; negative values are above terrain and belong to the air region.
Return type:: ndarray

Examples

>>> grid = GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1)
>>> surface = TopographicSurface(grid, [10, 9], 10)
>>> surface.local_depth_m().shape
(2, 2)

earth_mask()[source]

Return cells on or below the local terrain surface.

Returns:: Read-only physics-facing earth mask shaped like grid.
Return type:: ndarray of bool

Examples

>>> grid = GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1)
>>> TopographicSurface(grid, [10, 9], 10).earth_mask().dtype == bool
True

air_mask()[source]

Return cells strictly above the local terrain surface.

Returns:: Logical complement of earth_mask().
Return type:: ndarray of bool

Examples

>>> grid = GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1)
>>> surface = TopographicSurface(grid, [10, 9], 10)
>>> np.array_equal(surface.air_mask(), ~surface.earth_mask())
True

slope_degrees()[source]

Return terrain slope magnitude at horizontal cell centres.

Returns:: Slope angle in degrees, shaped like elevation_m.
Return type:: ndarray

Examples

>>> grid = GeologyGrid.regular_2d(nx=3, nz=2, dx_m=10, dz_m=1)
>>> slope = TopographicSurface(grid, [0, 10, 20], 20).slope_degrees()
>>> np.allclose(slope, 45)
True

summary()[source]

Return compact JSON-compatible terrain diagnostics.

Returns:: Elevation range, relief, slopes, air fraction, source, and datum.
Return type:: dict

Examples

>>> grid = GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1)
>>> TopographicSurface(grid, [10, 9], 10).summary()["relief_m"]
1.0

provenance()[source]

Return complete JSON-compatible surface provenance.

Returns:: Grid, datum, source, interpolation, samples, and station names.
Return type:: dict

Examples

>>> grid = GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1)
>>> TopographicSurface(grid, [10, 9], 10).provenance()["source"]
'array'

to_npz(path)[source]

Persist terrain and provenance in a pickle-free NPZ archive.

Parameters:: path (str or pathlib.Path) – Destination archive.
Returns:: Requested destination.
Return type:: pathlib.Path

Examples

>>> from tempfile import TemporaryDirectory
>>> surface = TopographicSurface(
...     GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1), [10, 9], 10
... )
>>> with TemporaryDirectory() as directory:
...     restored = TopographicSurface.from_npz(
...         surface.to_npz(Path(directory) / "topo.npz")
...     )
>>> restored.surface_hash == surface.surface_hash
True

classmethod from_npz(path)[source]

Load and validate a topographic archive without enabling pickle.

Parameters:: path (str or pathlib.Path) – Archive written by to_npz().
Returns:: Immutable restored surface.
Return type:: TopographicSurface

Examples

>>> from tempfile import TemporaryDirectory
>>> surface = TopographicSurface(
...     GeologyGrid.regular_2d(nx=2, nz=2, dx_m=1, dz_m=1), [10, 9], 10
... )
>>> with TemporaryDirectory() as directory:
...     restored = TopographicSurface.from_npz(
...         surface.to_npz(Path(directory) / "t.npz")
...     )
>>> np.array_equal(restored.elevation_m, surface.elevation_m)
True

pycsamt.ai.geology.topography.interpolate_topography(grid, sample_coordinates_m, sample_elevation_m, *, interpolation_method='linear', reference_elevation_m=None, vertical_datum='metres above sea level', source='array', station_names=())[source]

Interpolate projected elevation samples onto a geological grid.

Parameters:

grid (GeologyGrid) – Target 2-D or 3-D grid.
sample_coordinates_m (ndarray) – Shape (n,) or (n, 1) x positions in 2-D; shape (n, 2) x/y projected coordinates in 3-D.
sample_elevation_m (ndarray, shape (n,)) – Finite elevations in metres above the declared datum.
interpolation_method ({"linear", "nearest", "cubic"}, default="linear") – Interpolation method. Values outside the convex hull are filled from nearest samples rather than extrapolated polynomials.
reference_elevation_m (float or None, optional) – Geological depth-zero elevation. Default is maximum raster elevation.
vertical_datum (str, optional) – Provenance labels.
source (str, optional) – Provenance labels.
station_names (sequence of str, optional) – Unique labels for every sample.

Returns:

Immutable raster and original samples.

Return type:

TopographicSurface

Examples

>>> grid = GeologyGrid.regular_2d(nx=4, nz=3, dx_m=100, dz_m=50)
>>> surface = interpolate_topography(grid, [0, 400], [100, 120])
>>> surface.elevation_m.shape
(4,)

pycsamt.ai.geology.topography.topography_from_sites(sites, grid, *, station_names=None, coordinates_m=None, interpolation_method='linear', reference_elevation_m=None, vertical_datum='metres above sea level', profile_origin_m=None, allow_all_zero=False)[source]

Extract station elevations with pycsamt.topo and rasterize them.

Parameters:

sites (Sites or EDI-like collection) – Any container accepted by pycsamt.topo.extract_elevation().
grid (GeologyGrid) – Target geology grid.
station_names (sequence of str or None, optional) – Requested station order/subset. Matching is case-insensitive.
coordinates_m (ndarray or None, optional) – Required projected x/y coordinates for a 3-D grid. For 2-D, when supplied as (n, 2), cumulative projected chainage replaces the latitude/longitude-derived chainage from pycsamt.topo.
interpolation_method (str) – Forwarded to interpolate_topography().
reference_elevation_m (float | None) – Forwarded to interpolate_topography().
vertical_datum (str) – Forwarded to interpolate_topography().
profile_origin_m (float or None, optional) – X coordinate assigned to zero chainage in 2-D. Default is the grid’s minimum outer x edge.
allow_all_zero (bool, default=False) – Permit an all-zero elevation collection. False rejects the ambiguity between genuine sea-level terrain and missing EDI elevations.

Returns:

Physics-facing raster with original station samples retained.

Return type:

TopographicSurface

Raises:

ValueError – If station alignment, elevations, chainage, or required 3-D projected coordinates are invalid.

Examples

>>> from types import SimpleNamespace
>>> stations = [
...     SimpleNamespace(
...         Head=SimpleNamespace(dataid="S00", elev=100, lat=5, lon=-3)
...     ),
...     SimpleNamespace(
...         Head=SimpleNamespace(dataid="S01", elev=110, lat=5, lon=-2.999)
...     ),
... ]
>>> grid = GeologyGrid.regular_2d(nx=3, nz=2, dx_m=50, dz_m=20)
>>> surface = topography_from_sites(stations, grid)
>>> surface.source, surface.station_names
('sites', ('S00', 'S01'))