from __future__ import annotations
from typing import Any
import matplotlib.colors as mcolors
import matplotlib.pyplot as plt
from matplotlib.transforms import Affine2D
import numpy as np
import pandas as pd
from matplotlib.cm import ScalarMappable
from matplotlib.colors import BoundaryNorm, ListedColormap, Normalize
from matplotlib.patches import Ellipse
from mpl_toolkits.axes_grid1 import make_axes_locatable
from ..api._rose_style import (
_UNSET,
RoseStyle,
resolve_rose_style,
)
from ..api.labels import LOG10_PERIOD_LABEL, PERIOD_LABEL
from ..api.station import PYCSAMT_STATION_RENDERING
from ..api.style import PYCSAMT_STYLE
from ..site import compute as _compute
from ..site import edit as _edit
from ..site.base import Sites
from ..z import utils as zutils
from ._core import (
_apply_each,
_axes_list,
_get_t_block,
_get_z_block,
_iter_items,
_name,
ensure_sites,
hide_polar_radius_labels,
)
_BACKWARD_SINCE = "2.0.0"
_BACKWARD_REMOVE = "2.17.0"
#
[docs]
def rotate_to_strike(
sites,
*,
method: str = "swift",
inplace: bool = False,
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
) -> Sites:
S = ensure_sites(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
def _one(Si: Sites) -> Sites:
# _compute.strike_estimate only recognises objects exposing a
# ``.Z`` section (its own is_edi_file() duck-type check); a
# Site wrapper only exposes ``.z`` and is silently skipped,
# which made this always fall through to ang=0.0 regardless of
# data. Passing the raw EDI item makes it return a real float
# (per its own documented single-site contract) instead of an
# empty per-station DataFrame with no ``.angle`` column.
ed = next(_iter_items(Si))
edi = getattr(ed, "edi", None)
edi = edi if edi is not None else ed
ang = 0.0
try:
r = _compute.strike_estimate(edi, method=method)
if isinstance(r, (int, float)):
ang = float(r)
elif isinstance(r, dict):
ang = float(r.get("angle", 0.0))
else:
ang = float(getattr(r, "angle", 0.0))
except Exception:
ang = 0.0
if not np.isfinite(ang):
# No strike (e.g. every row masked upstream): rotating by NaN
# used to turn the whole tensor -- and its errors -- into NaN.
return Si
# _edit.rotate looks for a ``.Z`` section (the raw EDI layout);
# a Site wrapper only exposes ``.z`` directly, so calling it on
# `Si` (a Sites collection) is a silent no-op for every real
# station. Site.z reads through to edi.Z.z, so rotating the
# underlying EDI in place also updates the Site wrapper.
# _apply_each already handles inplace-vs-copy at the collection
# level, so this always mutates in place here regardless of
# the outer *inplace* flag.
_edit.rotate(edi, ang, inplace=True)
return Si
return _apply_each(S, _one, inplace=inplace, verbose=verbose)
# --- 2) rotate by angle or map ------------------------------------------- #
[docs]
def rotate(
sites,
angle: float,
*,
inplace: bool = False,
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
) -> Sites:
S = ensure_sites(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
# _edit.rotate expects a single EDI-like object with a ``.Z``
# section; passing the whole Sites collection silently rotates
# nothing. rotate_all() is the broadcast variant that correctly
# unwraps each item to its raw EDI first.
return _edit.rotate_all(S, angle, inplace=inplace)
[docs]
def rotate_by_map(
sites,
angle_by_station: dict,
*,
inplace: bool = False,
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
) -> Sites:
S = ensure_sites(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
def _one(Si: Sites) -> Sites:
ed = next(_iter_items(Si))
name = getattr(ed, "station", None) or getattr(ed, "name", None)
ang = float(angle_by_station.get(name, 0.0))
# See rotate_to_strike._one above: rotate the raw EDI item,
# not the Sites wrapper, and always in place (_apply_each
# handles the outer inplace-vs-copy semantics).
edi = getattr(ed, "edi", None)
_edit.rotate(edi if edi is not None else ed, ang, inplace=True)
return Si
return _apply_each(S, _one, inplace=inplace, verbose=verbose)
# --- 3) enforce off-diagonal antisymmetry -------------------------------- #
[docs]
def antisymmetrize(
sites,
*,
how: str = "rms",
inplace: bool = False,
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
) -> Sites:
S = ensure_sites(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
def _one(Si: Sites) -> Sites:
ed = next(_iter_items(Si))
Z, z, _ = _get_z_block(ed)
if Z is None or z is None:
return Si
ze = getattr(Z, "z_err", None)
z2, ze2 = zutils.enforce_offdiag_antisymmetry(z, ze)
try:
Z.z = z2
if ze is not None and ze2 is not None:
Z.z_err = ze2
except Exception:
pass
return Si
return _apply_each(S, _one, inplace=inplace, verbose=verbose)
# --- 4) invert impedance -------------------------------------------------- #
[docs]
def invert(
sites,
*,
inplace: bool = False,
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
) -> Sites:
S = ensure_sites(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
def _one(Si: Sites) -> Sites:
ed = next(_iter_items(Si))
Z, z, _ = _get_z_block(ed)
if Z is None or z is None:
return Si
ze = getattr(Z, "z_err", None)
z2, ze2 = zutils.invert_z(z, ze)
try:
Z.z = z2
if ze is not None and ze2 is not None:
Z.z_err = ze2
except Exception:
pass
return Si
return _apply_each(S, _one, inplace=inplace, verbose=verbose)
# --- 5) sensor-orientation correction ------------------------------------ #
[docs]
def orient_from_sensors(
sites,
ex: float,
ey: float,
bx: float,
by: float,
*,
degrees: bool = True,
inplace: bool = False,
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
) -> Sites:
S = ensure_sites(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
# zutils.correct_for_sensor_orientation always expects degrees and
# has no ``degrees=``/``z_err=`` keywords (it takes ``z_prime_err``);
# calling it with those names raised TypeError unconditionally, so
# this function could never succeed. Convert here when the caller
# passed radians instead.
if not degrees:
ex_, ey_, bx_, by_ = (float(np.degrees(v)) for v in (ex, ey, bx, by))
else:
ex_, ey_, bx_, by_ = float(ex), float(ey), float(bx), float(by)
def _one(Si: Sites) -> Sites:
ed = next(_iter_items(Si))
Z, z, _ = _get_z_block(ed)
if Z is None or z is None:
return Si
ze = getattr(Z, "z_err", None)
z2, ze2 = zutils.correct_for_sensor_orientation(
z,
ex=ex_,
ey=ey_,
bx=bx_,
by=by_,
z_prime_err=ze,
)
try:
Z.z = z2
if ze is not None and ze2 is not None:
Z.z_err = ze2
except Exception:
pass
return Si
return _apply_each(S, _one, inplace=inplace, verbose=verbose)
# --- 6) sigma-clip outliers in Z ----------------------------------------- #
[docs]
def sigma_clip_z(
sites,
*,
sigma: float = 3.0,
inplace: bool = False,
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
) -> Sites:
S = ensure_sites(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
def _one(Si: Sites) -> Sites:
ed = next(_iter_items(Si))
Z, z, _ = _get_z_block(ed)
if Z is None or z is None:
return Si
mask = zutils.sigma_clip_mask(z, nsigma=sigma)
z2 = z.copy()
z2[~mask] = np.nan
try:
Z.z = z2
except Exception:
pass
return Si
return _apply_each(S, _one, inplace=inplace, verbose=verbose)
# --- 7) balance off-diagonal magnitudes ---------------------------------- #
[docs]
def balance_offdiag(
sites,
*,
mode: str = "avgabs",
inplace: bool = False,
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
) -> Sites:
S = ensure_sites(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
def _one(Si: Sites) -> Sites:
ed = next(_iter_items(Si))
Z, z, _ = _get_z_block(ed)
if Z is None or z is None:
return Si
zx = z[:, 0, 1]
zy = z[:, 1, 0]
if mode == "avgabs":
m = 0.5 * (np.abs(zx) + np.abs(zy)) + 1e-12
sx = m / (np.abs(zx) + 1e-12)
sy = m / (np.abs(zy) + 1e-12)
z2 = z.copy()
z2[:, 0, 1] = zx * sx
z2[:, 1, 0] = zy * sy
else:
z2 = z
try:
Z.z = z2
except Exception:
pass
return Si
return _apply_each(S, _one, inplace=inplace, verbose=verbose)
def _is_z_tensor_shaped(z: np.ndarray) -> bool:
"""Return whether *z* looks like a proper (n, 2, 2) impedance tensor."""
return (z.ndim == 3 and z.shape[1:3] == (2, 2)) or (
z.ndim == 2 and z.shape == (2, 2)
)
def _get_z(ed: Any) -> tuple[np.ndarray | None, np.ndarray | None]:
# ``getattr(ed, "Z", None) or getattr(ed, "z", None)`` silently drops
# valid data for Site-wrapped objects that only expose a lowercase
# ``.z`` (raw array, no ``.freq`` sibling) with the real ``Z`` wrapper
# living at ``ed.edi.Z`` -- the same pattern documented in
# rotate_to_strike(). That path returned (None, None) for such sites
# even though `_get_z_block` (this module's sibling in `_core.py`)
# resolves them fine, making build_phase_tensor_table silently empty
# for real EDI batches that construct Site wrappers this way.
#
# The `ed.edi.Z` fallback below reaches past `Site.z`'s own shape
# validation straight to the raw EDI wrapper, so a malformed
# ``Z.z`` (not (n,2,2)) that `Site.z` would normally reject as None
# must be re-validated here explicitly with `_is_z_tensor_shaped`,
# or callers like `_finite_z_period_rows` crash on it instead of
# gracefully treating the station as "no impedance data".
Zobj = getattr(ed, "Z", None)
if Zobj is None:
edi = getattr(ed, "edi", None)
if edi is not None:
Zobj = getattr(edi, "Z", None)
if Zobj is None:
Zobj = getattr(ed, "z", None)
if Zobj is None:
return None, None
z = getattr(Zobj, "z", None)
fr = getattr(Zobj, "freq", None)
if (
isinstance(z, np.ndarray)
and isinstance(fr, np.ndarray)
and _is_z_tensor_shaped(z)
):
return z, fr
if isinstance(Zobj, np.ndarray) and _is_z_tensor_shaped(Zobj):
fr2 = getattr(ed, "freq", None)
if isinstance(fr2, np.ndarray):
return Zobj, fr2
return None, None
def _period(fr: np.ndarray) -> np.ndarray:
with np.errstate(divide="ignore", invalid="ignore"):
p = 1.0 / fr
return p
def _angles_deg(
a: np.ndarray, b: np.ndarray, c: np.ndarray, d: np.ndarray
) -> tuple[np.ndarray, np.ndarray]:
"""Phase-tensor orientation angle alpha and skew angle beta (degrees).
Caldwell, Bibby & Brown (2004) convention, Phi = [[a, b], [c, d]]:
alpha = 0.5 * atan2(b + c, a - d) -- rotation-variant orientation
beta = 0.5 * atan2(b - c, a + d) -- rotation-invariant skew
alpha shifts with coordinate rotation (like the strike angle); beta
does not, which is what makes it usable as a distortion/dimensionality
diagnostic. A previous version of this function swapped the two
formulas (and additionally sign-flipped the skew one), so the
dataframe's "beta" column tracked coordinate rotation while "alpha"
stayed invariant -- exactly backwards. Confirmed by rotating a
synthetic 2-D tensor: only the corrected "beta" below stays fixed.
"""
alpha = 0.5 * np.degrees(np.arctan2(b + c, a - d))
beta = 0.5 * np.degrees(np.arctan2(b - c, a + d))
return alpha, beta
def _phi_from_z(z: np.ndarray) -> np.ndarray:
# z: (n, 2, 2) complex
X = z.real
Y = z.imag
n = z.shape[0]
phi = np.empty((n, 2, 2), dtype=float)
for i in range(n):
Xi = X[i]
Yi = Y[i]
try:
invX = np.linalg.inv(Xi)
except np.linalg.LinAlgError:
invX = np.linalg.pinv(Xi)
phi[i] = invX @ Yi
return phi
def _finite_z_period_rows(
z: np.ndarray,
fr: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
"""Return rows with finite impedance tensors and positive frequency."""
z = np.asarray(z, dtype=complex)
fr = np.asarray(fr, dtype=float)
n = min(z.shape[0], fr.size)
if n == 0:
return z[:0], fr[:0]
z = z[:n]
fr = fr[:n]
finite_z = np.isfinite(z.real) & np.isfinite(z.imag)
keep = np.isfinite(fr) & (fr > 0.0) & np.all(finite_z, axis=(1, 2))
return z[keep], fr[keep]
def _svd_axes(phi: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
# returns s (n,2) and angle of major axis (deg), shape (n,)
n = phi.shape[0]
s = np.empty((n, 2), dtype=float)
ang = np.empty(n, dtype=float)
for i in range(n):
U, S, VT = np.linalg.svd(phi[i])
s[i] = S
vx = U[0, 0]
vy = U[1, 0]
ang[i] = np.degrees(np.arctan2(vy, vx))
return s, ang
def _ellipticity(s: np.ndarray) -> np.ndarray:
num = s[:, 0] - s[:, 1]
den = s[:, 0] + s[:, 1] + 1e-12
return num / den
def _df_phase_tensor_single(
st: str,
fr: np.ndarray,
phi: np.ndarray,
s: np.ndarray,
theta: np.ndarray,
alpha: np.ndarray,
beta: np.ndarray,
) -> pd.DataFrame:
p = _period(fr)
e = _ellipticity(s)
rec = dict(
station=[st] * len(fr),
freq=fr,
period=p,
s1=s[:, 0],
s2=s[:, 1],
theta=theta,
alpha=alpha,
beta=beta,
skew=beta,
ellipt=e,
)
return pd.DataFrame(rec)
# ------------------------ public computation ----------------------------- #
[docs]
def build_phase_tensor_table(
sites: Any,
*,
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
) -> pd.DataFrame:
S = ensure_sites(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
dfs: list[pd.DataFrame] = []
for i, ed in enumerate(_iter_items(S)):
nm = _name(ed, i)
z, fr = _get_z(ed)
if z is None or fr is None:
continue
z, fr = _finite_z_period_rows(z, fr)
if z.size == 0 or fr.size == 0:
continue
phi = _phi_from_z(z)
finite_phi = np.all(np.isfinite(phi), axis=(1, 2))
if not np.any(finite_phi):
continue
phi = phi[finite_phi]
fr = fr[finite_phi]
a = phi[:, 0, 0]
b = phi[:, 0, 1]
c = phi[:, 1, 0]
d = phi[:, 1, 1]
alpha, beta = _angles_deg(a, b, c, d)
s, theta = _svd_axes(phi)
df = _df_phase_tensor_single(nm, fr, phi, s, theta, alpha, beta)
dfs.append(df)
if not dfs:
return pd.DataFrame(
columns=[
"station",
"freq",
"period",
"s1",
"s2",
"theta",
"alpha",
"beta",
"skew",
"ellipt",
]
)
return pd.concat(dfs, ignore_index=True)
# ─── Plotting helpers ─────────────────────────────────────────────────────────
# Human-readable labels for each c_by specifier
_C_LABEL: dict[str, str] = {
"skew": "Skew β (°)",
"beta": "Skew β (°)",
"alpha": "α (°)",
"theta": "Strike θ (°)",
"ellipt": "Ellipticity",
"s1": "φ_max (s.v.)",
"s2": "φ_min (s.v.)",
"|skew|": "|β| (°)",
"|beta|": "|β| (°)",
"|theta|": "|θ| (°)",
"|alpha|": "|α| (°)",
"phi_mean": "Φ̄ = (φ_max+φ_min)/2",
"phi_max": "φ_max (tan units)",
"phi_min": "φ_min (tan units)",
"phimax_deg": r"$\phi_{max}$ (°)",
"phimin_deg": r"$\phi_{min}$ (°)",
}
# Symmetric c_by quantities (vmin = -vmax enforced by default).
# Absolute-value variants ("|skew|", "|beta|") are deliberately excluded:
# they are >= 0 by construction, so forcing vmin = -vmax would burn half
# the colour range on values that can never occur and desaturate the
# colours actually used for the 0..vmax data (see also
# PhaseTensorEllipseStyle's per-c_by natural_sym flag in api/style.py,
# which already gets this right and is the source of truth when
# symmetric_clim is left at its _UNSET default).
_SYMMETRIC_C: frozenset = frozenset(("skew", "beta", "alpha"))
def _resolve_cvals(df: pd.DataFrame, c_by: str) -> tuple[np.ndarray, str]:
"""Return ``(colour_array, colorbar_label)`` for a *c_by* specifier."""
if c_by.startswith("|") and c_by.endswith("|"):
inner = c_by[1:-1]
col = (df[inner] if inner in df.columns else df["skew"]).to_numpy(
float
)
return np.abs(col), _C_LABEL.get(c_by, f"|{inner}|")
if c_by == "phi_mean":
return (
0.5 * (df["s1"].to_numpy(float) + df["s2"].to_numpy(float)),
_C_LABEL["phi_mean"],
)
if c_by == "phi_max":
return df["s1"].to_numpy(float), _C_LABEL["phi_max"]
if c_by == "phi_min":
return df["s2"].to_numpy(float), _C_LABEL["phi_min"]
if c_by == "phimax_deg":
# s1/s2 are the phase-tensor singular values (tan units); the
# classic MTpy-style 0-90 deg phimin/phimax colouring needs the
# arctan of those, not the raw tan-unit values "phi_max"/"phi_min"
# above return (those saturate a 0-90 colour scale near 0).
return (
np.degrees(np.arctan(df["s1"].to_numpy(float))),
_C_LABEL["phimax_deg"],
)
if c_by == "phimin_deg":
return (
np.degrees(np.arctan(df["s2"].to_numpy(float))),
_C_LABEL["phimin_deg"],
)
col_name = c_by if c_by in df.columns else "skew"
return df[col_name].to_numpy(float), _C_LABEL.get(c_by, c_by)
def _attach_cbar(
ax: plt.Axes,
mappable,
label: str,
*,
size: str = "3%",
pad: float = 0.06,
labelsize: int = 9,
ticksize: int = 8,
**colorbar_kws,
):
"""Attach a right-side colorbar to *ax* via make_axes_locatable."""
div = make_axes_locatable(ax)
cax = div.append_axes("right", size=size, pad=pad)
cb = ax.get_figure().colorbar(mappable, cax=cax, **colorbar_kws)
cb.set_label(label, fontsize=labelsize)
cb.ax.tick_params(labelsize=ticksize)
return cb
# ------------------------------ plotting --------------------------------- #
[docs]
def plot_phase_tensor_psection(
sites: Any,
*,
# ── data selection ────────────────────────────────────────────────────
stations: list[str] | None = None,
period_range: tuple[float, float] | None = None,
# ── y-axis ────────────────────────────────────────────────────────────
axis_y: str = "logperiod",
period_up: bool = False,
frame_pct: tuple[float, float] | None = (1.0, 99.0),
# ── ellipse sizing ────────────────────────────────────────────────────
scale=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.scale
normalise_by=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.normalise_by
s1_ref: float | None = None,
min_aspect=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.min_aspect
# ── colour ────────────────────────────────────────────────────────────
c_by=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.c_by
cmap=_UNSET, # default: auto from c_by via resolve_cmap()
clim: tuple[float, float] | None = None,
clim_pct=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.clim_pct
symmetric_clim=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.resolve_symmetric_clim()
color_mode: str = "continuous",
segment_bounds: tuple[float, float] | None = None,
segment_colors: tuple[str, str, str] = ("#2166ac", "#f7f7f7", "#b2182b"),
# ── aesthetics ────────────────────────────────────────────────────────
edgecolor=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.edgecolor
linewidth=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.linewidth
alpha=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.alpha
ellipse_kws: dict[str, Any] | None = None,
cb_kws: dict[str, Any] | None = None,
colorbar: bool = True,
# ── overlays ─────────────────────────────────────────────────────────
skew_threshold=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.skew_threshold
mark_3d=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.mark_3d
ref_ellipse=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.show_ref
legend_fontsize: float = 8.0,
annotations: bool = True,
# ── labels & layout ───────────────────────────────────────────────────
title: str = "",
xlabel: str = "",
ylabel: str = "",
tick_label_rotation: float = 45.0,
figsize: tuple[float, float] = (10.0, 5.5),
# ── standard emtools ──────────────────────────────────────────────────
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
ax: plt.Axes | None = None,
) -> plt.Axes:
"""
Phase-tensor ellipse pseudo-section (Caldwell et al. 2004 style).
Each cell ``(station × period)`` is drawn as an ellipse whose shape,
orientation, and fill colour encode phase-tensor invariants:
* **major axis** ∝ φ_max (s1 eigenvalue — mean phase level)
* **minor axis** ∝ φ_min (s2 eigenvalue — minimum phase)
* **aspect ratio** = φ_min / φ_max (1 = isotropic / 1-D)
* **rotation** = θ (geoelectric strike, CCW from E)
* **fill colour** controlled by *c_by* (default: skew angle β)
Parameters
----------
sites : any
EDI paths, glob pattern, :class:`~pycsamt.site.base.Sites`, or any
input accepted by :func:`~pycsamt.emtools.ensure_sites`.
stations : list of str or None
Restrict to a subset of station names.
period_range : (T_min, T_max) in seconds or None
Restrict the period range plotted.
axis_y : ``"logperiod"`` | ``"logfreq"``
Y-axis quantity. ``"logperiod"`` is the MT convention.
period_up : bool, default ``False``
When ``False`` (default; standard pseudo-section convention, as
used elsewhere in :mod:`pycsamt.emtools`, e.g.
:func:`~pycsamt.emtools.tensor.plot_dimensionality_psection`),
short period (high frequency, shallow sensitivity) is at the top
of the figure and long period (low frequency, deep sensitivity)
is at the bottom -- matching a geological cross-section with
depth increasing downward. ``True`` flips the y-axis.
frame_pct : (lo, hi) or None, default ``(1, 99)``
Robust percentiles of finite log-period/log-frequency values used for
the visible y-frame. This prevents a few isolated frequency rows from
creating large empty bands. Pass ``None`` to show the absolute range,
or use *period_range* for explicit scientific limits.
scale : float, default ``0.85``
Fraction of each cell occupied by the *reference* ellipse
(the one with s1 = s1_ref). Values above 1 cause overlap.
normalise_by : ``"shape"`` | ``"cell"`` | ``"unity"`` | ``"abs"``
Sizing strategy:
``"shape"``
MTpy-style display: every major axis has the same physical size
and the minor/major ratio is ``s2 / s1``. This emphasizes tensor
shape without allowing absolute phase magnitude or unequal plot
axes to turn ellipses into misleading needles.
``"cell"``
Sizes are normalised so the 90th-percentile s1 fills *scale*
of its cell. Preserves relative size information.
``"unity"``
s1_ref = 1.0 (tan 45° = 1-D half-space reference). A 45°
phase tensor is drawn as a circle filling *scale* of its cell.
``"abs"``
Raw: width = scale × s1, height = scale × s2 in data units
(legacy behaviour).
s1_ref : float or None
Override the reference s1 for ``"cell"`` and ``"unity"`` modes.
min_aspect : float, default ``0.18``
Minimum height/width ratio enforced on every ellipse so that
near-degenerate (φ_min ≈ 0) cells stay visible as ovals instead
of collapsing into an invisible hairline. Set to ``0`` to plot
raw ellipticity.
c_by : str, default ``"skew"``
Column name or derived quantity to map to fill colour. Supported
values: ``"skew"``, ``"beta"``, ``"alpha"``, ``"theta"``,
``"ellipt"``, ``"s1"``, ``"s2"``, ``"|skew|"``, ``"|beta|"``,
``"|theta|"``, ``"phi_mean"``, ``"phi_max"``, ``"phi_min"``.
cmap : str, default ``"RdBu_r"``
Matplotlib colormap name.
clim : (vmin, vmax) or None
Explicit colour limits. When ``None``, derived from *clim_pct*.
clim_pct : (lo, hi), default ``(5.0, 95.0)``
Percentile limits used when *clim* is ``None``.
color_mode : ``"continuous"`` | ``"segmented"``
Continuous colour interpolation or three discrete classes. Segmented
mode follows the usual MTpy-style skew presentation.
segment_bounds : (lower, upper) or None
Class boundaries for segmented colour. By default these are
``(-skew_threshold, skew_threshold)``.
segment_colors : sequence of three colours
Colours below, inside, and above *segment_bounds*, respectively.
symmetric_clim : bool, default ``True``
Enforce vmin = −vmax for skew-like quantities.
edgecolor : str, default ``"k"``
Ellipse border colour. Set to ``"none"`` to suppress borders.
linewidth : float, default ``0.2``
Ellipse border width (pts). 3-D cells receive 3 × this width when
*mark_3d* is ``True``.
alpha : float, default ``0.92``
Ellipse fill opacity.
ellipse_kws : dict or None
Additional ellipse styling, such as ``edgecolor`` and ``linewidth``.
Geometry and transform keys remain controlled by this function.
cb_kws : dict or None
Colourbar customization. ``size``, ``pad``, ``labelsize`` and
``ticksize`` control layout; other entries go to ``Figure.colorbar``.
colorbar : bool, default ``True``
Draw the per-axes colourbar. Set ``False`` when composing several
calls onto a shared grid (e.g. one panel per line) and drawing a
single shared colourbar separately, using the same *cmap*/*clim*
passed to every call so the mapping is identical across panels.
skew_threshold : float or None, default ``3.0``
|β| threshold (degrees) separating 1-D/2-D from 3-D structure.
It does not set the color limits; *clim_pct* remains data-driven.
The threshold only highlights 3-D cells with a thicker border when
*mark_3d* is ``True``.
mark_3d : bool, default ``True``
Draw a thicker border on ellipses where |β| > *skew_threshold*.
ref_ellipse : bool, default ``True``
Draw a labelled reference circle (φ_max = φ_min = s1_ref, β = 0°)
inside the axes, without expanding the data frame.
legend_fontsize : float, default ``8.0``
Font size for the reference-circle label and the 1-D/2-D vs 3-D
annotation shown when *ref_ellipse* / *skew_threshold* are active.
annotations : bool, default ``True``
Draw the in-figure legend strip (the |β| 1-D/2-D vs 3-D labels and
the size-reference ellipse). ``False`` hides it and gives the
ellipses the full height.
title, xlabel, ylabel : str
Axes title and axis labels. Sensible defaults are used when empty.
tick_label_rotation : float, default ``45.0``
Station-name tick rotation in degrees.
figsize : (width, height), default ``(10.0, 5.5)``
Figure size (ignored when *ax* is provided).
recursive, on_dup, strict, verbose : see :func:`ensure_sites`.
ax : :class:`~matplotlib.axes.Axes` or None
If provided, draw into this axes and return it; otherwise create a
new figure.
Returns
-------
ax : :class:`~matplotlib.axes.Axes`
See Also
--------
build_phase_tensor_table : Underlying data computation.
plot_phase_tensor_summary : 3-panel figure combining ellipses,
dimensionality, and skew distribution.
plot_phase_tensor_map : Geographic map view.
"""
if isinstance(sites, pd.DataFrame):
df = sites
else:
df = build_phase_tensor_table(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
# ── resolve visual style from PYCSAMT_STYLE.pt_ellipse ──────────────
_es = PYCSAMT_STYLE.pt_ellipse
if scale is _UNSET:
scale = _es.scale
if normalise_by is _UNSET:
normalise_by = _es.normalise_by
if min_aspect is _UNSET:
min_aspect = _es.min_aspect
if c_by is _UNSET:
c_by = _es.c_by
if cmap is _UNSET:
cmap = _es.copy(c_by=c_by).resolve_cmap()
if clim_pct is _UNSET:
clim_pct = _es.clim_pct
if symmetric_clim is _UNSET:
symmetric_clim = _es.copy(c_by=c_by).resolve_symmetric_clim()
if edgecolor is _UNSET:
edgecolor = _es.edgecolor
if linewidth is _UNSET:
linewidth = _es.linewidth
if alpha is _UNSET:
alpha = _es.alpha
if skew_threshold is _UNSET:
skew_threshold = _es.skew_threshold
if mark_3d is _UNSET:
mark_3d = _es.mark_3d
if ref_ellipse is _UNSET:
ref_ellipse = _es.show_ref
# ── data selection ────────────────────────────────────────────────────
if not df.empty and stations is not None:
df = df[df["station"].isin(stations)]
if not df.empty and period_range is not None:
lo, hi = float(period_range[0]), float(period_range[1])
df = df[(df["period"] >= lo) & (df["period"] <= hi)]
if ax is None:
_, ax = plt.subplots(figsize=figsize)
if df.empty:
ax.text(
0.5,
0.5,
"no phase tensor data",
ha="center",
va="center",
transform=ax.transAxes,
)
return ax
# ── station ordering & y-axis ─────────────────────────────────────────
st_list = list(dict.fromkeys(df["station"].tolist())) # preserve order
st_map: dict[str, int] = {s: i for i, s in enumerate(st_list)}
n_st = len(st_list)
if axis_y == "logperiod":
df = df.copy()
df["_y"] = np.log10(df["period"].to_numpy(float))
ylab = ylabel or r"$\log_{10}(T)$ (s)"
else:
df = df.copy()
df["_y"] = np.log10(np.maximum(df["freq"].to_numpy(float), 1e-30))
ylab = ylabel or r"$\log_{10}(f)$ (Hz)"
# ── cell dimensions in data coordinates ───────────────────────────────
y_unique = np.unique(df["_y"].to_numpy(float))
cell_y = (
float(np.nanmedian(np.abs(np.diff(y_unique))))
if len(y_unique) > 1
else 1.0
)
cell_x = 1.0 # one station spacing
# ── ellipse-size normalisation reference ──────────────────────────────
s1_vals = df["s1"].to_numpy(float)
if normalise_by == "shape":
ref = None
elif normalise_by == "cell":
if s1_ref is not None:
ref = float(s1_ref)
else:
finite_s1 = s1_vals[np.isfinite(s1_vals)]
ref = (
float(np.nanpercentile(finite_s1, 90))
if len(finite_s1)
else 1.0
)
ref = max(ref, 1e-9)
elif normalise_by == "unity":
ref = float(s1_ref) if s1_ref is not None else 1.0
ref = max(ref, 1e-9)
else: # "abs" — legacy: w = scale × s1 in raw data units
ref = None
# ── colour values & limits ────────────────────────────────────────────
cvals, cbar_label = _resolve_cvals(df, c_by)
finite_c = cvals[np.isfinite(cvals)]
if clim is not None:
vmin, vmax = float(clim[0]), float(clim[1])
else:
lo_pct, hi_pct = clim_pct
vmin = (
float(np.nanpercentile(finite_c, lo_pct))
if len(finite_c)
else -1.0
)
vmax = (
float(np.nanpercentile(finite_c, hi_pct)) if len(finite_c) else 1.0
)
if symmetric_clim and c_by in _SYMMETRIC_C:
vlim = max(abs(vmin), abs(vmax))
vmin, vmax = -vlim, vlim
if not np.isfinite(vmin) or not np.isfinite(vmax):
vmin, vmax = -1.0, 1.0
elif np.isclose(vmin, vmax):
pad = max(abs(vmin) * 0.05, 1e-6)
vmin, vmax = vmin - pad, vmax + pad
color_mode = str(color_mode).lower()
if color_mode not in {"continuous", "segmented"}:
raise ValueError("color_mode must be 'continuous' or 'segmented'")
segment_ticks = segment_ticklabels = None
if color_mode == "segmented":
if len(segment_colors) != 3:
raise ValueError("segment_colors must contain exactly three colours")
if segment_bounds is None:
if skew_threshold is not None:
segment_bounds = (-float(skew_threshold), float(skew_threshold))
else:
step = (vmax - vmin) / 3.0
segment_bounds = (vmin + step, vmax - step)
lower, upper = map(float, segment_bounds)
if not lower < upper:
raise ValueError("segment_bounds must satisfy lower < upper")
outer_lo = min(vmin, lower - max(abs(lower) * 0.05, 1e-6))
outer_hi = max(vmax, upper + max(abs(upper) * 0.05, 1e-6))
cm = ListedColormap(segment_colors, name="phase_tensor_segments")
norm = BoundaryNorm([outer_lo, lower, upper, outer_hi], cm.N, clip=True)
segment_ticks = [0.5 * (outer_lo + lower), 0.5 * (lower + upper),
0.5 * (upper + outer_hi)]
segment_ticklabels = [f"< {lower:g}", f"{lower:g} to {upper:g}",
f"> {upper:g}"]
else:
norm = Normalize(vmin=vmin, vmax=vmax)
cm = plt.get_cmap(cmap)
# ── extract arrays for the drawing loop ───────────────────────────────
x_arr = df["station"].map(st_map).to_numpy(float)
y_arr = df["_y"].to_numpy(float)
s1_arr = df["s1"].to_numpy(float)
s2_arr = df["s2"].to_numpy(float)
th_arr = df["theta"].to_numpy(float)
# skew values for mark_3d (always needed even if c_by != "skew")
skew_arr = (
df["skew"].to_numpy(float)
if "skew" in df.columns
else np.zeros(len(df))
)
# Resolve the finite-data frame before drawing so ``shape`` mode can use
# physical display units, independent of station/period axis aspect.
y_all = y_arr[np.isfinite(y_arr)]
if len(y_all) and frame_pct is not None:
frame_lo, frame_hi = map(float, frame_pct)
if not (0.0 <= frame_lo < frame_hi <= 100.0):
raise ValueError("frame_pct must satisfy 0 <= lo < hi <= 100")
y_lo, y_hi = map(float, np.nanpercentile(y_all, [frame_lo, frame_hi]))
else:
y_lo = float(np.nanmin(y_all)) if len(y_all) else 0.0
y_hi = float(np.nanmax(y_all)) if len(y_all) else 1.0
margin = max(0.5 * cell_y, 0.02 * (y_hi - y_lo + 1e-9))
target_ylim = (
(y_lo - margin, y_hi + margin)
if period_up
else (y_hi + margin, y_lo - margin)
)
ax.set_xlim(-0.6, n_st - 0.4)
ax.set_ylim(*target_ylim)
PYCSAMT_STATION_RENDERING.apply(
ax,
np.arange(n_st, dtype=float),
st_list,
preset="pseudosection",
xlim=(-0.6, n_st - 0.4),
)
sm = ScalarMappable(cmap=cm, norm=norm)
sm.set_array([])
if colorbar:
cb = _attach_cbar(ax, sm, cbar_label, **dict(cb_kws or {}))
if segment_ticks is not None:
cb.set_ticks(segment_ticks)
cb.set_ticklabels(segment_ticklabels)
ax.figure.canvas.draw()
station_px = abs(
ax.transData.transform((1.0, y_lo))[0]
- ax.transData.transform((0.0, y_lo))[0]
)
# ── draw ellipses ─────────────────────────────────────────────────────
ellipse_options = dict(ellipse_kws or {})
forbidden = {"xy", "width", "height", "angle", "transform"}
bad_keys = forbidden.intersection(ellipse_options)
if bad_keys:
raise ValueError(
"ellipse_kws cannot override geometry keys: "
+ ", ".join(sorted(bad_keys))
)
ellipse_options.setdefault("edgecolor", edgecolor)
ellipse_options.setdefault("alpha", alpha)
ellipse_options.setdefault("zorder", 3)
for xi, yi, s1i, s2i, thi, ci, beti in zip(
x_arr, y_arr, s1_arr, s2_arr, th_arr, cvals, skew_arr
):
if not np.all(np.isfinite([xi, yi, s1i, s2i, thi, ci])):
continue
is_3d = (
skew_threshold is not None
and mark_3d
and np.isfinite(beti)
and abs(beti) > skew_threshold
)
if normalise_by == "shape":
ratio = abs(s2i) / max(abs(s1i), 1e-12)
ratio = min(1.0, max(float(min_aspect), ratio))
x0, x1 = ax.get_xlim()
y0, y1 = ax.get_ylim()
x_frac = (xi - x0) / (x1 - x0)
y_frac = (yi - y0) / (y1 - y0)
width_px = scale * station_px
shape_transform = (
Affine2D()
.scale(0.5 * width_px, 0.5 * width_px * ratio)
.rotate_deg(thi)
.scale(1.0 / ax.bbox.width, 1.0 / ax.bbox.height)
.translate(x_frac, y_frac)
+ ax.transAxes
)
patch_kws = dict(ellipse_options)
patch_kws.setdefault("facecolor", cm(norm(ci)))
base_lw = float(patch_kws.get("linewidth", linewidth))
patch_kws["linewidth"] = base_lw * 3.0 if is_3d else base_lw
ellipse = Ellipse(
(0.0, 0.0),
width=2.0,
height=2.0,
angle=0.0,
transform=shape_transform,
**patch_kws,
)
ellipse.set_clip_path(ax.patch)
ax.add_patch(ellipse)
continue
if ref is not None:
w = min(scale * cell_x * s1i / ref, 0.99 * cell_x)
h = min(scale * cell_y * s2i / ref, 0.99 * cell_y)
else:
w = scale * s1i
h = scale * s2i
if min_aspect > 0:
h = max(h, min_aspect * w)
patch_kws = dict(ellipse_options)
patch_kws.setdefault("facecolor", cm(norm(ci)))
base_lw = float(patch_kws.get("linewidth", linewidth))
patch_kws["linewidth"] = base_lw * 3.0 if is_3d else base_lw
ax.add_patch(
Ellipse(
(xi, yi),
width=w,
height=h,
angle=thi,
**patch_kws,
)
)
# ── axes limits and ticks ─────────────────────────────────────────────
ax.set_xlim(-0.6, n_st - 0.4)
# Reserve a dedicated legend strip beyond the data so the reference
# ellipse and 1-D/2-D vs 3-D annotation never overlap real ellipses.
want_legend = bool(annotations) and (bool(ref_ellipse) or (
skew_threshold is not None and c_by in ("skew", "beta")
))
legend_depth = 0.0
if period_up:
top_edge, top_sign = y_hi, 1.0
else:
top_edge, top_sign = y_lo, -1.0
ax.set_ylim(*target_ylim)
# y-ticks: integer log10 values
y_int = np.arange(int(np.floor(y_lo)), int(np.ceil(y_hi)) + 1)
ax.set_yticks(y_int)
ax.set_yticklabels([str(int(v)) for v in y_int], fontsize=8)
ax.set_ylabel(ylab, fontsize=9)
if title:
ax.set_title(title, fontsize=10)
ax.grid(True, ls=":", lw=0.4, color="#cccccc", alpha=0.7, zorder=0)
# ── colorbar ─────────────────────────────────────────────────────────
# ── optional overlays: dedicated legend strip beyond the data ─────────
if want_legend:
# boundary rule separating data from the legend strip
bound = top_edge + top_sign * margin
ax.plot(
[-0.6, n_st - 0.4], [bound, bound], color="0.82", lw=0.7, zorder=1
)
va_out = "bottom" if top_sign > 0 else "top"
if skew_threshold is not None and c_by in ("skew", "beta"):
y_near = 0.91
y_far = 0.975
ax.text(
0.01,
y_far,
f"|β| < {skew_threshold:.0f}° → 1-D/2-D",
fontsize=legend_fontsize,
va="center",
ha="left",
color="#2166ac",
transform=ax.transAxes,
bbox=dict(fc="white", ec="0.75", lw=0.5, alpha=0.9, pad=2.0),
)
ax.text(
0.01,
y_near,
f"|β| ≥ {skew_threshold:.0f}° → 3-D",
fontsize=legend_fontsize,
va="center",
ha="left",
color="#b2182b",
transform=ax.transAxes,
bbox=dict(fc="white", ec="0.75", lw=0.5, alpha=0.9, pad=2.0),
)
# Reference ellipse: the "1-D" size/shape used to normalise every
# data ellipse (φ_max = φ_min = ref), drawn at full data size so
# it is directly comparable to the ellipses in the plot.
if ref_ellipse and ref is not None:
rx, ry = 0.92, 0.94
rw, rh = 0.055, 0.03
ax.add_patch(
Ellipse(
(rx, ry),
width=rw,
height=rh,
angle=0.0,
facecolor="0.9",
edgecolor="k",
linewidth=0.9,
zorder=5,
transform=ax.transAxes,
)
)
ref_label = {
"unity": "1-D reference\n(Φ=45°)",
"cell": "size reference\n(90th-pct Φ_max)",
}.get(normalise_by, "size reference")
ax.text(
rx,
0.895,
ref_label,
ha="center",
va="top",
fontsize=legend_fontsize,
linespacing=1.3,
transform=ax.transAxes,
bbox=dict(fc="white", ec="0.75", lw=0.5, alpha=0.9, pad=2.0),
)
# Artists using ``ax.transAxes`` can still update Matplotlib's data limits
# when they are added as patches. Re-apply the finite-data frame last so
# annotations never create empty scientific-axis ranges.
ax.set_ylim(*target_ylim)
tick_lo = int(np.ceil(y_lo))
tick_hi = int(np.floor(y_hi))
if tick_lo <= tick_hi:
y_ticks = np.arange(tick_lo, tick_hi + 1)
else:
y_ticks = np.array([0.5 * (y_lo + y_hi)])
ax.set_yticks(y_ticks)
ax.set_yticklabels([f"{v:g}" for v in y_ticks], fontsize=8)
return ax
[docs]
def plot_phase_tensor_skewmap(
sites: Any,
*,
axis_y: str = "logperiod",
agg: str = "median",
figsize: tuple[float, float] = (9.0, 4.0),
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
ax: plt.Axes | None = None,
) -> plt.Axes:
df = build_phase_tensor_table(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
if df.empty:
if ax is None:
_, ax = plt.subplots(figsize=figsize)
ax.text(0.5, 0.5, "no phase tensor", ha="center", va="center")
return ax
if axis_y == "logperiod":
df = df.copy()
df["y"] = np.log10(df["period"].to_numpy())
ylab = LOG10_PERIOD_LABEL
else:
df = df.copy()
df["y"] = df["period"].to_numpy()
ylab = PERIOD_LABEL
piv = df.pivot_table(
index="y",
columns="station",
values="skew",
aggfunc=agg,
)
piv = piv.sort_index()
Z = piv.to_numpy(dtype=float) # (n_logp, n_st) — no transpose
if ax is None:
_, ax = plt.subplots(figsize=figsize)
im = ax.imshow(
Z,
aspect="auto",
origin="lower",
interpolation="nearest",
)
ax.set_ylabel(ylab)
PYCSAMT_STATION_RENDERING.apply(
ax,
np.arange(len(piv.columns), dtype=float),
list(piv.columns),
preset="pseudosection",
xlim=(-0.5, len(piv.columns) - 0.5),
)
yt = np.linspace(0, Z.shape[0] - 1, num=min(8, Z.shape[0]))
yvals = np.linspace(
piv.index.min(), piv.index.max(), num=min(8, len(piv.index))
)
ax.set_yticks(yt)
ax.set_yticklabels([f"{v:.3g}" for v in yvals])
if not ax.yaxis_inverted():
ax.invert_yaxis()
cb = plt.colorbar(im, ax=ax)
cb.set_label("skew")
return ax
# --- more plots ----------------------------------------------------------- #
def _coords_of_sites(sites) -> dict:
m = {}
for i, ed in enumerate(_iter_items(sites)):
nm = _name(ed, i)
lat = getattr(ed, "lat", None) or getattr(ed, "latitude", None)
lon = getattr(ed, "lon", None) or getattr(ed, "longitude", None)
try:
if lat is None or lon is None:
continue
m[nm] = (float(lat), float(lon))
except Exception:
continue
return m
[docs]
def plot_theta_vs_period(
sites: Any,
*,
figsize: tuple[float, float] = (8.0, 4.0),
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
ax: plt.Axes | None = None,
) -> plt.Axes:
df = build_phase_tensor_table(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
if ax is None:
_, ax = plt.subplots(figsize=figsize)
if df.empty:
ax.text(0.5, 0.5, "no phase tensor", ha="center", va="center")
return ax
for st, sdf in df.groupby("station"):
x = sdf["period"].to_numpy()
y = sdf["theta"].to_numpy()
ax.plot(x, y, ".", ms=3, label=st)
ax.set_xscale("log")
ax.set_xlabel("Period (s)")
ax.set_ylabel("theta (deg)")
ax.legend(ncol=2, fontsize=8)
return ax
def _strike_streamlines(ax, x, y, th, w, n_stations):
"""Overlay smoothed geoelectric-strike streamlines.
Strike is axial (mod 180 deg), so it is averaged through the
double-angle ``(cos 2theta, sin 2theta)`` (weighted by ``w``, the 2-D
strength) onto a regular grid, then halved back before integrating the
flow. Silently degrades to a no-op when SciPy is unavailable or the
grid is too sparse to integrate.
"""
if n_stations < 3 or y.size < 12:
return
try:
from scipy.interpolate import griddata
except Exception:
return
# Restrict the flow to the *populated* period band so it never
# extrapolates into empty short/long-period regions (e.g. long-period
# MT lines with no high-frequency data).
y_lo, y_hi = np.nanpercentile(y, [2.0, 98.0])
if not np.isfinite(y_lo) or y_hi - y_lo < 1e-6:
return
gx = np.arange(n_stations, dtype=float)
gy = np.linspace(y_lo, y_hi, 44)
GX, GY = np.meshgrid(gx, gy)
def _fill(val):
g = griddata((x, y), val, (GX, GY), method="linear")
nn = griddata((x, y), val, (GX, GY), method="nearest")
bad = ~np.isfinite(g)
g[bad] = nn[bad]
return g
c2 = _fill(w * np.cos(2.0 * th))
s2 = _fill(w * np.sin(2.0 * th))
if not (np.isfinite(c2).any() and np.isfinite(s2).any()):
return
thg = 0.5 * np.arctan2(s2, c2)
try:
ax.streamplot(
GX,
GY,
np.sin(thg),
np.cos(thg),
density=1.1,
color="0.2",
linewidth=0.7,
arrowstyle="-",
)
except Exception:
pass
[docs]
def plot_strike_director_field(
sites: Any,
*,
color_by: str = "skew",
length_by: str | None = "ellipt",
streamlines: bool = True,
skew_max: float = 6.0,
cmap: str | None = None,
period_subsample: int | None = None,
bar_scale: float = 26.0,
show_legend: bool = True,
title: str | None = None,
figsize: tuple[float, float] = (12.0, 5.2),
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
ax: plt.Axes | None = None,
) -> plt.Axes:
r"""Geoelectric-strike **director field** over station and period.
A supplement to :func:`plot_theta_vs_period`. Because the phase-tensor
strike ``theta`` is an *axial* angle (defined mod 180 deg), the correct
glyph is not a point on a linear axis but a **head-less bar** pointing
along the strike. This draws one director per ``(station, period)`` cell
on a station x log-period grid, encoding three more channels at once:
* **orientation** -- the strike ``theta``;
* **length** -- ``length_by`` (default ellipticity), the 2-D strength:
near-1-D cells get a short bar because there strike is ill-defined;
* **colour** -- ``color_by`` (default ``|skew|``), the departure from
2-D: green = low-skew / reliable, red = high-skew / 3-D or galvanic
distortion where the strike should not be trusted.
An optional smoothed **streamline** overlay (``streamlines=True``)
integrates the director field into a strike "flow", so lateral and
vertical coherence read at a glance.
Interpretation
--------------
* long, aligned, green bars flowing in a laminar bundle -> robust,
depth-consistent 2-D strike; read the azimuth with confidence;
* bars rotating smoothly with depth -> strike varies with depth
(dipping structure or layered anisotropy);
* short and/or red, swirling bars -> 1-D, 3-D, or noise: do not
over-interpret the direction there.
Parameters
----------
sites : path, EDI object, APISurvey, Sites, or iterable of sites
Anything accepted by :func:`build_phase_tensor_table`.
color_by : {'skew', 'ellipt', 's1', 's2', ...}, default 'skew'
Table column mapped to bar colour (its absolute value is used).
length_by : str or None, default 'ellipt'
Table column mapped to bar length (absolute value, 95th-percentile
normalised). ``None`` draws uniform-length bars.
streamlines : bool, default True
Overlay smoothed strike streamlines (needs SciPy).
skew_max : float, default 6.0
Upper clip of the ``|skew|`` colour scale, in degrees (only used
when ``color_by='skew'``). Skew above a few degrees already flags
3-D behaviour.
cmap : str, optional
Override the colour map (default ``'RdYlGn_r'`` for skew, else
``'viridis'``).
period_subsample : int, optional
Keep at most this many periods (evenly along the log axis) to thin
a very dense grid.
bar_scale : float, default 26.0
Matplotlib quiver ``scale`` -- larger makes shorter bars.
show_legend : bool, default True
Draw the director / streamline legend.
title : str, optional
Axes title.
figsize, recursive, on_dup, strict, verbose, ax
Standard ``emtools`` plotting arguments.
Returns
-------
matplotlib.axes.Axes
See Also
--------
plot_theta_vs_period : the linear scatter this supplements.
plot_phase_tensor_psection : per-cell ellipses (magnitudes as well).
"""
df = build_phase_tensor_table(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
if ax is None:
_, ax = plt.subplots(figsize=figsize)
fig = ax.get_figure()
if df.empty:
ax.text(0.5, 0.5, "no phase tensor", ha="center", va="center")
return ax
df = df.copy()
df["logp"] = np.log10(df["period"].to_numpy())
stations = list(dict.fromkeys(df["station"]))
sidx = {s: i for i, s in enumerate(stations)}
df["xi"] = df["station"].map(sidx).astype(float)
if period_subsample and period_subsample > 0:
periods = np.sort(df["period"].unique())
step = max(1, len(periods) // period_subsample)
keep = set(periods[::step])
df = df[df["period"].isin(keep)]
x = df["xi"].to_numpy(float)
y = df["logp"].to_numpy(float)
th = np.radians(df["theta"].to_numpy(float))
if length_by and length_by in df.columns:
lv = np.abs(df[length_by].to_numpy(float))
denom = np.nanpercentile(lv, 95)
denom = denom if np.isfinite(denom) and denom > 0 else 1.0
length = np.clip(lv / denom, 0.12, 1.0)
else:
length = np.ones(len(df))
if color_by == "skew":
cval = np.abs(df["skew"].to_numpy(float))
vmin, vmax = 0.0, float(skew_max)
cmap_name = cmap or "RdYlGn_r"
clabel = "|skew| (deg) - green: 2-D / reliable, red: 3-D / distorted"
elif color_by in df.columns:
cval = np.abs(df[color_by].to_numpy(float))
vmin = float(np.nanmin(cval))
vmax = float(np.nanpercentile(cval, 95))
cmap_name = cmap or "viridis"
clabel = str(color_by)
else:
cval = np.zeros(len(df))
vmin, vmax, cmap_name, clabel = 0.0, 1.0, cmap or "viridis", ""
ok = (
np.isfinite(x) & np.isfinite(y) & np.isfinite(th) & np.isfinite(length)
)
x, y, th, length, cval = x[ok], y[ok], th[ok], length[ok], cval[ok]
if x.size == 0:
ax.text(0.5, 0.5, "no finite strike", ha="center", va="center")
return ax
# north-up director components (east = sin, north = cos); axial, so
# (theta) and (theta+180) give the same head-less bar under pivot='mid'.
U = np.sin(th) * length
V = np.cos(th) * length
q = ax.quiver(
x,
y,
U,
V,
np.clip(cval, vmin, vmax),
cmap=cmap_name,
angles="uv",
pivot="mid",
headlength=0,
headaxislength=0,
headwidth=1,
scale=bar_scale,
width=0.0032,
alpha=0.9,
clim=(vmin, vmax),
)
if clabel:
cb = fig.colorbar(q, ax=ax, pad=0.01, shrink=0.9)
cb.set_label(clabel, fontsize=8)
if streamlines:
_strike_streamlines(ax, x, y, th, length, len(stations))
ax.set_ylabel(LOG10_PERIOD_LABEL)
PYCSAMT_STATION_RENDERING.apply(
ax,
np.arange(len(stations)),
stations,
preset="pseudosection",
xlim=(-0.8, len(stations) - 0.2),
)
if not ax.yaxis_inverted():
ax.invert_yaxis()
ax.set_title(
title if title is not None else "Geoelectric strike director field",
fontsize=10,
)
if show_legend:
from matplotlib.lines import Line2D
handles = [
Line2D(
[0],
[0],
color="#4a7a1f",
lw=3,
label="strike director (length = 2-D strength)",
),
]
if streamlines:
handles.append(
Line2D(
[0],
[0],
color="0.2",
lw=1.2,
label="smoothed strike flow",
),
)
ax.legend(
handles=handles, loc="upper right", fontsize=7, framealpha=0.85
)
return ax
[docs]
def plot_ellipticity_psection(
sites: Any,
*,
figsize: tuple[float, float] = (8.5, 4.0),
agg: str = "median",
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
ax: plt.Axes | None = None,
) -> plt.Axes:
df = build_phase_tensor_table(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
if ax is None:
_, ax = plt.subplots(figsize=figsize)
if df.empty:
ax.text(0.5, 0.5, "no phase tensor", ha="center", va="center")
return ax
df = df.copy()
df["logp"] = np.log10(df["period"].to_numpy())
piv = df.pivot_table(
index="logp",
columns="station",
values="ellipt",
aggfunc=agg,
)
piv = piv.sort_index()
Z = piv.to_numpy(dtype=float) # (n_logp, n_st) — no transpose
im = ax.imshow(
Z,
aspect="auto",
origin="lower",
interpolation="nearest",
)
ax.set_ylabel(LOG10_PERIOD_LABEL)
PYCSAMT_STATION_RENDERING.apply(
ax,
np.arange(Z.shape[1], dtype=float),
list(piv.columns),
preset="pseudosection",
xlim=(-0.5, Z.shape[1] - 0.5),
)
yt = np.linspace(
0, Z.shape[0] - 1, num=min(8, Z.shape[0])
) # shape[0] = n_logp
yvals = np.linspace(
piv.index.min(), piv.index.max(), num=min(8, len(piv.index))
)
ax.set_yticks(yt)
ax.set_yticklabels([f"{v:.2g}" for v in yvals])
if not ax.yaxis_inverted():
ax.invert_yaxis()
cb = plt.colorbar(im, ax=ax)
cb.set_label("ellipticity")
return ax
#: Traffic-light dimensionality palette (1D=safe green, 2D=caution amber,
#: 3D=warning red), matching :func:`pycsamt.emtools.skew.plot_skew_traffic_psection`.
DIMENSIONALITY_TRAFFIC_COLORS = (
(0.20, 0.60, 0.20), # 1D — green
(0.95, 0.70, 0.20), # 2D — amber
(0.85, 0.25, 0.20), # 3D — red
)
[docs]
def plot_dimensionality_psection(
sites: Any,
*,
skew_th: float = 3.0,
ellipt_th: float = 0.2,
segmented_colors: bool = True,
cmap: str | Any = "viridis",
colorbar: bool = True,
figsize: tuple[float, float] = (8.5, 4.0),
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
ax: plt.Axes | None = None,
) -> plt.Axes:
df = build_phase_tensor_table(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
if ax is None:
_, ax = plt.subplots(figsize=figsize)
if df.empty:
ax.text(0.5, 0.5, "no phase tensor", ha="center", va="center")
return ax
df = df.copy()
df["logp"] = np.log10(df["period"].to_numpy())
# 0=1D, 1=2D, 2=3D (simple rule)
lab = np.zeros(len(df), dtype=int)
lab[
(np.abs(df["skew"]) <= skew_th) & (df["ellipt"].abs() <= ellipt_th)
] = 0
lab[(np.abs(df["skew"]) <= skew_th) & (df["ellipt"].abs() > ellipt_th)] = 1
lab[(np.abs(df["skew"]) > skew_th)] = 2
df["dim"] = lab
piv = df.pivot_table(
index="logp",
columns="station",
values="dim",
aggfunc="median",
)
piv = piv.sort_index()
Z = piv.to_numpy(dtype=float) # (n_logp, n_st) — no transpose
if segmented_colors:
plot_cmap = ListedColormap(DIMENSIONALITY_TRAFFIC_COLORS)
plot_norm = BoundaryNorm([-0.5, 0.5, 1.5, 2.5], plot_cmap.N)
else:
plot_cmap = cmap
plot_norm = None
im = ax.imshow(
Z,
aspect="auto",
origin="lower",
interpolation="nearest",
cmap=plot_cmap,
norm=plot_norm,
)
ax.set_ylabel(LOG10_PERIOD_LABEL)
PYCSAMT_STATION_RENDERING.apply(
ax,
np.arange(Z.shape[1]),
list(piv.columns),
preset="pseudosection",
xlim=(-0.5, Z.shape[1] - 0.5),
)
if not ax.yaxis_inverted():
ax.invert_yaxis()
yt = np.linspace(
0, Z.shape[0] - 1, num=min(8, Z.shape[0])
) # shape[0] = n_logp
yvals = np.linspace(
piv.index.min(), piv.index.max(), num=min(8, len(piv.index))
)
ax.set_yticks(yt)
ax.set_yticklabels([f"{v:.2g}" for v in yvals])
if colorbar:
cb = plt.colorbar(im, ax=ax)
if segmented_colors:
cb.set_ticks([0, 1, 2])
cb.set_ticklabels(["1D", "2D", "3D"])
cb.set_label("Dimensionality")
else:
cb.set_label("dim (0=1D,1=2D,2=3D)")
return ax
[docs]
def plot_phase_tensor_rose(
sites: Any,
*,
# ── visual style ─────────────────────────────────────────────────────
style: str | RoseStyle | None = "pycsamt",
# ── data selection ───────────────────────────────────────────────────
band: tuple[float, float] | None = None,
freq_bands: list[tuple[float, float]] | None = None,
band_labels: list[str] | None = None,
band_colors: list | None = None,
# ── histogram ────────────────────────────────────────────────────────
bins: int = 36,
# ── visual overrides (sentinel → taken from *style*) ─────────────────
bar_style=_UNSET,
bar_color=_UNSET,
bar_alpha=_UNSET,
bar_edgecolor=_UNSET,
bar_edgelw=_UNSET,
cmap=_UNSET,
outer_ring_lw=_UNSET,
outer_ring_color=_UNSET,
n_rings=_UNSET,
ring_color=_UNSET,
ring_ls=_UNSET,
ring_lw=_UNSET,
ring_labels=_UNSET,
ring_label_angle=_UNSET,
ring_label_fontsize=_UNSET,
ring_label_color=_UNSET,
ring_label_fmt=_UNSET,
spoke_every=_UNSET,
spoke_color=_UNSET,
spoke_ls=_UNSET,
spoke_lw=_UNSET,
compass_labels=_UNSET,
compass_fontsize=_UNSET,
compass_color=_UNSET,
compass_fontweight=_UNSET,
show_mean=_UNSET,
mean_color=_UNSET,
mean_lw=_UNSET,
mean_ls=_UNSET,
show_secondary=_UNSET,
secondary_ls=_UNSET,
secondary_lw=_UNSET,
secondary_color=_UNSET,
show_annotation=_UNSET,
annotation_pos=_UNSET,
annotation_fontsize=_UNSET,
annotation_bg=_UNSET,
annotation_ec=_UNSET,
show_n=_UNSET,
# ── layout ───────────────────────────────────────────────────────────
figsize: tuple[float, float] = (5.5, 5.5),
title: str = "",
title_fontsize: float = 10.0,
# ── core ─────────────────────────────────────────────────────────────
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
ax: plt.Axes | None = None,
) -> plt.Axes:
"""Publication-quality phase-tensor θ rose diagram.
Bars are drawn with axial symmetry — 0–180° mirrored to 180–360° —
reflecting the inherent 180° ambiguity of the phase-tensor principal
axis direction.
Parameters
----------
sites : any
EDI paths, objects, or collection accepted by
:func:`~pycsamt.emtools._core.ensure_sites`.
band : (lo_s, hi_s) or None
Period window in seconds. ``None`` uses all available periods.
freq_bands : list of (lo_s, hi_s), optional
Sub-bands for ``bar_style="bands"`` — one stacked bar colour per
band.
band_labels : list[str], optional
Legend labels for each entry in *freq_bands*.
band_colors : list, optional
Colours for each *freq_bands* entry. Defaults to ``tab10``.
bins : int
Number of bins over 0–180°, mirrored to 360°. Default ``36``
gives 5° bins.
bar_style : {"gradient", "solid", "bands"}
``"gradient"`` colours bars by height via *cmap*;
``"solid"`` uses *bar_color* uniformly;
``"bands"`` stacks one colour per period sub-band.
bar_color : str
Bar fill colour for ``bar_style="solid"``.
bar_edgecolor : str
Bar edge colour (``"none"`` → no edge).
bar_edgelw : float
Bar edge line-width.
bar_alpha : float
Bar opacity (0–1).
cmap : str
Colormap for ``bar_style="gradient"``.
outer_ring_lw : float
Line-width of the bold outer bounding circle.
outer_ring_color : str
Colour of the outer circle.
n_rings : int
Number of concentric reference rings.
ring_color : str
Colour of concentric rings and spokes.
ring_ls : str
Line-style of concentric rings.
ring_lw : float
Line-width of concentric rings.
ring_labels : list[float], optional
Explicit count values to annotate on the rings (e.g.
``[25, 50, 75, 100]``). ``None`` → evenly spaced from
``rmax / n_rings`` to ``rmax``.
ring_label_angle : float
Clockwise angle from North (degrees) at which ring count labels
are placed. Default ``22.5``.
ring_label_fontsize : float
Font size for ring count labels.
ring_label_color : str
Colour for ring count labels.
ring_label_fmt : str
Format string for ring labels, e.g. ``"{:.0f}"``.
spoke_every : float
Angular spacing (degrees) between radial spokes.
spoke_color : str
Colour of radial spokes.
spoke_ls : str
Line-style of radial spokes.
spoke_lw : float
Line-width of radial spokes.
compass_labels : {"NESW", "degrees", "none"}
Perimeter labels. ``"NESW"`` shows cardinal directions;
``"degrees"`` shows degree values; ``"none"`` hides all.
compass_fontsize : float
Font size for perimeter labels.
compass_color : str
Colour for perimeter labels.
compass_fontweight : str
Font weight for perimeter labels.
show_mean : bool
Draw the axial mean direction as a line through the centre.
mean_color : str
Colour of the mean-direction line.
mean_lw : float
Line-width of the mean-direction line.
mean_ls : str
Line-style of the mean-direction line.
show_secondary : bool
Draw the 180°-conjugate mean line.
secondary_ls : str
Line-style of the conjugate line.
secondary_lw : float, optional
Line-width of the conjugate line; defaults to *mean_lw*.
secondary_color : str, optional
Colour of the conjugate line; defaults to *mean_color*.
show_annotation : bool
Show a text box with the mean θ and station/pair count.
annotation_pos : (float, float)
Axes-fraction ``(x, y)`` for the annotation box.
annotation_fontsize : float
Font size of the annotation text.
annotation_bg : str
Background colour of the annotation box.
annotation_ec : str
Edge colour of the annotation box.
show_n : bool
Append ``n = N`` to the annotation text.
figsize : (float, float)
Figure size in inches.
title : str
Axes title (set via ``ax.set_title``).
title_fontsize : float
Font size for *title*.
recursive, on_dup, strict, verbose
Passed to :func:`~pycsamt.emtools._core.ensure_sites`.
ax : matplotlib.axes.Axes, optional
Pre-existing polar axes to draw into. Created when ``None``.
Returns
-------
matplotlib.axes.Axes
The polar axes containing the rose diagram.
Examples
--------
Default gradient style, all periods:
>>> from pycsamt.emtools import plot_phase_tensor_rose
>>> ax = plot_phase_tensor_rose("path/to/edis/", figsize=(6, 6))
Frequency-band decomposition (stacked):
>>> ax = plot_phase_tensor_rose(
... sites,
... bar_style="bands",
... freq_bands=[(1e-4, 1e-2), (1e-2, 1e0)],
... band_labels=["Short period", "Long period"],
... )
Custom ring count labels:
>>> ax = plot_phase_tensor_rose(
... sites,
... ring_labels=[25, 50, 75, 100],
... ring_label_angle=15.0,
... )
"""
# ── resolve style → fill _UNSET visual params ─────────────────────────
rs = resolve_rose_style(style)
def _v(val, attr):
return getattr(rs, attr) if val is _UNSET else val
bar_style = _v(bar_style, "bar_style")
bar_color = _v(bar_color, "bar_color")
bar_alpha = _v(bar_alpha, "bar_alpha")
bar_edgecolor = _v(bar_edgecolor, "bar_edgecolor")
bar_edgelw = _v(bar_edgelw, "bar_edgelw")
cmap = _v(cmap, "cmap")
outer_ring_lw = _v(outer_ring_lw, "outer_ring_lw")
outer_ring_color = _v(outer_ring_color, "outer_ring_color")
n_rings = _v(n_rings, "n_rings")
ring_color = _v(ring_color, "ring_color")
ring_ls = _v(ring_ls, "ring_ls")
ring_lw = _v(ring_lw, "ring_lw")
ring_labels = _v(ring_labels, "ring_labels")
ring_label_angle = _v(ring_label_angle, "ring_label_angle")
ring_label_fontsize = _v(ring_label_fontsize, "ring_label_fontsize")
ring_label_color = _v(ring_label_color, "ring_label_color")
ring_label_fmt = _v(ring_label_fmt, "ring_label_fmt")
spoke_every = _v(spoke_every, "spoke_every")
spoke_color = _v(spoke_color, "spoke_color")
spoke_ls = _v(spoke_ls, "spoke_ls")
spoke_lw = _v(spoke_lw, "spoke_lw")
compass_labels = _v(compass_labels, "compass_labels")
compass_fontsize = _v(compass_fontsize, "compass_fontsize")
compass_color = _v(compass_color, "compass_color")
compass_fontweight = _v(compass_fontweight, "compass_fontweight")
show_mean = _v(show_mean, "show_mean")
mean_color = _v(mean_color, "mean_color")
mean_lw = _v(mean_lw, "mean_lw")
mean_ls = _v(mean_ls, "mean_ls")
show_secondary = _v(show_secondary, "show_secondary")
secondary_color = _v(secondary_color, "secondary_color")
secondary_ls = _v(secondary_ls, "secondary_ls")
secondary_lw = _v(secondary_lw, "secondary_lw")
show_annotation = _v(show_annotation, "show_annotation")
annotation_pos = _v(annotation_pos, "annotation_pos")
annotation_fontsize = _v(annotation_fontsize, "annotation_fontsize")
annotation_bg = _v(annotation_bg, "annotation_bg")
annotation_ec = _v(annotation_ec, "annotation_ec")
show_n = _v(show_n, "show_n")
# ── axial mean helper (local, avoids import cycle with strike.py) ──────
def _axial_mean(a_deg: np.ndarray) -> float:
"""Circular mean for axial (0–180°) data, returned in [0,180°)."""
rad = np.radians(2.0 * a_deg)
mu = (
np.degrees(
np.arctan2(np.nanmean(np.sin(rad)), np.nanmean(np.cos(rad)))
)
/ 2.0
)
return float(mu % 180.0)
# ── build phase-tensor table ───────────────────────────────────────────
df = build_phase_tensor_table(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
if ax is None:
_, ax = plt.subplots(subplot_kw=dict(polar=True), figsize=figsize)
if df.empty:
ax.text(
0.5,
0.5,
"no phase tensor data",
ha="center",
va="center",
transform=ax.transAxes,
)
return ax
# ── period filtering ───────────────────────────────────────────────────
if band is not None:
lo, hi = float(band[0]), float(band[1])
sel = (df["period"] >= lo) & (df["period"] <= hi)
band_label = f"{lo:.4g}–{hi:.4g} s"
else:
sel = np.ones(len(df), dtype=bool)
p = df["period"]
band_label = f"{p.min():.4g}–{p.max():.4g} s"
# ── histogram over 0–180° ─────────────────────────────────────────────
bins_ = int(max(12, bins))
edges_deg = np.linspace(0.0, 180.0, bins_ + 1)
dw = np.radians(180.0 / bins_)
use_bands = bar_style == "bands" and bool(freq_bands)
if use_bands:
n_fb = len(freq_bands) # type: ignore[arg-type]
_bc: list = (
list(band_colors)
if band_colors is not None
else list(plt.get_cmap("tab10")(np.linspace(0, 0.8, n_fb)))
)
_bl: list[str] = (
list(band_labels)
if band_labels is not None
else [f"{lo_:.4g}–{hi_:.4g} s" for lo_, hi_ in freq_bands]
) # type: ignore
hists: list[np.ndarray] = []
for fb in freq_bands: # type: ignore[union-attr]
lo_, hi_ = float(fb[0]), float(fb[1])
sel_fb = (df["period"] >= lo_) & (df["period"] <= hi_)
th_fb = df.loc[sel_fb, "theta"].to_numpy(float) % 180.0
h, _ = np.histogram(th_fb, bins=edges_deg)
hists.append(h)
hist_total = np.sum(hists, axis=0)
else:
th_raw = df.loc[sel, "theta"].to_numpy(float) % 180.0
hist_total, _ = np.histogram(th_raw, bins=edges_deg)
rmax = float(hist_total.max()) if hist_total.max() > 0 else 1.0
ang_mid = np.radians(0.5 * (edges_deg[1:] + edges_deg[:-1]))
# ── polar setup ───────────────────────────────────────────────────────
ax.set_theta_zero_location("N")
ax.set_theta_direction(-1)
# hide radial tick labels; keep custom rings without numeric radius text
ax.set_yticklabels([])
ax.set_yticks([])
hide_polar_radius_labels(ax)
ax.yaxis.grid(False)
ax.xaxis.grid(False)
ax.set_frame_on(False)
# ── draw bars (0–180° and mirror 180–360°) ────────────────────────────
def _draw_bars(
ang_c: np.ndarray, heights: np.ndarray, color, bottom: float = 0.0
) -> None:
for a, h in zip(ang_c, heights):
for a_plot in (a, a + np.pi): # axial symmetry
ax.bar(
a_plot,
h,
width=dw,
bottom=bottom,
color=color,
edgecolor=bar_edgecolor,
linewidth=bar_edgelw,
alpha=bar_alpha,
)
if use_bands:
bottoms = np.zeros(bins_)
for h, col in zip(hists, _bc):
_draw_bars(ang_mid, h, col, bottom=0.0)
# stacked: accumulate actual bottom for next band
bottoms += h
# legend
from matplotlib.patches import Patch as _Patch
ax.legend(
handles=[
_Patch(fc=c, label=l, alpha=bar_alpha)
for c, l in zip(_bc, _bl)
],
loc="lower left",
bbox_to_anchor=(1.05, 0.0),
fontsize=max(6, annotation_fontsize - 1),
framealpha=0.9,
)
elif bar_style == "gradient":
cm_ = plt.get_cmap(cmap)
for a, h in zip(ang_mid, hist_total):
col = cm_(h / rmax) if rmax > 0 else cm_(0.5)
for a_plot in (a, a + np.pi):
ax.bar(
a_plot,
h,
width=dw,
color=col,
edgecolor=bar_edgecolor,
linewidth=bar_edgelw,
alpha=bar_alpha,
)
else: # solid
_draw_bars(ang_mid, hist_total, bar_color)
# ── concentric rings ──────────────────────────────────────────────────
if ring_labels is not None:
r_levels = [float(v) for v in ring_labels]
else:
step = rmax / n_rings
r_levels = [step * k for k in range(1, n_rings + 1)]
theta_full = np.linspace(0, 2 * np.pi, 360)
for rv in r_levels:
ax.plot(
theta_full,
np.full_like(theta_full, rv),
color=ring_color,
ls=ring_ls,
lw=ring_lw,
zorder=0,
)
# ── radial spokes ─────────────────────────────────────────────────────
n_spokes = int(round(360.0 / spoke_every))
spoke_angles = np.radians(np.arange(n_spokes) * spoke_every)
for sa in spoke_angles:
ax.plot(
[sa, sa],
[0, rmax],
color=spoke_color,
ls=spoke_ls,
lw=spoke_lw,
zorder=0,
)
# ── outer ring ────────────────────────────────────────────────────────
ax.spines["polar"].set_linewidth(outer_ring_lw)
ax.spines["polar"].set_edgecolor(outer_ring_color)
ax.set_ylim(0, rmax * 1.08)
# ── compass labels ────────────────────────────────────────────────────
n_spk = int(round(360.0 / spoke_every))
spoke_degs = np.arange(n_spk) * spoke_every
if compass_labels == "NESW":
_card = {0: "N", 90: "E", 180: "S", 270: "W"}
lbl_list = [_card.get(int(d) % 360, "") for d in spoke_degs]
elif compass_labels == "degrees":
lbl_list = [f"{int(d)}°" for d in spoke_degs]
else:
lbl_list = [""] * n_spk
ax.set_thetagrids(
spoke_degs,
labels=lbl_list,
fontsize=compass_fontsize,
)
for lbl in ax.get_xticklabels():
lbl.set_color(compass_color)
lbl.set_fontweight(compass_fontweight)
# ── mean direction ────────────────────────────────────────────────────
th_all = df.loc[sel, "theta"].to_numpy(float) % 180.0
mu = _axial_mean(th_all) if len(th_all) > 0 else 0.0
sec_lw = secondary_lw if secondary_lw is not None else mean_lw
sec_col = secondary_color if secondary_color is not None else mean_color
if show_mean and len(th_all) > 0:
mu_r = np.radians(mu)
# draw full diameter: two radii (center→rim) in opposite directions
for ang_pt in (mu_r, mu_r + np.pi):
ax.plot(
[ang_pt, ang_pt],
[0, rmax],
color=mean_color,
lw=mean_lw,
ls=mean_ls,
solid_capstyle="round",
zorder=10,
)
if show_secondary:
mu90 = np.radians(mu + 90.0)
for ang_pt in (mu90, mu90 + np.pi):
ax.plot(
[ang_pt, ang_pt],
[0, rmax],
color=sec_col,
lw=sec_lw,
ls=secondary_ls,
solid_capstyle="round",
zorder=9,
)
# ── annotation box ────────────────────────────────────────────────────
if show_annotation and len(th_all) > 0:
txt = f"θ̄ = {mu:.1f}°"
if show_n:
txt += f"\nn = {len(th_all)}"
ax.text(
annotation_pos[0],
annotation_pos[1],
txt,
transform=ax.transAxes,
fontsize=annotation_fontsize,
va="top",
ha="left",
bbox=dict(
boxstyle="round,pad=0.3",
fc=annotation_bg,
ec=annotation_ec,
lw=0.8,
),
zorder=20,
)
# ── title ─────────────────────────────────────────────────────────────
if title:
ax.set_title(title, fontsize=title_fontsize, pad=14)
else:
ax.set_title(
f"Phase-tensor θ rose ({band_label})",
fontsize=title_fontsize,
pad=14,
)
return ax
def _topography_bg_grid(
topography: Any,
sites: Any,
coords: dict[str, tuple[float, float]],
*,
cmap: str = "terrain",
alpha: float = 0.55,
label: str = "Elevation (m)",
n: int = 90,
) -> dict[str, Any] | None:
"""Build a ``bg_grid`` dict from an elevation source.
Parameters
----------
topography : dict | bool | str | pathlib.Path
``dict`` — returned unchanged (already a ``bg_grid`` spec).
``True`` — station elevations read from the EDI headers
(``ed.coords[2]`` / ``ed.elev`` / ``ed.meta``).
path to a CSV — either ``lon,lat,elev`` rows (a scattered DEM)
or ``station,elev`` rows (joined to *coords*).
sites : any
The already-loaded sites/collection.
coords : dict[str, (lat, lon)]
Station coordinates in the map frame.
cmap, alpha, label, n :
Background colormap, opacity, colorbar label, and grid density.
Returns
-------
dict or None
``{"lons", "lats", "values", "cmap", "alpha", "label"}`` ready
for ``plot_phase_tensor_map(bg_grid=...)``, or ``None`` when
fewer than 3 finite elevation points are available.
"""
if isinstance(topography, dict):
return topography
tri: list[tuple[float, float, float]] = []
if topography is True:
for i, ed in enumerate(_iter_items(sites)):
st = _name(ed, i)
if st not in coords:
continue
z = None
c = getattr(ed, "coords", None)
if (
c is not None
and len(c) >= 3
and c[2] is not None
and np.isfinite(c[2])
):
z = float(c[2])
else:
e = getattr(ed, "elev", None)
if e is None:
e = (getattr(ed, "meta", {}) or {}).get("elev")
try:
if e is not None and np.isfinite(float(e)):
z = float(e)
except (TypeError, ValueError):
z = None
latlon = coords[st]
if z is None and len(latlon) >= 3 and np.isfinite(latlon[2]):
z = float(latlon[2]) # elevation carried in coords
if z is not None:
tri.append((float(latlon[1]), float(latlon[0]), z))
else:
try:
df = pd.read_csv(topography)
except Exception: # noqa: BLE001 — best effort
return None
low = {str(c).lower(): c for c in df.columns}
def _col(*names: str) -> Any:
return next((low[k] for k in names if k in low), None)
z_c = _col(
"elev", "elevation", "elevation_m", "z", "alt", "altitude",
"height",
)
lon_c = _col("lon", "longitude", "long", "x", "easting")
lat_c = _col("lat", "latitude", "y", "northing")
st_c = _col("station", "site", "id", "name", "sta")
if z_c is None:
return None
if lon_c is not None and lat_c is not None:
for _, r in df.iterrows():
try:
tri.append(
(float(r[lon_c]), float(r[lat_c]), float(r[z_c]))
)
except (TypeError, ValueError):
pass
elif st_c is not None:
emap: dict[str, float] = {}
for _, r in df.iterrows():
try:
emap[str(r[st_c])] = float(r[z_c])
except (TypeError, ValueError):
pass
for st, (lat, lon) in coords.items():
if st in emap:
tri.append((lon, lat, emap[st]))
pts = np.asarray(tri, dtype=float)
if pts.ndim != 2 or pts.shape[0] < 3:
return None
pts = pts[np.isfinite(pts).all(axis=1)]
if pts.shape[0] < 3:
return None
lon = np.linspace(pts[:, 0].min(), pts[:, 0].max(), n)
lat = np.linspace(pts[:, 1].min(), pts[:, 1].max(), n)
grid_lon, grid_lat = np.meshgrid(lon, lat)
try:
from scipy.interpolate import griddata
# linear interpolation only — leaves NaN (→ transparent) outside
# the data hull so the background does not extrapolate into
# unsurveyed ground
vals = griddata(
pts[:, :2], pts[:, 2], (grid_lon, grid_lat), method="linear"
)
except Exception: # noqa: BLE001 — SciPy missing / degenerate points
return None
if not np.isfinite(vals).any():
return None
return {
"lons": lon,
"lats": lat,
"values": vals,
"cmap": cmap,
"alpha": float(alpha),
"label": label,
"cbar_side": "bottom",
}
[docs]
def plot_phase_tensor_map(
sites: Any,
*,
period: float = 10.0,
# ── ellipse style (honours PhaseTensorEllipseStyle) ───────────────
scale=_UNSET,
ellipse_scale: float = 1.0,
normalise_by=_UNSET,
s1_ref=_UNSET,
min_aspect=_UNSET,
c_by=_UNSET,
cmap=_UNSET,
clim: tuple[float, float] | None = None,
clim_pct=_UNSET,
symmetric_clim=_UNSET,
alpha=_UNSET,
edgecolor=_UNSET,
linewidth=_UNSET,
skew_threshold=_UNSET,
mark_3d=_UNSET,
lw_3d_factor=_UNSET,
ref_ellipse=_UNSET,
# ── tipper arrows ─────────────────────────────────────────────────
show_tipper: bool = True,
tipper_convention: str = "parkinson",
tipper_component: str = "real",
tipper_scale: float | None = None,
tipper_color: str = "k",
tipper_lw: float = 1.4,
# ── optional background grid ──────────────────────────────────────
bg_grid: dict[str, Any] | None = None,
topography: Any = None,
topography_cmap: str = "terrain",
topography_alpha: float = 0.55,
# ── map decoration ────────────────────────────────────────────────
station_labels: bool = True,
station_marker: str = "v",
station_ms: float = 5.0,
station_color: str = "k",
label_fontsize: float = 7.0,
title: str = "",
colorbar_label: str | None = None,
show_colorbar: bool = True,
# ── coordinate override (station → (lat, lon)) ────────────────────
coords: dict[str, tuple[float, float]] | None = None,
# ── layout ────────────────────────────────────────────────────────
figsize: tuple[float, float] = (9.0, 7.0),
# ── core ──────────────────────────────────────────────────────────
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
ax: plt.Axes | None = None,
) -> plt.Axes:
"""Publication-quality phase-tensor map at a single period.
Renders each station as a phase-tensor ellipse positioned at its
geographic coordinates. The ellipse shape encodes the phase-tensor
principal values (φ_max, φ_min) and orientation (θ); fill colour
encodes an additional scalar (default: skewness β). Induction arrows
can be overlaid when tipper data are available, and an optional
background field (gravity, resistivity, …) may be drawn behind the
ellipses.
Parameters
----------
sites : any
EDI paths, objects, or collection accepted by
:func:`~pycsamt.emtools._core.ensure_sites`.
period : float
Target period (s). The nearest available period in the data is
used for each station independently.
scale : float or _UNSET
Maximum ellipse semi-axis in geographic units (degrees).
``_UNSET`` → auto-derived from the median inter-station spacing.
ellipse_scale : float, default 1.0
Multiplier applied on top of *scale* — a quick "make the
ellipses bigger/smaller" knob (``1.5`` = 50 % larger).
normalise_by : ``"cell"`` | ``"unity"`` | ``"abs"``
Size normalisation strategy, see
:class:`~pycsamt.api.style.PhaseTensorEllipseStyle`.
s1_ref : float or None
Manual reference s1 used with ``normalise_by="cell"``/``"unity"``.
min_aspect : float or _UNSET
Minimum height/width ratio enforced on every ellipse; see
:func:`plot_phase_tensor_psection`. Default:
``PYCSAMT_STYLE.pt_ellipse.min_aspect``.
c_by : str
Scalar quantity to map to fill colour. Recognised values:
``"skew"``, ``"beta"``, ``"|skew|"``, ``"theta"``,
``"ellipt"``, ``"phi_mean"``, ``"phi_max"``, ``"phi_min"``,
``"s1"``, ``"s2"``, ``"alpha"``. Defaults to
``PYCSAMT_STYLE.pt_ellipse.c_by``.
cmap : str or None
Matplotlib colormap. Auto-selected from *c_by* when ``None``
(same logic as :class:`~pycsamt.api.style.PhaseTensorEllipseStyle`).
clim : (vmin, vmax) or None
Explicit colour limits. ``None`` → derived from *clim_pct*.
clim_pct : (lo, hi)
Percentile window for automatic colour limits.
symmetric_clim : bool
Force vmin = −vmax (useful for diverging quantities).
alpha : float
Ellipse fill opacity.
edgecolor : str
Ellipse border colour.
linewidth : float
Normal ellipse border width (pts).
skew_threshold : float or None
|β| above which a cell is flagged as 3-D.
mark_3d : bool
Draw thicker borders on 3-D flagged ellipses.
lw_3d_factor : float
Line-width multiplier for 3-D flagged borders.
ref_ellipse : bool
Draw a reference circle in the lower-left corner as a scale bar.
show_tipper : bool
Overlay induction arrows when tipper data are found.
tipper_convention : ``"parkinson"`` | ``"wiese"``
Parkinson: arrow toward anomaly (negated real T);
Wiese: arrow along real T.
tipper_component : ``"real"`` | ``"imag"`` | ``"both"``
Which tipper component to draw.
tipper_scale : float or None
Arrow length = tipper_scale × *scale*. Auto-derived when ``None``.
tipper_color : str
Arrow colour (real component). Imaginary component uses
``mcolors.to_rgba(tipper_color, 0.55)``.
tipper_lw : float
Arrow line-width.
bg_grid : dict or None
Optional background field drawn behind ellipses::
bg_grid = dict(
lons = 1D or 2D longitude array,
lats = 1D or 2D latitude array,
values = 2D array (shape matches meshgrid of lons × lats),
cmap = "RdYlGn", # colormap
clim = (vmin, vmax), # or None → auto
alpha = 0.55, # opacity
label = "Gravity (gu)",
)
topography : bool | str | dict | None, default None
Convenience elevation background (ignored when *bg_grid* is
given). ``True`` interpolates the station elevations from the
EDI headers; a CSV path supplies either ``lon,lat,elev`` rows
(a scattered DEM) or ``station,elev`` rows; a dict is used
directly as *bg_grid*.
topography_cmap : str, default ``"terrain"``
Colormap for the *topography* background.
topography_alpha : float, default 0.55
Opacity for the *topography* background.
station_labels : bool
Annotate each station position with its name.
station_marker : str
Marker style for station positions (default ``"v"``).
station_ms : float
Marker size.
station_color : str
Marker and label colour.
label_fontsize : float
Font size for station labels.
title : str
Axes title.
colorbar_label : str or None
Override the automatic colorbar label derived from *c_by*.
show_colorbar : bool, default True
Draw the fill-colour colorbar. Set ``False`` when tiling several
maps that share one figure-level colorbar (see
:func:`plot_phase_tensor_map_grid`).
coords : dict[str, (lat, lon)] or None
Explicit station coordinates. When ``None`` the function reads
``ed.coords`` (a ``(lat, lon, elev)`` tuple) from each Site object.
figsize : (float, float)
Figure size in inches.
recursive, on_dup, strict, verbose
Passed to :func:`~pycsamt.emtools._core.ensure_sites`.
ax : matplotlib.axes.Axes or None
Pre-existing axes to draw into.
Returns
-------
matplotlib.axes.Axes
Examples
--------
Default skew-coloured map:
>>> from pycsamt.emtools import plot_phase_tensor_map
>>> ax = plot_phase_tensor_map(sites, period=10.0)
Ellipticity coloured, with tipper arrows:
>>> ax = plot_phase_tensor_map(
... sites,
... period=10.0,
... c_by="ellipt",
... cmap="viridis",
... show_tipper=True,
... tipper_convention="parkinson",
... )
With a gravity background:
>>> ax = plot_phase_tensor_map(
... sites,
... period=10.0,
... bg_grid=dict(
... lons=g_lon,
... lats=g_lat,
... values=g_bouguer,
... cmap="RdYlGn",
... alpha=0.45,
... label="Gravity (gu)",
... ),
... )
"""
import matplotlib.colors as _mc
from ..api.plot import add_colorbar
# ── 1. resolve style ──────────────────────────────────────────────────
_es = PYCSAMT_STYLE.pt_ellipse
def _sv(v, attr):
return getattr(_es, attr) if v is _UNSET else v
normalise_by = _sv(normalise_by, "normalise_by")
s1_ref_ = _sv(s1_ref, "s1_ref")
min_aspect_ = _sv(min_aspect, "min_aspect")
c_by = _sv(c_by, "c_by")
cmap_ = _sv(cmap, "cmap")
clim_pct = _sv(clim_pct, "clim_pct")
symmetric_clim = _sv(symmetric_clim, "symmetric_clim")
alpha_ = _sv(alpha, "alpha")
edgecolor_ = _sv(edgecolor, "edgecolor")
linewidth_ = _sv(linewidth, "linewidth")
skew_threshold_ = _sv(skew_threshold, "skew_threshold")
mark_3d_ = _sv(mark_3d, "mark_3d")
lw_3d_factor_ = _sv(lw_3d_factor, "lw_3d_factor")
ref_ellipse_ = _sv(ref_ellipse, "show_ref")
# auto-resolve cmap from c_by when None
if cmap_ is None:
cmap_ = _es.copy(c_by=c_by).resolve_cmap()
# honour natural symmetric_clim for diverging quantities
sym = symmetric_clim and _es.copy(c_by=c_by).resolve_symmetric_clim()
# ── 2. load data ──────────────────────────────────────────────────────
S = ensure_sites(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
df = build_phase_tensor_table(
S,
recursive=False,
on_dup=on_dup,
strict=False,
verbose=verbose,
)
if ax is None:
fig, ax = plt.subplots(figsize=figsize, constrained_layout=True)
if df.empty:
ax.text(
0.5,
0.5,
"no phase tensor data",
ha="center",
va="center",
transform=ax.transAxes,
)
return ax
# ── 3. station coordinates ────────────────────────────────────────────
if coords is None:
coords = {}
for i, ed in enumerate(_iter_items(S)):
st = _name(ed, i)
# preferred: ed.coords → (lat, lon, elev)
# A missing EDI ``>HEAD`` LAT/LONG (e.g. KAP03, which only
# carries REFLAT/REFLONG in >=DEFINEMEAS) surfaces as
# ``.coords`` returning NaN, not None; the "is not None"
# checks alone let NaN through into `coords`, and NaN then
# crashes ax.set_xlim() below instead of reaching the
# intended "no geographic coordinates" message.
c = getattr(ed, "coords", None)
if (
c is not None
and len(c) >= 2
and c[0] is not None
and c[1] is not None
and np.isfinite(c[0])
and np.isfinite(c[1])
):
coords[st] = (float(c[0]), float(c[1]))
continue
# fallback: ed.meta
meta = getattr(ed, "meta", {}) or {}
lat = meta.get("lat") or meta.get("latitude")
lon = meta.get("lon") or meta.get("longitude") or meta.get("long")
if (
lat is not None
and lon is not None
and np.isfinite(float(lat))
and np.isfinite(float(lon))
):
coords[st] = (float(lat), float(lon))
if not coords:
ax.text(
0.5,
0.5,
"no geographic coordinates in EDI data\n"
"pass coords={station: (lat, lon)} to override",
ha="center",
va="center",
transform=ax.transAxes,
fontsize=9,
style="italic",
color="0.45",
)
return ax
# ── 4. select nearest period per station ─────────────────────────────
rows: list[Any] = []
for st, sdf in df.groupby("station"):
if st not in coords:
continue
p = sdf["period"].to_numpy()
idx = int(np.nanargmin(np.abs(p - period)))
rows.append(sdf.iloc[idx])
if not rows:
ax.text(
0.5,
0.5,
f"no data near T = {period:g} s",
ha="center",
va="center",
transform=ax.transAxes,
)
return ax
st_names = [r["station"] for r in rows]
lats = np.array([coords[s][0] for s in st_names], float)
lons = np.array([coords[s][1] for s in st_names], float)
s1s = np.array([float(r["s1"]) for r in rows], float)
s2s = np.array([float(r["s2"]) for r in rows], float)
ths = np.array([float(r["theta"]) for r in rows], float)
np.array(
[float(r.get("skew", r.get("beta", np.nan))) for r in rows], float
)
# ── 5. auto scale from grid spacing ──────────────────────────────────
if scale is _UNSET:
if len(lons) > 1:
# Use unique grid spacings (cross-profile > within-profile)
u_lons = np.unique(np.round(lons, 3))
u_lats = np.unique(np.round(lats, 3))
d_lon = (
float(np.median(np.diff(u_lons)))
if len(u_lons) > 1
else float(np.ptp(lons) or 0.01)
)
d_lat = (
float(np.median(np.diff(u_lats)))
if len(u_lats) > 1
else float(np.ptp(lats) or 0.01)
)
scale_ = 0.55 * max(abs(d_lon) + 1e-9, abs(d_lat) + 1e-9)
else:
scale_ = 0.01
else:
scale_ = float(scale)
scale_ = scale_ * max(float(ellipse_scale), 1e-6)
# ── 6. ellipse size normalisation ────────────────────────────────────
# Clip extreme s1 outliers (robust 95th pct) before normalising
s1s_clip = np.clip(s1s, 0.0, float(np.nanpercentile(s1s, 95)))
if normalise_by == "unity":
_s1ref = (
s1_ref_ if s1_ref_ is not None else float(np.tan(np.radians(45.0)))
)
elif normalise_by == "cell":
_s1ref = (
s1_ref_
if s1_ref_ is not None
else float(np.nanpercentile(s1s_clip, 90))
)
else: # "abs"
_s1ref = s1_ref_ if s1_ref_ is not None else 1.0
_s1ref = max(_s1ref, 1e-9)
# ── 7. colour mapping ─────────────────────────────────────────────────
c_vals = np.array([float(r.get(c_by, np.nan)) for r in rows], float)
if clim is not None:
v0, v1 = float(clim[0]), float(clim[1])
else:
lo_pct, hi_pct = float(clim_pct[0]), float(clim_pct[1])
v0 = float(np.nanpercentile(c_vals, lo_pct))
v1 = float(np.nanpercentile(c_vals, hi_pct))
if sym and v0 != v1:
vmax = max(abs(v0), abs(v1))
v0, v1 = -vmax, vmax
v0 = v0 if v0 != v1 else v0 - 1.0
cm_obj = plt.get_cmap(cmap_)
norm_c = plt.Normalize(vmin=v0, vmax=v1)
# ── 8. optional background grid ───────────────────────────────────────
if bg_grid is None and topography is not None and topography is not False:
bg_grid = _topography_bg_grid(
topography,
S,
coords,
cmap=topography_cmap,
alpha=topography_alpha,
)
bg_mappable = None
if bg_grid is not None:
g_lon = np.asarray(bg_grid.get("lons", []), float)
g_lat = np.asarray(bg_grid.get("lats", []), float)
g_val = np.asarray(bg_grid.get("values", []), float)
g_cm = bg_grid.get("cmap", "RdYlGn")
g_cl = bg_grid.get("clim", None)
g_al = float(bg_grid.get("alpha", 0.5))
g_lbl = bg_grid.get("label", "")
g_side = bg_grid.get("cbar_side", "bottom")
if g_val.size > 0:
if g_val.ndim == 1:
g_val = g_val.reshape(len(g_lat), len(g_lon))
gv0 = (
float(np.nanpercentile(g_val, 5))
if g_cl is None
else float(g_cl[0])
)
gv1 = (
float(np.nanpercentile(g_val, 95))
if g_cl is None
else float(g_cl[1])
)
bg_mappable = ax.pcolormesh(
g_lon,
g_lat,
g_val,
cmap=g_cm,
vmin=gv0,
vmax=gv1,
alpha=g_al,
shading="auto",
zorder=0,
)
if g_lbl and g_side == "bottom":
cbg = ax.figure.colorbar(
bg_mappable,
ax=ax,
location="bottom",
fraction=0.045,
shrink=0.7,
pad=0.12,
aspect=40,
)
cbg.set_label(g_lbl, fontsize=8)
cbg.ax.tick_params(labelsize=7)
elif g_lbl:
add_colorbar(
bg_mappable,
ax,
label=g_lbl,
side=g_side,
size="3.5%",
pad=0.08,
)
# ── 9. draw ellipses ──────────────────────────────────────────────────
for row, lat, lon, s1, s2, th in zip(rows, lats, lons, s1s_clip, s2s, ths):
# clamp to 2×scale_ so outlier sites don't swamp the map
w = min(2.0 * scale_ * (s1 / _s1ref), 2.0 * scale_)
h = min(2.0 * scale_ * (s2 / _s1ref), 2.0 * scale_)
if min_aspect_ > 0:
h = max(h, min_aspect_ * w)
cv = float(row.get(c_by, np.nan))
fc = cm_obj(norm_c(cv)) if np.isfinite(cv) else (0.8, 0.8, 0.8, 0.8)
# 3-D flagging
beta = float(row.get("skew", row.get("beta", 0.0)))
is_3d = (
skew_threshold_ is not None
and mark_3d_
and abs(beta) > skew_threshold_
)
lw = linewidth_ * lw_3d_factor_ if is_3d else linewidth_
ec = edgecolor_
e = Ellipse(
xy=(lon, lat),
width=w,
height=h,
angle=th,
facecolor=fc,
edgecolor=ec,
linewidth=lw,
alpha=alpha_,
zorder=2,
)
ax.add_patch(e)
# ── 10. tipper arrows ─────────────────────────────────────────────────
if show_tipper:
_tip_data: dict[str, tuple[float, float]] = {}
for i, ed in enumerate(_iter_items(S)):
st = _name(ed, i)
if st not in coords:
continue
_T, t, fr = _get_t_block(ed)
if t is None or fr is None:
continue
per_arr = 1.0 / np.where(fr == 0, np.nan, fr)
j = int(np.nanargmin(np.abs(per_arr - period)))
tx, ty = t[j, 0], t[j, 1]
if tipper_component == "real":
u, v = float(np.real(ty)), float(np.real(tx)) # E, N
elif tipper_component == "imag":
u, v = float(np.imag(ty)), float(np.imag(tx))
else:
u, v = float(np.real(ty)), float(np.real(tx))
if tipper_convention == "parkinson":
u, v = -u, -v
_tip_data[st] = (u, v)
if _tip_data:
# normalise arrow lengths
magnitudes = np.array(
[np.hypot(u, v) for u, v in _tip_data.values()], float
)
mag_max = float(np.nanpercentile(magnitudes, 95)) + 1e-12
t_scale = (
tipper_scale if tipper_scale is not None else scale_ * 1.6
)
t_lons, t_lats, t_u, t_v = [], [], [], []
for st_t, (u, v) in _tip_data.items():
lat_t, lon_t = coords[st_t]
t_lons.append(lon_t)
t_lats.append(lat_t)
t_u.append(u / mag_max * t_scale)
t_v.append(v / mag_max * t_scale)
ax.quiver(
np.array(t_lons),
np.array(t_lats),
np.array(t_u),
np.array(t_v),
color=tipper_color,
scale=1.0,
scale_units="xy",
angles="xy",
width=0.003,
headwidth=4,
headlength=5,
linewidth=tipper_lw,
zorder=4,
)
# imaginary component if "both"
if tipper_component == "both":
_tip_im: dict[str, tuple[float, float]] = {}
for i, ed in enumerate(_iter_items(S)):
st = _name(ed, i)
if st not in coords:
continue
_T, t, fr = _get_t_block(ed)
if t is None or fr is None:
continue
per_arr = 1.0 / np.where(fr == 0, np.nan, fr)
j = int(np.nanargmin(np.abs(per_arr - period)))
tx, ty = t[j, 0], t[j, 1]
u_i, v_i = float(np.imag(ty)), float(np.imag(tx))
if tipper_convention == "parkinson":
u_i, v_i = -u_i, -v_i
_tip_im[st] = (u_i, v_i)
if _tip_im:
im_lons, im_lats, im_u, im_v = [], [], [], []
for st_t, (u, v) in _tip_im.items():
lat_t, lon_t = coords[st_t]
im_lons.append(lon_t)
im_lats.append(lat_t)
im_u.append(u / mag_max * t_scale)
im_v.append(v / mag_max * t_scale)
ax.quiver(
np.array(im_lons),
np.array(im_lats),
np.array(im_u),
np.array(im_v),
color=_mc.to_rgba(tipper_color, 0.55),
scale=1.0,
scale_units="xy",
angles="xy",
width=0.002,
headwidth=3,
headlength=4,
linestyle="dashed",
zorder=3,
)
# ── 11. station markers + labels ──────────────────────────────────────
ax.scatter(
lons,
lats,
marker=station_marker,
c=station_color,
s=station_ms**2,
zorder=5,
linewidths=0,
)
if station_labels:
offset_y = scale_ * 0.8
for st, lat, lon in zip(st_names, lats, lons):
ax.annotate(
st,
(lon, lat),
xytext=(0, offset_y * 111000 * 0.00001),
textcoords="offset points",
fontsize=label_fontsize,
color=station_color,
ha="center",
va="bottom",
zorder=6,
)
# ── 12. reference ellipse ─────────────────────────────────────────────
if ref_ellipse_:
pad = scale_ * 0.7
ref_lon = float(lons.min()) - pad
ref_lat = float(lats.min()) - pad
ref_w = 2.0 * scale_
ref_h = 2.0 * scale_
ref_e = Ellipse(
(ref_lon, ref_lat),
width=ref_w,
height=ref_h,
facecolor="none",
edgecolor=_es.ref_edgecolor,
linewidth=_es.ref_lw,
zorder=5,
)
ax.add_patch(ref_e)
ax.text(
ref_lon,
ref_lat - scale_ * 1.3,
"ref",
ha="center",
va="top",
fontsize=_es.ref_fontsize,
color=_es.ref_edgecolor,
)
# ── 13. colorbar ──────────────────────────────────────────────────────
if show_colorbar:
sm = plt.cm.ScalarMappable(cmap=cm_obj, norm=norm_c)
sm.set_array([])
cb_label = colorbar_label if colorbar_label is not None else c_by
add_colorbar(
sm, ax, label=cb_label, side="right", size="3.5%", pad=0.06
)
# ── 14. axes decoration ───────────────────────────────────────────────
from matplotlib.ticker import ScalarFormatter
pad_xy = scale_ * 2.2
ax.set_xlim(lons.min() - pad_xy, lons.max() + pad_xy)
ax.set_ylim(lats.min() - pad_xy, lats.max() + pad_xy)
# geographic aspect ratio: 1° lon = cos(lat) × 1° lat in physical distance
lat_mid = float(np.mean(lats))
ax.set_aspect(1.0 / max(np.cos(np.radians(lat_mid)), 1e-6))
ax.set_xlabel("Longitude", fontsize=9)
ax.set_ylabel("Latitude", fontsize=9)
ax.tick_params(labelsize=8)
# suppress matplotlib offset notation (e.g. "+1.1912e2")
for _axis in (ax.xaxis, ax.yaxis):
_fmt = ScalarFormatter(useOffset=False)
_fmt.set_scientific(False)
_axis.set_major_formatter(_fmt)
ax.grid(True, alpha=0.18, linewidth=0.5)
period_txt = (
f"{period * 1e3:.4g} ms" if period < 1.0 else f"{period:.4g} s"
)
ax.set_title(
title or f"Phase Tensor Map for {period_txt}",
fontsize=10,
pad=8,
fontweight="bold",
)
return ax
def _grid_station_lonlat(sites, coords):
"""Return ``{station: (lat, lon)}`` for a phase-tensor map grid.
Mirrors the coordinate resolution of :func:`plot_phase_tensor_map`
(``ed.coords`` then ``ed.meta``) but only for the shared ellipse-size
estimate; the per-panel call resolves coordinates itself.
"""
if coords is not None:
return {k: (float(v[0]), float(v[1])) for k, v in coords.items()}
out: dict[str, tuple[float, float]] = {}
for i, ed in enumerate(_iter_items(sites)):
st = _name(ed, i)
c = getattr(ed, "coords", None)
if (
c is not None
and len(c) >= 2
and c[0] is not None
and c[1] is not None
and np.isfinite(c[0])
and np.isfinite(c[1])
):
out[st] = (float(c[0]), float(c[1]))
continue
meta = getattr(ed, "meta", {}) or {}
lat = meta.get("lat") or meta.get("latitude")
lon = meta.get("lon") or meta.get("longitude") or meta.get("long")
if lat is not None and lon is not None:
try:
if np.isfinite(float(lat)) and np.isfinite(float(lon)):
out[st] = (float(lat), float(lon))
except (TypeError, ValueError):
pass
return out
_MAP_GRID_NCOLS = {1: 1, 2: 2, 3: 3, 4: 2, 5: 3, 6: 3}
[docs]
def plot_phase_tensor_map_grid(
sites: Any,
*,
frequencies: list[float] | None = None,
periods: list[float] | None = None,
n_cols: int | None = None,
panel_labels: bool | list[str] = True,
panel_titles: list[str] | None = None,
tidy_labels: bool = True,
label_by: str = "frequency",
share_scale: bool = True,
share_color: bool = True,
shared_colorbar: bool = True,
station_labels: str | bool = "first",
ref_ellipse: str | bool = "first",
bg_grid: dict[str, Any] | None = None,
topography: Any = None,
topography_cmap: str = "terrain",
topography_alpha: float = 0.55,
suptitle: str = "",
figsize: tuple[float, float] | None = None,
panel_size: tuple[float, float] = (4.6, 4.0),
# ── shared ellipse / colour style ────────────────────────────────
scale=_UNSET,
ellipse_scale: float = 1.0,
s1_ref=_UNSET,
normalise_by=_UNSET,
c_by=_UNSET,
abs_skew: bool = False,
cmap=_UNSET,
clim: tuple[float, float] | None = None,
clim_pct=_UNSET,
symmetric_clim=_UNSET,
skew_threshold=_UNSET,
colorbar_label: str | None = None,
# ── core ────────────────────────────────────────────────────────
coords: dict[str, tuple[float, float]] | None = None,
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
axes: Any | None = None,
**map_kwargs: Any,
) -> plt.Figure:
"""Multi-frequency phase-tensor map — one panel per frequency.
Tiles :func:`plot_phase_tensor_map` across a list of frequencies (or
periods) in a subplot grid, with a **single shared colorbar** and,
by default, a shared ellipse-size reference and colour scale so the
panels are directly comparable. This reproduces the classic
"phase-tensor ellipses filled with skew, induction arrows in
Parkinson convention, at 30 / 3 / 0.3 / 0.03 Hz" figure.
Any number of panels is accepted (1, 2, 3, 4, …); the layout adapts.
Parameters
----------
sites : any
EDI paths, objects, or a collection accepted by
:func:`~pycsamt.emtools._core.ensure_sites`.
frequencies : list of float, optional
Panel frequencies in **hertz**, e.g. ``[30, 3, 0.3, 0.03]``.
Exactly one of *frequencies* / *periods* must be given.
periods : list of float, optional
Panel periods in **seconds**. Alternative to *frequencies*.
n_cols : int, optional
Columns in the panel grid. Default: 1→1, 2→2, 3→3, 4→2×2,
then 3 per row.
panel_labels : bool or list of str, default True
``True`` stamps ``"(a)"``, ``"(b)"`` … in each panel corner; a
list sets the text explicitly; ``False`` disables it.
panel_titles : list of str, optional
Per-panel titles. Default: ``"30 Hz"`` … (or ``"0.033 s"`` …
when ``label_by="period"``).
tidy_labels : bool, default True
Drop the x-axis label/ticks on non-bottom panels and the y-axis
label/ticks on non-left panels (every panel shares the same
extent), which tightens the grid.
label_by : {"frequency", "period"}, default "frequency"
Units for the default panel titles.
share_scale : bool, default True
Use one ellipse-size reference for every panel (ellipse areas
become comparable across frequencies).
share_color : bool, default True
Use one colour scale for every panel. Ignored when *clim* is
given (that is always shared).
shared_colorbar : bool, default True
Draw a single figure-level colorbar instead of one per panel.
station_labels : {"first", "all", "none"} or bool, default "first"
Which panels annotate station names.
ref_ellipse : {"first", "all", "none"} or bool, default "first"
Which panels draw the reference-circle scale bar.
bg_grid : dict, optional
Explicit background field (gravity, resistivity, …) drawn behind
every panel; see :func:`plot_phase_tensor_map`. Takes priority
over *topography*.
topography : bool | str | dict | None, default None
Convenience elevation background shared by all panels. ``True``
interpolates the station elevations from the EDI headers; a CSV
path supplies ``lon,lat,elev`` or ``station,elev`` rows; a dict
is used as *bg_grid*. A single horizontal "Elevation (m)"
colorbar is added.
topography_cmap : str, default ``"terrain"``
topography_alpha : float, default 0.55
Colormap / opacity for the *topography* background.
suptitle : str, optional
Figure-level title.
figsize : tuple of float, optional
Overrides ``panel_size``-derived sizing.
panel_size : tuple of float, default (4.6, 4.0)
Per-panel ``(width, height)`` in inches.
ellipse_scale : float, default 1.0
Multiplier on the ellipse size for every panel — a quick
"bigger/smaller ellipses" knob (``1.5`` = 50 % larger).
abs_skew : bool, default False
Colour by the **absolute** skew ``|β|`` (sequential) instead of
signed ``β`` (diverging). Shorthand for ``c_by="|beta|"``.
scale, s1_ref, normalise_by, c_by, cmap, clim, clim_pct, \
symmetric_clim, skew_threshold, colorbar_label :
Ellipse-style controls forwarded to every panel; see
:func:`plot_phase_tensor_map`. When *share_scale* /
*share_color* is on and the value is left unset, a single shared
value is derived from the pooled data of every panel.
coords : dict[str, (lat, lon)], optional
Explicit station coordinates, forwarded to every panel.
recursive, on_dup, strict, verbose
Forwarded to :func:`~pycsamt.emtools._core.ensure_sites`.
axes : array of matplotlib.axes.Axes, optional
Pre-existing axes grid (flattened length ≥ number of panels).
**map_kwargs
Any other keyword (``show_tipper``, ``tipper_convention``,
``tipper_component``, ``tipper_scale``, ``tipper_color``,
``tipper_lw``, ``station_marker``, ``alpha``, ``edgecolor`` …)
is passed straight through to every
:func:`plot_phase_tensor_map` call.
Returns
-------
matplotlib.figure.Figure
Raises
------
ValueError
If neither or both of *frequencies* / *periods* are given, or the
list is empty.
Examples
--------
Broken Hill Fig. 3 style — signed skew β on the MTPy-style
``pt_skew`` colormap, Parkinson arrows:
>>> from pycsamt.emtools import plot_phase_tensor_map_grid
>>> fig = plot_phase_tensor_map_grid(
... "data/MT/broken-hill/edis",
... frequencies=[30, 3, 0.3, 0.03],
... c_by="skew", # signed β; auto "pt_skew" cmap
... tipper_convention="parkinson",
... suptitle="Broken Hill — phase tensor + induction arrows",
... )
Two panels only:
>>> fig = plot_phase_tensor_map_grid(sites, frequencies=[10, 0.1])
See Also
--------
plot_phase_tensor_map : Single-frequency phase-tensor map.
plot_phase_tensor_strip_grid : Ellipse strips tiled by profile.
"""
if (frequencies is None) == (periods is None):
raise ValueError(
"pass exactly one of frequencies= (Hz) or periods= (s)."
)
if periods is not None:
target_p = [float(p) for p in periods]
freqs = [1.0 / p if p else np.inf for p in target_p]
else:
freqs = [float(f) for f in frequencies]
target_p = [1.0 / f if f else np.inf for f in freqs]
n = len(target_p)
if n == 0:
raise ValueError("frequencies/periods list is empty.")
_es = PYCSAMT_STYLE.pt_ellipse
def _sv(v, attr):
return getattr(_es, attr) if v is _UNSET else v
c_by_ = _sv(c_by, "c_by")
if abs_skew and c_by_ in ("skew", "beta"):
c_by_ = "|beta|"
cmap_ = cmap if cmap is not _UNSET else None
if cmap_ is None:
cmap_ = _es.copy(c_by=c_by_).resolve_cmap()
clim_pct_ = _sv(clim_pct, "clim_pct")
symmetric_clim_ = symmetric_clim
if symmetric_clim_ is _UNSET:
symmetric_clim_ = _es.copy(c_by=c_by_).resolve_symmetric_clim()
skew_threshold_ = _sv(skew_threshold, "skew_threshold")
normalise_by_ = _sv(normalise_by, "normalise_by")
# ── load once ───────────────────────────────────────────────────
S = ensure_sites(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
df = build_phase_tensor_table(
S, recursive=False, on_dup=on_dup, strict=False, verbose=verbose
)
# ── rows nearest to each target period, per station ─────────────
sub_rows: list[Any] = []
if not df.empty:
for _, sdf in df.groupby("station"):
p = sdf["period"].to_numpy(float)
for tp in target_p:
sub_rows.append(sdf.iloc[int(np.nanargmin(np.abs(p - tp)))])
df_sub = pd.DataFrame(sub_rows) if sub_rows else df
# ── shared colour limits ────────────────────────────────────────
if clim is not None:
clim_shared: tuple[float, float] = (
float(clim[0]),
float(clim[1]),
)
cbar_label = colorbar_label or _resolve_cvals(df_sub, c_by_)[1]
elif share_color and not df_sub.empty:
cvals, cbar_label = _resolve_cvals(df_sub, c_by_)
finite = cvals[np.isfinite(cvals)]
if finite.size:
lo = float(np.nanpercentile(finite, clim_pct_[0]))
hi = float(np.nanpercentile(finite, clim_pct_[1]))
if symmetric_clim_ and c_by_ in _SYMMETRIC_C:
vlim = max(abs(lo), abs(hi))
# keep at least ±skew_threshold so a near-1-D survey does
# not saturate the whole diverging map at one hue
if skew_threshold_ and c_by_ in ("skew", "beta"):
vlim = max(vlim, float(skew_threshold_))
lo, hi = -vlim, vlim
clim_shared = (lo, hi if hi != lo else lo + 1.0)
else:
clim_shared = (-1.0, 1.0)
if colorbar_label:
cbar_label = colorbar_label
else:
clim_shared = None # each panel picks its own
cbar_label = colorbar_label or c_by_
# ── shared geographic scale ────────────────────────────────────
# Only the geographic ellipse/arrow *scale* is shared: each panel
# still normalises ellipse sizes to its own phase distribution
# (``s1_ref`` left unset), otherwise the small high-frequency phase
# tensors would collapse next to the large low-frequency ones. Pass
# an explicit ``s1_ref`` to compare magnitudes across frequency.
scale_arg = scale
s1_ref_arg = s1_ref
if share_scale and scale is _UNSET:
ll = _grid_station_lonlat(S, coords)
if len(ll) > 1:
lats = np.array([v[0] for v in ll.values()], float)
lons = np.array([v[1] for v in ll.values()], float)
u_lon = np.unique(np.round(lons, 3))
u_lat = np.unique(np.round(lats, 3))
d_lon = (
float(np.median(np.diff(u_lon)))
if len(u_lon) > 1
else float(np.ptp(lons) or 0.01)
)
d_lat = (
float(np.median(np.diff(u_lat)))
if len(u_lat) > 1
else float(np.ptp(lats) or 0.01)
)
scale_arg = 0.68 * max(abs(d_lon) + 1e-9, abs(d_lat) + 1e-9)
# ── shared topography background (built once, drawn on every panel)
topo_bg: dict[str, Any] | None = None
if (
bg_grid is None
and topography is not None
and topography is not False
):
topo_bg = _topography_bg_grid(
topography,
S,
_grid_station_lonlat(S, coords),
cmap=topography_cmap,
alpha=topography_alpha,
)
panel_bg = bg_grid
if panel_bg is None and topo_bg is not None:
panel_bg = dict(topo_bg)
panel_bg["label"] = "" # one shared elevation bar instead
# ── layout ─────────────────────────────────────────────────────
if n_cols is None:
n_cols = _MAP_GRID_NCOLS.get(n, 3)
n_cols = max(1, min(int(n_cols), n))
n_rows = int(np.ceil(n / n_cols))
if axes is None:
# Each panel is drawn with a latitude-corrected aspect, so match
# the subplot slot to that box (width/height) — otherwise the
# fixed-aspect axes leave wide gutters between columns.
box_wh = 1.0
_ll = _grid_station_lonlat(S, coords)
if len(_ll) > 1:
_la = np.array([v[0] for v in _ll.values()], float)
_lo = np.array([v[1] for v in _ll.values()], float)
_dlo = float(np.ptp(_lo)) or 1.0
_dla = float(np.ptp(_la)) or 1.0
box_wh = float(
np.clip(
_dlo * np.cos(np.radians(_la.mean())) / _dla,
0.45,
2.2,
)
)
panel_h = float(panel_size[1])
panel_w = panel_h * box_wh + 0.85 # + tick labels / y-axis title
fs = figsize or (
panel_w * n_cols + (0.95 if shared_colorbar else 0.15),
panel_h * n_rows + (0.55 if suptitle else 0.2),
)
fig, axgrid = plt.subplots(
n_rows, n_cols, figsize=fs, squeeze=False,
constrained_layout=True,
)
try:
fig.set_constrained_layout_pads(
w_pad=0.015, h_pad=0.02, wspace=0.015, hspace=0.03
)
except Exception: # noqa: BLE001 — older matplotlib
pass
else:
axgrid = np.atleast_2d(axes)
fig = axgrid.flat[0].figure
flat = list(np.asarray(axgrid).ravel())
def _panel_on(spec: str | bool, k: int) -> bool:
if isinstance(spec, bool):
return spec
spec = str(spec).lower()
if spec in ("all", "every"):
return True
if spec in ("none", "off", "never"):
return False
return k == 0 # "first"
for k, tp in enumerate(target_p):
ax = flat[k]
if panel_titles is not None and k < len(panel_titles):
ptitle = panel_titles[k]
elif label_by == "period":
ptitle = (
f"{tp * 1e3:.4g} ms" if tp < 1.0 else f"{tp:.4g} s"
)
else:
fk = freqs[k]
ptitle = f"{fk:g} Hz" if np.isfinite(fk) else "0 Hz"
plot_phase_tensor_map(
S,
period=tp,
ax=ax,
show_colorbar=False,
c_by=c_by_,
cmap=cmap_,
clim=clim_shared,
clim_pct=clim_pct_,
symmetric_clim=symmetric_clim_,
skew_threshold=skew_threshold_,
scale=scale_arg,
ellipse_scale=ellipse_scale,
s1_ref=s1_ref_arg,
normalise_by=normalise_by_,
bg_grid=panel_bg,
coords=coords,
station_labels=_panel_on(station_labels, k),
ref_ellipse=_panel_on(ref_ellipse, k),
title=ptitle,
recursive=False,
on_dup=on_dup,
strict=strict,
verbose=verbose,
**map_kwargs,
)
ax.locator_params(axis="x", nbins=4)
ax.locator_params(axis="y", nbins=5)
if tidy_labels and axes is None:
row, col = divmod(k, n_cols)
last_row_start = (n_rows - 1) * n_cols
if k < last_row_start and (k + n_cols) < n:
ax.set_xlabel("")
ax.tick_params(labelbottom=False)
if col != 0:
ax.set_ylabel("")
ax.tick_params(labelleft=False)
if panel_labels is not False:
lbl = (
panel_labels[k]
if isinstance(panel_labels, (list, tuple))
and k < len(panel_labels)
else f"({chr(97 + k)})"
)
ax.text(
0.02,
0.97,
lbl,
transform=ax.transAxes,
ha="left",
va="top",
fontsize=11,
fontweight="bold",
bbox={
"boxstyle": "square,pad=0.15",
"fc": "white",
"ec": "none",
"alpha": 0.75,
},
)
for k in range(n, len(flat)):
flat[k].set_visible(False)
if shared_colorbar and clim_shared is not None:
sm = ScalarMappable(
cmap=plt.get_cmap(cmap_),
norm=Normalize(vmin=clim_shared[0], vmax=clim_shared[1]),
)
sm.set_array([])
cb = fig.colorbar(
sm,
ax=flat[:n],
fraction=0.018,
pad=0.006,
shrink=0.92,
aspect=32,
)
cb.set_label(cbar_label, fontsize=9)
cb.ax.tick_params(labelsize=8)
if shared_colorbar and topo_bg is not None:
tv = np.asarray(topo_bg["values"], float)
tv = tv[np.isfinite(tv)]
if tv.size:
sm2 = ScalarMappable(
cmap=plt.get_cmap(topo_bg.get("cmap", topography_cmap)),
norm=Normalize(vmin=float(tv.min()), vmax=float(tv.max())),
)
sm2.set_array([])
cb2 = fig.colorbar(
sm2,
ax=flat[:n],
orientation="horizontal",
fraction=0.03,
shrink=0.55,
pad=0.04,
aspect=45,
)
cb2.set_label(topo_bg.get("label") or "Elevation (m)", fontsize=9)
cb2.ax.tick_params(labelsize=8)
if suptitle:
fig.suptitle(suptitle, fontsize=12)
return fig
[docs]
def plot_phase_tensor_summary(
sites: Any,
*,
stations: list[str] | None = None,
period_range: tuple[float, float] | None = None,
scale=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.scale
c_by=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.c_by
cmap=_UNSET, # default: auto from c_by via resolve_cmap()
clim: tuple[float, float] | None = None,
skew_threshold=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.skew_threshold
ellipt_threshold: float = 0.2,
axes=None,
figsize: tuple[float, float] = (12.0, 9.0),
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
) -> plt.Figure:
"""
Three-panel phase-tensor summary figure.
Panels
------
**(a)** Ellipse pseudo-section — the full :func:`plot_phase_tensor_psection`
with all cell-size normalisations.
**(b)** Dimensionality classification grid — 1-D (white), 2-D (grey),
3-D (red), derived from *skew_threshold* and *ellipt_threshold*.
**(c)** Skew–ellipticity density hexbin — joint distribution over all
``(station, period)`` cells, with contours at 20 / 50 / 80 % of the
peak density.
Parameters
----------
sites : any
Passed directly to :func:`build_phase_tensor_table`.
stations : list of str or None
Restrict to a station subset (applied to all panels).
period_range : (T_min, T_max) or None
Period window in seconds (applied to all panels).
scale : float, default ``0.85``
Passed to panel (a) — fraction of each cell to fill.
c_by : str, default ``"skew"``
Colour quantity for panel (a); see :func:`plot_phase_tensor_psection`.
cmap : str, default ``"RdBu_r"``
Colourmap for panel (a).
clim : (vmin, vmax) or None
Explicit colour limits for panel (a).
skew_threshold : float, default ``3.0``
|β| threshold (°) for 1-D/2-D vs 3-D in panels (a) and (b).
ellipt_threshold : float, default ``0.2``
Ellipticity threshold for 1-D vs 2-D in panel (b).
figsize : (width, height), default ``(12.0, 9.0)``
recursive, on_dup, strict, verbose : see :func:`ensure_sites`.
Returns
-------
fig : :class:`~matplotlib.figure.Figure`
"""
# ── resolve visual style from PYCSAMT_STYLE.pt_ellipse ──────────────
_es = PYCSAMT_STYLE.pt_ellipse
if scale is _UNSET:
scale = _es.scale
if c_by is _UNSET:
c_by = _es.c_by
if cmap is _UNSET:
cmap = _es.copy(c_by=c_by).resolve_cmap()
if skew_threshold is _UNSET:
skew_threshold = _es.skew_threshold
df = build_phase_tensor_table(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
# Apply filters once; pass the already-filtered data to sub-functions
if not df.empty and stations is not None:
df = df[df["station"].isin(stations)]
if not df.empty and period_range is not None:
lo, hi = float(period_range[0]), float(period_range[1])
df = df[(df["period"] >= lo) & (df["period"] <= hi)]
# ── layout: 2 rows × 2 cols; row-0 spans both columns ────────────────
axes_given = _axes_list(axes, 3) if axes is not None else None
if axes_given is None:
import matplotlib.gridspec as gridspec
fig = plt.figure(figsize=figsize)
gs = gridspec.GridSpec(
2,
2,
figure=fig,
height_ratios=[1.5, 1.0],
hspace=0.38,
wspace=0.35,
)
ax_ell = fig.add_subplot(gs[0, :]) # (a) full-width ellipse section
ax_dim = fig.add_subplot(gs[1, 0]) # (b) dimensionality grid
ax_den = fig.add_subplot(gs[1, 1]) # (c) skew-ellipticity density
else:
ax_ell, ax_dim, ax_den = axes_given
fig = ax_ell.figure
# ── panel (a): ellipse pseudo-section ────────────────────────────────
if df.empty:
ax_ell.text(
0.5,
0.5,
"no phase tensor data",
ha="center",
va="center",
transform=ax_ell.transAxes,
)
else:
plot_phase_tensor_psection(
df, # pass pre-filtered DataFrame directly
scale=scale,
c_by=c_by,
cmap=cmap,
clim=clim,
skew_threshold=skew_threshold,
mark_3d=True,
ref_ellipse=True,
title="(a) Phase-tensor ellipse pseudo-section",
recursive=False,
ax=ax_ell,
)
# ── panel (b): dimensionality grid ────────────────────────────────────
if df.empty:
ax_dim.text(
0.5,
0.5,
"no data",
ha="center",
va="center",
transform=ax_dim.transAxes,
)
else:
_draw_dim_grid(
df,
ax_dim,
skew_threshold=skew_threshold,
ellipt_threshold=ellipt_threshold,
)
ax_dim.set_title("(b) Dimensionality (1-D/2-D/3-D)", fontsize=9)
# ── panel (c): skew–ellipticity density ───────────────────────────────
if df.empty:
ax_den.text(
0.5,
0.5,
"no data",
ha="center",
va="center",
transform=ax_den.transAxes,
)
else:
_draw_skew_ellipt_density(
df,
ax_den,
skew_threshold=skew_threshold,
ellipt_threshold=ellipt_threshold,
)
ax_den.set_title("(c) Skew–ellipticity distribution", fontsize=9)
return fig
def _draw_dim_grid(
df: pd.DataFrame,
ax: plt.Axes,
*,
skew_threshold: float = 3.0,
ellipt_threshold: float = 0.2,
) -> None:
"""Internal helper: draw the dimensionality heatmap onto *ax*."""
df = df.copy()
df["_logp"] = np.log10(df["period"].to_numpy(float))
a = np.abs(df["skew"].to_numpy(float))
e = np.abs(df["ellipt"].to_numpy(float))
lab = np.full(len(df), 2, dtype=int) # default: 3-D
lab[(a <= skew_threshold) & (e <= ellipt_threshold)] = 0 # 1-D
lab[(a <= skew_threshold) & (e > ellipt_threshold)] = 1 # 2-D
df["_dim"] = lab
st_list = list(dict.fromkeys(df["station"].tolist()))
{s: i for i, s in enumerate(st_list)}
piv = df.pivot_table(
index="_logp",
columns="station",
values="_dim",
aggfunc="median",
)
piv = piv.reindex(columns=st_list).sort_index()
# piv has shape (n_logp, n_stations); do NOT transpose so that
# imshow maps x → stations and y → log-periods.
Z = piv.to_numpy(dtype=float) # (n_logp, n_st)
cmap_d = mcolors.ListedColormap(["#f7f7f7", "#969696", "#d6604d"])
bounds = [-0.5, 0.5, 1.5, 2.5]
norm_d = mcolors.BoundaryNorm(bounds, cmap_d.N)
ax.imshow(
Z,
aspect="auto",
origin="lower",
interpolation="nearest",
cmap=cmap_d,
norm=norm_d,
)
ax.set_ylabel(LOG10_PERIOD_LABEL, fontsize=8)
n_st = len(st_list)
PYCSAMT_STATION_RENDERING.apply(
ax,
np.arange(n_st, dtype=float),
st_list,
preset="pseudosection",
xlim=(-0.5, n_st - 0.5),
)
n_logp = Z.shape[0] # shape[0] = n_logp now that .T is removed
y_ticks = np.linspace(0, n_logp - 1, num=min(6, n_logp))
y_vals = np.linspace(piv.index.min(), piv.index.max(), num=min(6, n_logp))
ax.set_yticks(y_ticks)
ax.set_yticklabels([f"{v:.1f}" for v in y_vals], fontsize=7)
if not ax.yaxis_inverted():
ax.invert_yaxis()
from matplotlib.patches import Patch
ax.legend(
handles=[
Patch(fc="#f7f7f7", ec="k", lw=0.5, label="1-D"),
Patch(fc="#969696", ec="k", lw=0.5, label="2-D"),
Patch(fc="#d6604d", ec="k", lw=0.5, label="3-D"),
],
fontsize=7,
framealpha=0.9,
loc="lower right",
title=f"|β|≤{skew_threshold:.0f}°,ε≤{ellipt_threshold:.2f}",
title_fontsize=6,
)
def _draw_skew_ellipt_density(
df: pd.DataFrame,
ax: plt.Axes,
*,
skew_threshold: float = 3.0,
ellipt_threshold: float = 0.2,
gridsize: int = 30,
) -> None:
"""Internal helper: draw |β|–ellipticity hexbin density onto *ax*."""
a = np.abs(df["skew"].to_numpy(float))
e = np.abs(df["ellipt"].to_numpy(float))
fin = np.isfinite(a) & np.isfinite(e)
if not fin.any():
ax.text(
0.5,
0.5,
"no data",
ha="center",
va="center",
transform=ax.transAxes,
)
return
ax.hexbin(
a[fin],
e[fin],
gridsize=gridsize,
mincnt=1,
cmap="YlOrRd",
linewidths=0.0,
)
# contour at 20 / 50 / 80 % of peak density
nx = ny = gridsize
H, xe, ye = np.histogram2d(a[fin], e[fin], bins=[nx, ny])
if H.max() > 0:
Xc = 0.5 * (xe[1:] + xe[:-1])
Yc = 0.5 * (ye[1:] + ye[:-1])
levs = [H.max() * f for f in (0.20, 0.50, 0.80)]
ax.contour(Xc, Yc, H.T, levels=levs, colors="k", linewidths=0.6)
# threshold lines
ax.axvline(
skew_threshold,
color="#b2182b",
lw=1.0,
ls="--",
label=f"|β|={skew_threshold:.0f}°",
)
ax.axhline(
ellipt_threshold,
color="#2166ac",
lw=1.0,
ls="--",
label=f"ε={ellipt_threshold:.2f}",
)
ax.set_xlabel("|β| (°)", fontsize=8)
ax.set_ylabel("Ellipticity", fontsize=8)
ax.legend(fontsize=7, framealpha=0.9)
[docs]
def phase_tensor_legend(
*,
size: float = 1.0,
ellipt: float = 0.45,
theta_deg: float = 20.0,
figsize: tuple[float, float] = (2.8, 2.8),
ax: plt.Axes | None = None,
) -> plt.Axes:
r"""Draw a labeled reference phase-tensor ellipse.
Explains the ellipse convention shared by every phase-tensor plot in
this module -- :func:`plot_phase_tensor_psection`,
:func:`plot_phase_tensor_map`, :func:`plot_phase_tensor_strip`, and
friends -- by drawing one annotated example: a major axis
:math:`\phi_{\max}`, a minor axis :math:`\phi_{\min}`, and an
orientation angle ``theta`` measured counterclockwise from the
horizontal (dashed reference line), matching the ``width``,
``height``, ``angle`` convention Matplotlib's own ``Ellipse`` patch
uses.
Parameters
----------
size : float, default 1.0
Semi-major axis length (:math:`\phi_{\max}`) of the example
ellipse, in axis data units.
ellipt : float, default 0.45
Example ellipticity used only to make the minor axis visibly
shorter than the major axis; :math:`\phi_{\min} = \text{size}
\times (1-\text{ellipt})`. Purely illustrative -- it does not
come from real data.
theta_deg : float, default 20.0
Example orientation angle, in degrees, used to draw the
rotated ellipse and the ``theta`` arc.
figsize : (float, float), default (2.8, 2.8)
Figure size when *ax* is not supplied.
ax : matplotlib.axes.Axes or None
Axes to draw on; created if *None*.
Returns
-------
matplotlib.axes.Axes
"""
if ax is None:
_, ax = plt.subplots(figsize=figsize)
phi_max = float(size)
phi_min = phi_max * max(0.05, 1.0 - float(ellipt))
theta = np.radians(float(theta_deg))
ax.add_patch(
Ellipse(
(0.0, 0.0),
width=2.0 * phi_max,
height=2.0 * phi_min,
angle=float(theta_deg),
fill=False,
lw=1.4,
edgecolor="crimson",
)
)
major = np.array([np.cos(theta), np.sin(theta)]) * phi_max
minor = np.array([-np.sin(theta), np.cos(theta)]) * phi_min
ax.plot([-major[0], major[0]], [-major[1], major[1]], "-", lw=1.0, color="0.15")
ax.plot([-minor[0], minor[0]], [-minor[1], minor[1]], "-", lw=1.0, color="0.15")
label_pad = 1.2
ax.annotate(
r"$\phi_{\max}$",
xy=tuple(major * label_pad),
ha="center",
va="center",
fontsize=9,
)
ax.annotate(
r"$\phi_{\min}$",
xy=tuple(minor * label_pad),
ha="center",
va="center",
fontsize=9,
)
ax.plot([0.0, phi_max * 1.05], [0.0, 0.0], "--", lw=0.7, color="0.6")
arc_r = 0.4 * phi_max
arc = np.linspace(0.0, theta, 30)
ax.plot(arc_r * np.cos(arc), arc_r * np.sin(arc), "-", lw=0.8, color="0.4")
ax.annotate(
r"$\theta$",
xy=(
arc_r * 1.4 * np.cos(theta / 2.0),
arc_r * 1.4 * np.sin(theta / 2.0),
),
ha="center",
va="center",
fontsize=9,
color="0.3",
)
lim = phi_max * 1.6
ax.set_xlim(-lim, lim)
ax.set_ylim(-lim, lim)
ax.set_aspect("equal", adjustable="box")
ax.axis("off")
return ax
# -------------------- small utilities (local) --------------------------- #
def _rho_det_from_z(z: np.ndarray, fr: np.ndarray) -> np.ndarray:
zx = z[:, 0, 1]
zy = z[:, 1, 0]
rx = 0.2 * (np.abs(zx) ** 2) / (fr + 1e-24)
ry = 0.2 * (np.abs(zy) ** 2) / (fr + 1e-24)
rdet = np.sqrt(rx * ry)
return rdet
def _phi_det_from_z(z: np.ndarray) -> np.ndarray:
d = z[:, 0, 0] * z[:, 1, 1] - z[:, 0, 1] * z[:, 1, 0]
return np.degrees(np.angle(d))
def _attach_det_metrics(
df: pd.DataFrame,
S,
) -> pd.DataFrame:
if df.empty:
return df
df = df.copy()
df["rho_det"] = np.nan
df["phi_det"] = np.nan
for i, ed in enumerate(_iter_items(S)):
st = _name(ed, i)
Z, z, fr = _get_z_block(ed)
if Z is None:
continue
rho = _rho_det_from_z(z, fr)
phi = _phi_det_from_z(z)
per = 1.0 / fr
m = df["station"] == st
if not m.any():
continue
p = df.loc[m, "period"].to_numpy()
idx = np.searchsorted(per, p)
idx = np.clip(idx, 0, len(per) - 1)
df.loc[m, "rho_det"] = rho[idx]
df.loc[m, "phi_det"] = phi[idx]
return df
def _color_from_hsv(h: np.ndarray, s: np.ndarray, v: np.ndarray) -> np.ndarray:
hsv = np.stack([h, s, v], axis=-1)
rgb = mcolors.hsv_to_rgb(hsv)
return rgb
# ---------------- 2) dimensionality grid (1D/2D/3D) --------------------- #
[docs]
def plot_dimensionality_grid(
sites: Any,
*,
skew_th: float = 3.0,
ellipt_th: float = 0.2,
figsize: tuple[float, float] = (8.5, 4.0),
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
ax: plt.Axes | None = None,
) -> plt.Axes:
df = build_phase_tensor_table(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
if df.empty:
if ax is None:
_, ax = plt.subplots(figsize=figsize)
ax.text(0.5, 0.5, "no phase tensor", ha="center", va="center")
return ax
df = df.copy()
df["logp"] = np.log10(df["period"].to_numpy())
lab = np.zeros(len(df), dtype=int)
a = np.abs(df["skew"].to_numpy())
e = np.abs(df["ellipt"].to_numpy())
lab[:] = 2
lab[(a <= skew_th) & (e <= ellipt_th)] = 0
lab[(a <= skew_th) & (e > ellipt_th)] = 1
df["dim"] = lab
piv = df.pivot_table(
index="logp",
columns="station",
values="dim",
aggfunc="median",
)
piv = piv.sort_index()
Z = piv.to_numpy(dtype=float) # (n_logp, n_st) — no transpose
if ax is None:
_, ax = plt.subplots(figsize=figsize)
im = ax.imshow(
Z,
aspect="auto",
origin="lower",
interpolation="nearest",
)
ax.set_ylabel(LOG10_PERIOD_LABEL)
PYCSAMT_STATION_RENDERING.apply(
ax,
np.arange(len(piv.columns), dtype=float),
list(piv.columns),
preset="pseudosection",
xlim=(-0.5, len(piv.columns) - 0.5),
)
yt = np.linspace(
0, Z.shape[0] - 1, num=min(8, Z.shape[0])
) # shape[0] = n_logp
yv = np.linspace(
piv.index.min(), piv.index.max(), num=min(8, len(piv.index))
)
ax.set_yticks(yt)
ax.set_yticklabels([f"{v:.2g}" for v in yv])
if not ax.yaxis_inverted():
ax.invert_yaxis()
cb = plt.colorbar(im, ax=ax)
cb.set_label("dim (0=1D,1=2D,2=3D)")
return ax
# ---------------- 3) theta stability stripe (HSV image) ----------------- #
[docs]
def plot_theta_stability_stripe(
sites: Any,
*,
win: int = 5,
figsize: tuple[float, float] = (9.0, 4.0),
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
ax: plt.Axes | None = None,
) -> plt.Axes:
df = build_phase_tensor_table(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
if df.empty:
if ax is None:
_, ax = plt.subplots(figsize=figsize)
ax.text(0.5, 0.5, "no phase tensor", ha="center", va="center")
return ax
df = df.copy()
df["logp"] = np.log10(df["period"].to_numpy())
sts = list(df["station"].unique())
X = []
H = []
for st in sts:
s = df[df["station"] == st].sort_values("logp")
th = s["theta"].to_numpy()
lp = s["logp"].to_numpy()
# hue from theta in [0,180) → [0,1)
h = (th % 180.0) / 180.0
# variance in sliding window
k = max(3, int(win))
if len(h) < k:
v = np.full_like(h, np.nan)
else:
v = np.convolve(
(h - np.nanmean(h)) ** 2,
np.ones(k) / k,
mode="same",
)
# saturation from 1 - norm(var)
v0 = (
np.nanpercentile(v[np.isfinite(v)], 5)
if np.isfinite(v).any()
else 0.0
)
v1 = (
np.nanpercentile(v[np.isfinite(v)], 95)
if np.isfinite(v).any()
else 1.0
)
s_sat = 1.0 - np.clip((v - v0) / (v1 - v0 + 1e-12), 0, 1)
H.append(np.vstack([h, s_sat, np.ones_like(h)]))
X.append(lp)
# align to common y-grid for imshow
yall = np.unique(np.concatenate(X))
Himg = np.zeros((len(yall), len(sts), 3))
for j, (lp, hs) in enumerate(zip(X, H)):
i = np.searchsorted(yall, lp)
i = np.clip(i, 0, len(yall) - 1)
h, s, v = hs
rgb = _color_from_hsv(h, s, v)
for r, k in enumerate(i):
Himg[k, j, :] = rgb[r]
if ax is None:
_, ax = plt.subplots(figsize=figsize)
ax.imshow(
Himg,
aspect="auto",
origin="lower",
interpolation="nearest",
)
ax.set_ylabel(LOG10_PERIOD_LABEL)
PYCSAMT_STATION_RENDERING.apply(
ax,
np.arange(len(sts)),
sts,
preset="pseudosection",
xlim=(-0.5, len(sts) - 0.5),
)
if not ax.yaxis_inverted():
ax.invert_yaxis()
yt = np.linspace(0, len(yall) - 1, num=min(8, len(yall)))
yv = np.linspace(yall.min(), yall.max(), num=min(8, len(yall)))
ax.set_yticks(yt)
ax.set_yticklabels([f"{v:.2g}" for v in yv])
return ax
# ---------------- 4) skew–ellipticity density (hexbin + contours) ------- #
[docs]
def plot_skew_ellipt_density(
sites: Any,
*,
band: tuple[float, float] | None = None,
gridsize: int = 40,
figsize: tuple[float, float] = (6.5, 5.5),
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
ax: plt.Axes | None = None,
) -> plt.Axes:
df = build_phase_tensor_table(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
if df.empty:
if ax is None:
_, ax = plt.subplots(figsize=figsize)
ax.text(0.5, 0.5, "no phase tensor", ha="center", va="center")
return ax
a = np.abs(df["beta"].to_numpy())
e = np.abs(df["ellipt"].to_numpy())
if band is not None:
lo, hi = band
m = (df["period"] >= lo) & (df["period"] <= hi)
a = a[m]
e = e[m]
if ax is None:
_, ax = plt.subplots(figsize=figsize)
ax.hexbin(a, e, gridsize=gridsize, mincnt=1, linewidths=0.0)
ax.set_xlabel("|beta| (deg)")
ax.set_ylabel("|ellipticity|")
# contours from 2D hist
nx = ny = gridsize
H, xedges, yedges = np.histogram2d(a, e, bins=[nx, ny])
Xc = 0.5 * (xedges[1:] + xedges[:-1])
Yc = 0.5 * (yedges[1:] + yedges[:-1])
lev = np.linspace(H.max() * 0.2, H.max() * 0.9, 4)
ax.contour(Xc, Yc, H.T, levels=lev, colors="k", linewidths=0.6)
return ax
# ---------------- 5) theta rose per band (small multiples) -------------- #
[docs]
def plot_theta_rose_grid(
sites: Any,
*,
n_bands: int = 6,
axes=None,
figsize: tuple[float, float] = (13.0, 3.8),
bins: int = 24,
style: str | RoseStyle | None = "pycsamt",
panel_title_fontsize: float = 7.5,
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
) -> plt.Figure:
"""Phase-tensor θ rose grid — one pycsamt-styled rose per frequency decade.
Draws *n_bands* polar rose diagrams side by side, each covering one
equal-log-width period band. Each rose applies axial symmetry (0–180°
mirrored to 180–360°) and is rendered using the active
:class:`~pycsamt.api._rose_style.RoseStyle`.
Parameters
----------
sites : any
Input accepted by :func:`~pycsamt.emtools._core.ensure_sites`.
n_bands : int, default ``6``
Number of equal-log-width period bands.
figsize : (width, height), default ``(13.0, 3.8)``
Figure size in inches. Constrained layout is used internally so
an external ``fig.suptitle`` fits without producing blank space.
bins : int, default ``24``
Number of bins over 0–180° (mirrored to 360°).
style : str, :class:`~pycsamt.api._rose_style.RoseStyle`, or None
Rose visual style. Strings resolved via
:func:`~pycsamt.api._rose_style.resolve_rose_style`.
panel_title_fontsize : float, default ``7.5``
Font size for the period-band label above each panel.
recursive, on_dup, strict, verbose
Passed to :func:`~pycsamt.emtools._core.ensure_sites`.
Returns
-------
:class:`matplotlib.figure.Figure`
"""
# ── resolve rose style ────────────────────────────────────────────────
rs = resolve_rose_style(style)
# ── axial mean helper ─────────────────────────────────────────────────
def _axial_mean(a_deg: np.ndarray) -> float:
rad = np.radians(2.0 * a_deg)
mu = (
np.degrees(
np.arctan2(np.nanmean(np.sin(rad)), np.nanmean(np.cos(rad)))
)
/ 2.0
)
return float(mu % 180.0)
df = build_phase_tensor_table(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
axes_given = _axes_list(axes, 1) if axes is not None else None
if df.empty:
if axes_given is None:
fig = plt.figure(figsize=figsize)
ax = fig.add_subplot(111, polar=True)
else:
ax = axes_given[0]
fig = ax.figure
ax.text(0.5, 0.5, "no phase tensor", ha="center", va="center")
return fig
p = df["period"].to_numpy()
lo_g = float(np.nanmin(p[p > 0])) if np.any(p > 0) else 1e-6
hi_g = float(np.nanmax(p))
edges = np.logspace(np.log10(lo_g), np.log10(hi_g), n_bands + 1)
bins_ = int(max(12, bins))
edges_deg = np.linspace(0.0, 180.0, bins_ + 1)
dw = np.radians(180.0 / bins_)
ang_mid = np.radians(0.5 * (edges_deg[1:] + edges_deg[:-1]))
# ── constrained_layout avoids blank space when suptitle is added ─────
axes_given = _axes_list(axes, n_bands) if axes is not None else None
fig = (
plt.figure(figsize=figsize, constrained_layout=True)
if axes_given is None
else axes_given[0].figure
)
for i in range(n_bands):
ax = (
axes_given[i]
if axes_given is not None
else fig.add_subplot(1, n_bands, i + 1, polar=True)
)
m = (p >= edges[i]) & (p < edges[i + 1])
th_deg = df.loc[m, "theta"].to_numpy(float) % 180.0
n_obs = int(np.sum(m))
# histogram — 0–180°, mirrored for axial symmetry
hist, _ = np.histogram(th_deg, bins=edges_deg)
rmax = float(hist.max()) if hist.max() > 0 else 1.0
# ── polar setup ───────────────────────────────────────────────────
ax.set_theta_zero_location("N")
ax.set_theta_direction(-1)
ax.set_yticklabels([])
ax.set_yticks([])
hide_polar_radius_labels(ax)
ax.yaxis.grid(False)
ax.xaxis.grid(False)
ax.set_frame_on(False)
# ── bars with axial symmetry (0–180° + mirror 180–360°) ──────────
if rs.bar_style == "gradient":
cm_ = plt.get_cmap(rs.cmap)
for a, h in zip(ang_mid, hist):
col = cm_(h / rmax) if rmax > 0 else cm_(0.5)
for a_plot in (a, a + np.pi):
ax.bar(
a_plot,
h,
width=dw,
color=col,
edgecolor=rs.bar_edgecolor,
linewidth=rs.bar_edgelw,
alpha=rs.bar_alpha,
)
else: # solid
for a, h in zip(ang_mid, hist):
for a_plot in (a, a + np.pi):
ax.bar(
a_plot,
h,
width=dw,
color=rs.bar_color,
edgecolor=rs.bar_edgecolor,
linewidth=rs.bar_edgelw,
alpha=rs.bar_alpha,
)
# ── concentric rings ──────────────────────────────────────────────
step = rmax / rs.n_rings
r_levels = [step * k for k in range(1, rs.n_rings + 1)]
theta_full = np.linspace(0, 2 * np.pi, 360)
for rv in r_levels:
ax.plot(
theta_full,
np.full_like(theta_full, rv),
color=rs.ring_color,
ls=rs.ring_ls,
lw=rs.ring_lw,
zorder=0,
)
# ── radial spokes ─────────────────────────────────────────────────
n_spokes = int(round(360.0 / rs.spoke_every))
spoke_angles = np.radians(np.arange(n_spokes) * rs.spoke_every)
for sa in spoke_angles:
ax.plot(
[sa, sa],
[0, rmax],
color=rs.spoke_color,
ls=rs.spoke_ls,
lw=rs.spoke_lw,
zorder=0,
)
# ── outer bold ring ────────────────────────────────────────────────
ax.spines["polar"].set_linewidth(rs.outer_ring_lw)
ax.spines["polar"].set_edgecolor(rs.outer_ring_color)
ax.set_ylim(0, rmax * 1.08)
# ── compass / degree labels ───────────────────────────────────────
spoke_degs = np.arange(n_spokes) * rs.spoke_every
if rs.compass_labels == "NESW":
_card = {0: "N", 90: "E", 180: "S", 270: "W"}
lbl_list = [_card.get(int(d) % 360, "") for d in spoke_degs]
elif rs.compass_labels == "degrees":
lbl_list = [f"{int(d)}°" for d in spoke_degs]
else:
lbl_list = [""] * n_spokes
ax.set_thetagrids(
spoke_degs, labels=lbl_list, fontsize=rs.compass_fontsize
)
for lbl in ax.get_xticklabels():
lbl.set_color(rs.compass_color)
lbl.set_fontweight(rs.compass_fontweight)
# ── mean direction line ───────────────────────────────────────────
if rs.show_mean and n_obs > 0:
mu = _axial_mean(th_deg)
mu_r = np.radians(mu)
sec_lw = rs.mean_lw if rs.secondary_lw is None else rs.secondary_lw
sec_col = (
rs.mean_color
if rs.secondary_color is None
else rs.secondary_color
)
for ang_pt in (mu_r, mu_r + np.pi):
ax.plot(
[ang_pt, ang_pt],
[0, rmax],
color=rs.mean_color,
lw=rs.mean_lw,
ls=rs.mean_ls,
zorder=6,
)
if rs.show_secondary:
for ang_pt in (mu_r + np.pi / 2, mu_r - np.pi / 2):
ax.plot(
[ang_pt, ang_pt],
[0, rmax],
color=sec_col,
lw=sec_lw,
ls=rs.secondary_ls,
zorder=6,
)
# ── panel annotation: period range + n ────────────────────────────
if rs.show_annotation and n_obs > 0:
mu_str = ""
if rs.show_mean and n_obs > 0:
mu = _axial_mean(th_deg)
mu_str = f"θ̄={mu:.0f}°"
n_str = f" n={n_obs}" if rs.show_n else ""
ax.text(
*rs.annotation_pos,
mu_str + n_str,
transform=ax.transAxes,
fontsize=rs.annotation_fontsize,
ha="left",
va="top",
bbox=dict(
fc=rs.annotation_bg,
ec=rs.annotation_ec,
alpha=0.80,
pad=1.5,
boxstyle="round,pad=0.3",
),
zorder=10,
)
# ── panel title: period band ───────────────────────────────────────
ax.set_title(
f"[{edges[i]:.2g}, {edges[i + 1]:.2g}]s",
fontsize=panel_title_fontsize,
pad=8,
)
return fig
# ─────────────────────────────────────────────────────────────────────────────
# Per-station ellipse strip + multi-profile small-multiples grid
# ─────────────────────────────────────────────────────────────────────────────
def _pt_size_reference(
s1_vals: np.ndarray,
normalise_by: str,
s1_ref: float | None,
) -> float | None:
"""Resolve the ellipse-size reference for ``"cell"``/``"unity"``/``"abs"``.
Mirrors the sizing rule used by :func:`plot_phase_tensor_psection` so
that all phase-tensor ellipse views scale identically for a given
*normalise_by* / *s1_ref* combination.
"""
if normalise_by == "cell":
if s1_ref is not None:
return max(float(s1_ref), 1e-9)
finite_s1 = s1_vals[np.isfinite(s1_vals)]
ref = float(np.nanpercentile(finite_s1, 90)) if len(finite_s1) else 1.0
return max(ref, 1e-9)
if normalise_by == "unity":
return max(float(s1_ref) if s1_ref is not None else 1.0, 1e-9)
return None # "abs" — raw data units
[docs]
def plot_phase_tensor_strip(
sites: Any,
*,
# ── data selection ────────────────────────────────────────────────────
station: str | None = None,
period_range: tuple[float, float] | None = None,
# ── ellipse sizing ────────────────────────────────────────────────────
scale=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.scale
normalise_by=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.normalise_by
s1_ref: float | None = None,
min_aspect=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.min_aspect
cells_per_decade: float = 7.0,
# ── colour ────────────────────────────────────────────────────────────
c_by=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.c_by
cmap=_UNSET, # default: auto from c_by via resolve_cmap()
clim: tuple[float, float] | None = None,
clim_pct=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.clim_pct
symmetric_clim=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.resolve_symmetric_clim()
# ── aesthetics ────────────────────────────────────────────────────────
edgecolor=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.edgecolor
linewidth=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.linewidth
alpha=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.alpha
# ── overlays ─────────────────────────────────────────────────────────
skew_threshold=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.skew_threshold
mark_3d=_UNSET, # default: PYCSAMT_STYLE.pt_ellipse.mark_3d
# ── decorative phase scale (y-axis) ───────────────────────────────────
phase_ticks: tuple[float, float, float] | None = (0.0, 45.0, 90.0),
ylabel: str = "",
# ── station corner label ──────────────────────────────────────────────
station_label: bool = True,
station_label_fontsize: float = 8.0,
# ── labels & layout ───────────────────────────────────────────────────
title: str = "",
xlabel: str = "",
axis_style: str = "period",
figsize: tuple[float, float] = (6.0, 1.4),
show_colorbar: bool = True,
colorbar_label: str | None = None,
# ── standard emtools ──────────────────────────────────────────────────
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
ax: plt.Axes | None = None,
) -> plt.Axes:
"""Single-station phase-tensor ellipse strip vs period.
Draws one horizontal row of phase-tensor ellipses for a single
station, one ellipse per period, in the classic "ellipse timeseries"
style used e.g. by WinGLink / mtpy single-station PT plots. The
y-axis carries no per-point data — it is a schematic 0°–90° phase
scale (*phase_ticks*) used only to calibrate ellipse size by eye; the
physical information (φ_max, φ_min, θ, fill colour) is identical to
:func:`plot_phase_tensor_psection`.
Combine several stations/profiles into the small-multiples layout
(columns = profiles, rows = stations, one shared colorbar) with
:func:`plot_phase_tensor_strip_grid`.
Parameters
----------
sites : any
EDI paths, glob pattern, :class:`~pycsamt.site.base.Sites`, any
input accepted by :func:`~pycsamt.emtools.ensure_sites`, or a
pre-built :class:`pandas.DataFrame` from
:func:`build_phase_tensor_table` (already filtered to one
station — in that case *station* may be left ``None``).
station : str or None
Station to plot. Required when *sites* resolves to more than
one station.
period_range : (T_min, T_max) in seconds or None
Restrict the period range plotted.
scale, normalise_by, s1_ref, min_aspect : see
:func:`plot_phase_tensor_psection`. Sizing controls;
``normalise_by="unity"`` (φ=45° 1-D reference fills *scale* of the
row) makes *phase_ticks* an exact scale. *min_aspect* floors the
ellipse height so near-degenerate cells stay visible.
cells_per_decade : float, default ``7.0``
Visual ellipse pitch along the period axis, expressed as the
number of ellipse-widths per log10 decade. Unlike
:func:`plot_phase_tensor_psection` (one ellipse per station
column), a period axis is typically sampled far more densely
than is useful for ellipse width — sizing each ellipse to its
*local* sample spacing would shrink it to an invisible sliver
wherever sampling is dense. This value is independent of how
many periods were actually measured, so ellipses overlap
(by design — the classic single-station "ellipse timeseries"
look) rather than shrinking as sampling gets denser.
c_by, cmap, clim, clim_pct, symmetric_clim :
Fill-colour controls; see :func:`plot_phase_tensor_psection`.
edgecolor, linewidth, alpha : ellipse border/opacity controls.
skew_threshold, mark_3d : 3-D cell highlighting, see
:func:`plot_phase_tensor_psection`.
phase_ticks : (lo, mid, hi) or None, default ``(0.0, 45.0, 90.0)``
Tick values drawn on the schematic y-scale. ``None`` hides the
y-axis entirely.
ylabel : str
Y-axis label; defaults to ``"Phase (°)"`` when *phase_ticks* is
not ``None``.
station_label : bool, default ``True``
Annotate the station name in the upper-left corner of the axes.
station_label_fontsize : float, default ``8.0``
title, xlabel : str
Axes title / x-label. *xlabel* defaults to ``"Period (s)"``
(or :data:`~pycsamt.api.labels.LOG10_PERIOD_LABEL` when
*axis_style* is ``"logperiod"``).
axis_style : {"period", "logperiod"}, default ``"period"``
Tick-label convention for the internally log10(period)-spaced
x-axis. ``"period"`` (default, unchanged) formats ticks as
``$10^{n}$`` under a ``"Period (s)"`` label. ``"logperiod"``
shows the raw log10 values (``-3``, ``-2``, ...) under
:data:`~pycsamt.api.labels.LOG10_PERIOD_LABEL`, matching
:func:`plot_phase_tensor_psection`'s own convention.
figsize : (float, float), default ``(6.0, 1.4)``
Figure size (ignored when *ax* is provided).
show_colorbar : bool, default ``True``
Attach a right-side colorbar. Set ``False`` when composing a
grid with a single shared colorbar (see
:func:`plot_phase_tensor_strip_grid`).
colorbar_label : str or None
Override the automatic colorbar label derived from *c_by*.
recursive, on_dup, strict, verbose : see :func:`ensure_sites`.
ax : :class:`~matplotlib.axes.Axes` or None
If provided, draw into this axes and return it; otherwise create
a new figure.
Returns
-------
ax : :class:`~matplotlib.axes.Axes`
Examples
--------
>>> from pycsamt.emtools import plot_phase_tensor_strip
>>> ax = plot_phase_tensor_strip(sites, station="S1")
Fixed colour scale (for visual consistency across several calls):
>>> ax = plot_phase_tensor_strip(
... sites,
... station="S1",
... c_by="skew",
... clim=(-9.0, 9.0),
... )
See Also
--------
plot_phase_tensor_strip_grid : Multi-station / multi-profile facet grid.
plot_phase_tensor_psection : Station × period pseudo-section.
"""
if isinstance(sites, pd.DataFrame):
df = sites.copy()
else:
df = build_phase_tensor_table(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
# ── resolve visual style from PYCSAMT_STYLE.pt_ellipse ──────────────
_es = PYCSAMT_STYLE.pt_ellipse
if scale is _UNSET:
scale = _es.scale
if normalise_by is _UNSET:
normalise_by = _es.normalise_by
if min_aspect is _UNSET:
min_aspect = _es.min_aspect
if c_by is _UNSET:
c_by = _es.c_by
if cmap is _UNSET:
cmap = _es.copy(c_by=c_by).resolve_cmap()
if clim_pct is _UNSET:
clim_pct = _es.clim_pct
if symmetric_clim is _UNSET:
symmetric_clim = _es.copy(c_by=c_by).resolve_symmetric_clim()
if edgecolor is _UNSET:
edgecolor = _es.edgecolor
if linewidth is _UNSET:
linewidth = _es.linewidth
if alpha is _UNSET:
alpha = _es.alpha
if skew_threshold is _UNSET:
skew_threshold = _es.skew_threshold
if mark_3d is _UNSET:
mark_3d = _es.mark_3d
# ── data selection ────────────────────────────────────────────────────
if not df.empty and period_range is not None:
lo, hi = float(period_range[0]), float(period_range[1])
df = df[(df["period"] >= lo) & (df["period"] <= hi)]
if not df.empty:
stations_present = list(dict.fromkeys(df["station"].tolist()))
if station is not None:
df = df[df["station"] == station]
elif len(stations_present) > 1:
raise ValueError(
"plot_phase_tensor_strip() needs a single station; "
f"got {len(stations_present)} — pass station=... "
"(or use plot_phase_tensor_strip_grid for many stations)."
)
else:
station = stations_present[0] if stations_present else station
if ax is None:
_, ax = plt.subplots(figsize=figsize)
if df.empty:
ax.text(
0.5,
0.5,
"no phase tensor data",
ha="center",
va="center",
transform=ax.transAxes,
)
return ax
df = df.sort_values("period")
# ── x-axis: log10(period); y-axis: single schematic row ──────────────
# cell_x is a *visual* pitch (ellipse-widths per decade), not the raw
# sample spacing — see `cells_per_decade` docstring for why.
x_all = np.log10(df["period"].to_numpy(float))
cell_x = 1.0 / max(cells_per_decade, 1e-6)
cell_y = 1.0
s1_vals = df["s1"].to_numpy(float)
ref = _pt_size_reference(s1_vals, normalise_by, s1_ref)
cvals, cbar_label = _resolve_cvals(df, c_by)
if colorbar_label:
cbar_label = colorbar_label
finite_c = cvals[np.isfinite(cvals)]
if clim is not None:
vmin, vmax = float(clim[0]), float(clim[1])
elif skew_threshold is not None and c_by in ("skew", "beta"):
vmin, vmax = -float(skew_threshold), float(skew_threshold)
else:
lo_pct, hi_pct = clim_pct
vmin = (
float(np.nanpercentile(finite_c, lo_pct))
if len(finite_c)
else -1.0
)
vmax = (
float(np.nanpercentile(finite_c, hi_pct)) if len(finite_c) else 1.0
)
if symmetric_clim and c_by in _SYMMETRIC_C:
vlim = max(abs(vmin), abs(vmax))
vmin, vmax = -vlim, vlim
norm = Normalize(vmin=vmin, vmax=vmax)
cm = plt.get_cmap(cmap)
s1_arr = df["s1"].to_numpy(float)
s2_arr = df["s2"].to_numpy(float)
th_arr = df["theta"].to_numpy(float)
skew_arr = (
df["skew"].to_numpy(float)
if "skew" in df.columns
else np.zeros(len(df))
)
for xi, s1i, s2i, thi, ci, beti in zip(
x_all, s1_arr, s2_arr, th_arr, cvals, skew_arr
):
if not np.all(np.isfinite([xi, s1i, s2i, thi, ci])):
continue
if ref is not None:
w = min(scale * cell_x * s1i / ref, 0.99 * cell_x)
h = min(scale * cell_y * s2i / ref, 0.99 * cell_y)
else:
w = scale * s1i
h = scale * s2i
if min_aspect > 0:
h = max(h, min_aspect * w)
is_3d = (
skew_threshold is not None
and mark_3d
and np.isfinite(beti)
and abs(beti) > skew_threshold
)
lw_i = linewidth * 3.0 if is_3d else linewidth
ax.add_patch(
Ellipse(
(xi, 0.0),
width=w,
height=h,
angle=thi,
facecolor=cm(norm(ci)),
edgecolor=edgecolor,
linewidth=lw_i,
alpha=alpha,
zorder=3,
)
)
# ── x-axis limits & ticks (log10 period; label convention per
# *axis_style*) ─────────────────────────────────────────────────────
if axis_style not in ("period", "logperiod"):
raise ValueError(
f"axis_style must be 'period' or 'logperiod'; got {axis_style!r}."
)
x_lo = float(np.nanmin(x_all)) if len(x_all) else 0.0
x_hi = float(np.nanmax(x_all)) if len(x_all) else 1.0
margin = max(0.5 * cell_x, 0.02 * (x_hi - x_lo + 1e-9))
ax.set_xlim(x_lo - margin, x_hi + margin)
x_int = np.arange(int(np.floor(x_lo)), int(np.ceil(x_hi)) + 1)
ax.set_xticks(x_int)
if axis_style == "logperiod":
ax.set_xticklabels([f"{int(v)}" for v in x_int], fontsize=8)
ax.set_xlabel(xlabel or LOG10_PERIOD_LABEL, fontsize=9)
else:
ax.set_xticklabels([f"$10^{{{int(v)}}}$" for v in x_int], fontsize=8)
ax.set_xlabel(xlabel or "Period (s)", fontsize=9)
# ── y-axis: schematic phase scale, no real data ───────────────────────
ax.set_ylim(-0.62, 0.62)
if phase_ticks is not None:
lo_t, mid_t, hi_t = phase_ticks
ax.set_yticks([-0.5, 0.0, 0.5])
ax.set_yticklabels(
[f"{lo_t:g}", f"{mid_t:g}", f"{hi_t:g}"], fontsize=8
)
ax.set_ylabel(ylabel or "Phase (°)", fontsize=9)
else:
ax.set_yticks([])
ax.set_ylabel(ylabel, fontsize=9)
if title:
ax.set_title(title, fontsize=9)
if station_label and station:
ax.text(
0.015,
0.92,
str(station),
transform=ax.transAxes,
fontsize=station_label_fontsize,
fontweight="bold",
ha="left",
va="top",
bbox=dict(fc="white", ec="none", alpha=0.65, pad=1.0),
)
ax.grid(True, ls=":", lw=0.4, color="#cccccc", alpha=0.7, zorder=0)
if show_colorbar:
sm = ScalarMappable(cmap=cm, norm=norm)
sm.set_array([])
_attach_cbar(ax, sm, cbar_label)
return ax
[docs]
def plot_phase_tensor_strip_grid(
sites: Any,
profiles: dict[str, list[str]],
*,
period_range: tuple[float, float] | None = None,
scale=_UNSET,
normalise_by=_UNSET,
s1_ref: float | None = None,
min_aspect=_UNSET,
cells_per_decade: float = 7.0,
c_by=_UNSET,
cmap=_UNSET,
clim: tuple[float, float] | None = None,
clim_pct=_UNSET,
symmetric_clim=_UNSET,
edgecolor=_UNSET,
linewidth=_UNSET,
alpha=_UNSET,
skew_threshold=_UNSET,
mark_3d=_UNSET,
phase_ticks: tuple[float, float, float] | None = (0.0, 45.0, 90.0),
col_titles: dict[str, str] | None = None,
xlabel: str = "Period (s)",
suptitle: str = "",
colorbar_label: str | None = None,
panel_size: tuple[float, float] = (4.4, 1.05),
recursive: bool = True,
on_dup: str = "replace",
strict: bool = False,
verbose: int = 0,
axes: Any | None = None,
) -> plt.Figure:
"""Phase-tensor ellipse strips for several stations, grouped by profile.
Reproduces the classic multi-panel figure — one
:func:`plot_phase_tensor_strip` row per selected station, tiled in a
grid where each *column* is a profile/line and each *row* is one of
the stations picked for that profile — with a **single shared
colorbar** for the whole figure (so fill colours are comparable
across every panel).
Parameters
----------
sites : any
Input accepted by :func:`~pycsamt.emtools._core.ensure_sites`
covering every station referenced in *profiles*.
profiles : dict[str, list[str]]
Mapping ``{profile_label: [station, ...]}``, e.g.::
{
"Profile L1": ["S1", "S15", "S30", "S45"],
"Profile L3": ["S1", "S4", "S7", "S11"],
}
Each list becomes one column; rows are aligned by position, not
by station identity (profiles may list a different number of
stations — empty cells are left blank).
period_range : (T_min, T_max) or None
Restrict the period range plotted in every panel.
scale, normalise_by, s1_ref, min_aspect, cells_per_decade, c_by, cmap,
clim, clim_pct, symmetric_clim, edgecolor, linewidth, alpha,
skew_threshold, mark_3d :
Forwarded to every :func:`plot_phase_tensor_strip` call; see its
docstring (and :func:`plot_phase_tensor_psection`) for details.
When *clim* is ``None`` it is derived once from the pooled data
of every listed station, so all panels share one colour scale.
phase_ticks : (lo, mid, hi) or None, default ``(0.0, 45.0, 90.0)``
Schematic y-scale ticks drawn on every panel.
col_titles : dict[str, str] or None
Override the column title text; defaults to the *profiles* keys.
xlabel : str, default ``"Period (s)"``
Shown once under the bottom-most panel of each column.
suptitle : str
Figure-level title.
colorbar_label : str or None
Override the automatic colorbar label derived from *c_by*.
panel_size : (width, height), default ``(4.4, 1.05)``
Per-panel size in inches; the figure size is
``(n_cols * width, n_rows * height)``.
recursive, on_dup, strict, verbose : see :func:`ensure_sites`.
axes : 2-D array of :class:`~matplotlib.axes.Axes` or None
Pre-existing ``(n_rows, n_cols)`` axes grid to draw into
(``n_rows = max(len(v) for v in profiles.values())``,
``n_cols = len(profiles)``). A new figure is created when
``None``.
Returns
-------
fig : :class:`~matplotlib.figure.Figure`
Examples
--------
>>> from pycsamt.emtools import plot_phase_tensor_strip_grid
>>> fig = plot_phase_tensor_strip_grid(
... sites,
... profiles={
... "Profile L1": ["S1", "S15", "S30", "S45"],
... "Profile L3": ["S1", "S4", "S7", "S11"],
... },
... c_by="skew",
... cmap="RdBu_r",
... )
See Also
--------
plot_phase_tensor_strip : Single-station ellipse strip (one panel).
plot_phase_tensor_psection : Station × period pseudo-section.
"""
# ── resolve visual style from PYCSAMT_STYLE.pt_ellipse ──────────────
_es = PYCSAMT_STYLE.pt_ellipse
if scale is _UNSET:
scale = _es.scale
if normalise_by is _UNSET:
normalise_by = _es.normalise_by
if min_aspect is _UNSET:
min_aspect = _es.min_aspect
if c_by is _UNSET:
c_by = _es.c_by
if cmap is _UNSET:
cmap = _es.copy(c_by=c_by).resolve_cmap()
if clim_pct is _UNSET:
clim_pct = _es.clim_pct
if symmetric_clim is _UNSET:
symmetric_clim = _es.copy(c_by=c_by).resolve_symmetric_clim()
if edgecolor is _UNSET:
edgecolor = _es.edgecolor
if linewidth is _UNSET:
linewidth = _es.linewidth
if alpha is _UNSET:
alpha = _es.alpha
if skew_threshold is _UNSET:
skew_threshold = _es.skew_threshold
if mark_3d is _UNSET:
mark_3d = _es.mark_3d
all_stations = list(
dict.fromkeys(st for st_list in profiles.values() for st in st_list)
)
df_all = build_phase_tensor_table(
sites,
recursive=recursive,
on_dup=on_dup,
strict=strict,
verbose=verbose,
)
if not df_all.empty:
df_all = df_all[df_all["station"].isin(all_stations)]
if period_range is not None:
lo, hi = float(period_range[0]), float(period_range[1])
df_all = df_all[
(df_all["period"] >= lo) & (df_all["period"] <= hi)
]
# ── shared colour scale across every panel ────────────────────────────
if clim is not None:
clim_shared = clim
elif not df_all.empty:
cvals_all, _ = _resolve_cvals(df_all, c_by)
finite_all = cvals_all[np.isfinite(cvals_all)]
if skew_threshold is not None and c_by in ("skew", "beta"):
clim_shared = (-float(skew_threshold), float(skew_threshold))
else:
lo_pct, hi_pct = clim_pct
vmin = (
float(np.nanpercentile(finite_all, lo_pct))
if len(finite_all)
else -1.0
)
vmax = (
float(np.nanpercentile(finite_all, hi_pct))
if len(finite_all)
else 1.0
)
if symmetric_clim and c_by in _SYMMETRIC_C:
vlim = max(abs(vmin), abs(vmax))
vmin, vmax = -vlim, vlim
clim_shared = (vmin, vmax)
else:
clim_shared = (-1.0, 1.0)
n_cols = len(profiles)
n_rows = max((len(v) for v in profiles.values()), default=1)
if axes is None:
fig, axes = plt.subplots(
n_rows,
n_cols,
squeeze=False,
sharex="col",
figsize=(panel_size[0] * n_cols, panel_size[1] * n_rows),
constrained_layout=True,
)
else:
axes = np.atleast_2d(axes)
fig = axes[0, 0].figure
col_titles = col_titles or {}
_, cbar_label = (
_resolve_cvals(df_all, c_by) if not df_all.empty else (None, c_by)
)
if colorbar_label:
cbar_label = colorbar_label
for j, (prof_name, st_list) in enumerate(profiles.items()):
for i in range(n_rows):
ax = axes[i, j]
if i >= len(st_list):
ax.axis("off")
continue
st = st_list[i]
plot_phase_tensor_strip(
df_all,
station=st,
scale=scale,
normalise_by=normalise_by,
s1_ref=s1_ref,
min_aspect=min_aspect,
cells_per_decade=cells_per_decade,
c_by=c_by,
cmap=cmap,
clim=clim_shared,
clim_pct=clim_pct,
symmetric_clim=symmetric_clim,
edgecolor=edgecolor,
linewidth=linewidth,
alpha=alpha,
skew_threshold=skew_threshold,
mark_3d=mark_3d,
phase_ticks=phase_ticks,
show_colorbar=False,
xlabel=xlabel if i == len(st_list) - 1 else "",
ax=ax,
)
if i != len(st_list) - 1:
ax.tick_params(labelbottom=False)
ax.set_xlabel("")
axes[0, j].set_title(col_titles.get(prof_name, prof_name), fontsize=10)
cm = plt.get_cmap(cmap)
norm = Normalize(vmin=clim_shared[0], vmax=clim_shared[1])
sm = ScalarMappable(cmap=cm, norm=norm)
sm.set_array([])
cb = fig.colorbar(sm, ax=axes, shrink=0.9, pad=0.015, aspect=30)
cb.set_label(cbar_label, fontsize=9)
cb.ax.tick_params(labelsize=8)
if suptitle:
fig.suptitle(suptitle, fontsize=11)
return fig
