3-D maps with topography#

Real surveys are not flat. pycsamt.map can drape the 3-D views over surface elevation, so cross-sections and volumes hang beneath the true terrain rather than a flat datum — essential for shallow targets where a few tens of metres of relief matter. WILLY_DATA carries per-station elevations (37–224 m), so topography here is real, not synthetic.

The scenes are interactive: orbit to see the terrain surface and the survey draped beneath it.

Load the survey#

import os

import numpy as np

from pycsamt.map import MapView


DATA = os.path.join(
    os.environ.get("PYCSAMT_DOCS_REPO_ROOT", "."), "data", "AMT", "WILLY_DATA"
)
mv = MapView.from_folder(DATA, recursive=True)
elev = mv.table()["Elevation"].to_numpy(dtype=float)
print(f"elevation range: {np.nanmin(elev):.0f}{np.nanmax(elev):.0f} m")
elevation range: 37–177 m

Baseline: a flat datum#

Turning topography off pins every panel to a flat z=0 datum — the reference to compare against. This is the pure cross-section geometry, ignoring surface relief.

fig = mv.map3d(mode="fence", topography=False, show_terrain=False)
fig.update_layout(height=640, scene_aspectmode="cube")
fig


Draped over real terrain#

Topography is on by default: each panel is shifted to sit beneath its station’s true elevation, and show_terrain adds the interpolated ground surface. Compare with the flat fence above — the panels now follow the relief.

fig = mv.map3d(mode="fence")  # topography=True is the default
fig.update_layout(height=660, scene_aspectmode="cube")
fig


Volume beneath the terrain#

The block view honours topography too — the resistivity body hangs under the real ground surface, so depths read from the surface down. A lighter terrain lets the volume show through.

fig = mv.map3d(
    mode="block", terrain_opacity=0.5, show_stations=True, station_size=3
)
fig.update_layout(height=660, scene_aspectmode="cube")
fig


Supplying your own elevations#

When EDI headers lack elevation (or you have a better DEM), attach one with MapView.with_elevations — a {station_id: elevation_m} mapping. MapView.fetch_elevations() can also pull a DEM online. Here we exaggerate the relief to show the effect.

boosted = {
    s: float(200 + 3 * (e - np.nanmean(elev)))
    for s, e in zip(mv.stations, elev)
}
fig = mv.with_elevations(boosted).map3d(mode="fence")
fig.update_layout(height=660, scene_aspectmode="cube")
fig


Next. The advanced compositions example tunes colour, slicing, iso-values, and the camera for presentation-quality figures.

Total running time of the script: (0 minutes 0.478 seconds)

Gallery generated by Sphinx-Gallery