"""Minimal plotting helpers for DASdata.
Two free functions (`imshow`, `wiggle`) and a bound `_PlotAccessor`
exposed via `DASdata.plot`. Matplotlib is imported lazily inside each
function body so a bare `import dasio` (or `import dasdata`) stays
free of the matplotlib startup cost — the load only fires when an
actual plot is requested.
`plot_xcorr` images a cross-correlation gather (the `cc` array of a
`dasio.noise.CCGather`); `CCGather.plot()` is the convenience wrapper
around it. `plot_spectrogram` images the `(freq, times, S)` result of
`psd.spectrogram` / `psd.spectrogram_cwt`, and `plot_fk` / `plot_fc`
the f-k and dispersion panels from `dasio.fk`.
"""
from typing import Optional, Tuple
import numpy as np
[docs]
def imshow(
d,
ax=None,
ch_range: Optional[Tuple[int, int]] = None,
t_range: Optional[Tuple] = None,
skip_ch: int = 1,
skip_t: int = 1,
style: str = 'seismic',
usedatetime: bool = False,
perc: float = 99.5,
vmin: Optional[float] = None,
vmax: Optional[float] = None,
cmap: str = 'seismic',
cbar: bool = True,
figsize: Optional[Tuple[float, float]] = None,
**kwargs,
):
"""Color-mapped image of a `DASdata` (regular grid).
Parameters
----------
d : DASdata
ax : matplotlib Axes, optional
Plot target; one is created if not given.
ch_range, t_range : tuples, optional
Forwarded to `DASdata.truncate` for the displayed window.
For an arbitrary channel subset, pre-select with
`d.select(ch_index=...)` and plot the result.
skip_ch, skip_t : int, default 1
Plot only every `skip_ch`-th channel / `skip_t`-th time sample —
a coarser, faster raster for big arrays (e.g. `skip_t=5` on a wide
window). Reuses `DASdata.skip_ch` / `.skip_t`. The percentile and
the image both run on the decimated array, so this genuinely speeds
the plot up; the channel/time axes still span the full range. `1`
(default) draws every sample.
style : {'seismic', 'normal'}
`'seismic'` (default, matches legacy DASutils) puts channels
on the x-axis and time on the y-axis with time growing
downward — the traditional seismic-record orientation.
`'normal'` puts time on the x-axis and channels on the y-axis
— easier to read as a stack of time-series.
usedatetime : bool
Show the time axis as `datetime64` ticks (True) instead of
seconds-from-`t0_sec` (False, default).
perc : float
Percentile (of `|data|`) used to set symmetric `vmin` / `vmax`
when the latter aren't provided. Default 99.5 trims out the
loudest 0.5 % so background noise stays visible.
vmin, vmax : floats, optional
Manual color-scale bounds; bypass `perc` when both are given.
cmap : str or Colormap
Default `'seismic'` (red-white-blue), matches legacy DASutils.
cbar : bool
figsize : (w, h), optional
Figure size when `ax` is created here. Defaults to (8, 4)
for `style='seismic'` and (4, 8) for `style='normal'`.
Ignored when `ax` is supplied.
**kwargs : forwarded to `ax.imshow`
Returns
-------
(ax, im) so callers can attach annotations / set titles / etc.
"""
import matplotlib.pyplot as plt
from matplotlib import dates as mdates
if style not in ('seismic', 'normal'):
raise ValueError(f"style must be 'seismic' or 'normal', got {style!r}")
sub = d if (ch_range is None and t_range is None) else d.truncate(
ch_range=ch_range, t_range=t_range,
)
# Channel count and time-axis extent come from the *full* window, before
# any skip, so the axes always label the true channel / time range; the
# skip below only coarsens the rasterized array, not the coordinates.
nx_full = sub.nx
# Channel numbers of the first and last row, captured before `skip_ch`
# below multiplies `dch`. On a min_ch=2000 read these are 2000 and up, so
# the axis labels the fiber rather than the array. From `channels()`, the
# authority — `ch0 + i*dch` only matches it while the rows are a ramp.
ch_axis = sub.channels()
ch_lo, ch_hi = int(ch_axis[0]), int(ch_axis[-1])
if usedatetime:
dt_axis = sub.times(type='datetime')
t_left = mdates.date2num(dt_axis[0].astype('datetime64[us]').astype(object))
t_right = mdates.date2num(dt_axis[-1].astype('datetime64[us]').astype(object))
else:
tax = sub.times()
t_left, t_right = float(tax[0]), float(tax[-1])
# Decimate the raster — every `skip_ch`-th channel and `skip_t`-th sample,
# reusing DASdata.skip_ch / .skip_t. The percentile and imshow then run on
# the smaller array (the real speedup); the extent above keeps the axes.
sub = sub.skip_ch(skip_ch).skip_t(skip_t)
if vmin is None or vmax is None:
v = float(np.nanpercentile(np.abs(sub.data), perc))
if vmin is None:
vmin = -v
if vmax is None:
vmax = v
if ax is None:
if figsize is None:
figsize = (8, 4) if style == 'seismic' else (4, 8)
_, ax = plt.subplots(figsize=figsize)
if style == 'normal':
# data shape (nx, nt) — imshow it directly: rows=channels, cols=time.
# extent = (left, right, bottom, top); first channel at top.
extent = (t_left, t_right, ch_hi, ch_lo)
im = ax.imshow(
sub.data, aspect='auto', cmap=cmap,
vmin=vmin, vmax=vmax, extent=extent, **kwargs,
)
ax.set_xlabel('time' if usedatetime else 't (s)')
ax.set_ylabel('channel')
date_axis = ax.xaxis if usedatetime else None
else:
# 'seismic': transpose so rows=time, cols=channels; flip y so
# time grows downward (small-time at top, large-time at bottom).
extent = (ch_lo, ch_hi, t_right, t_left)
im = ax.imshow(
sub.data.T, aspect='auto', cmap=cmap,
vmin=vmin, vmax=vmax, extent=extent, **kwargs,
)
ax.set_xlabel('channel')
ax.set_ylabel('time' if usedatetime else 't (s)')
date_axis = ax.yaxis if usedatetime else None
if date_axis is not None:
date_axis.axis_date() if hasattr(date_axis, 'axis_date') else None
date_axis.set_major_locator(mdates.AutoDateLocator())
date_axis.set_major_formatter(
mdates.AutoDateFormatter(date_axis.get_major_locator())
)
if cbar:
ax.figure.colorbar(im, ax=ax, label=getattr(d, 'units', None))
return ax, im
[docs]
def wiggle(
d,
ax=None,
ch_range: Optional[Tuple[int, int]] = None,
t_range: Optional[Tuple] = None,
n_max_ch: Optional[int] = 100,
max_nt: Optional[int] = 5000,
style: str = 'seismic',
scale: float = 1.0,
color: str = 'k',
lw: float = 0.6,
usedatetime: bool = False,
normalize: bool = True,
figsize: Optional[Tuple[float, float]] = None,
**kwargs,
):
"""Per-channel wiggle traces of a `DASdata`.
Each retained channel becomes a line offset by its (original)
channel index, optionally amplitude-normalized so all traces
share the same visual range.
`style='seismic'` (default) places channels on the x-axis and
time on the y-axis growing downward — the traditional seismic
record orientation. `style='normal'` places time on the x-axis
and channel on the y-axis (descending), which reads more like
a stack of time-series.
`n_max_ch` caps the number of traces drawn — when
`sub.nx > n_max_ch` the channels are decimated with stride
`ceil(sub.nx / n_max_ch)`. Defaults to 100; pass `None` to draw
every channel. This stops a casual `d.plot.wiggle()` on a
full-array DASdata from rendering thousands of overlapping lines.
Use `ch_range` first if you want a specific window.
`max_nt` caps the number of time samples drawn *per trace* — when
`sub.nt > max_nt` the time axis is strided by `ceil(sub.nt / max_nt)`.
Wiggle render cost scales with `n_traces * n_samples`, and at screen
resolution a few thousand points per trace is already pixel-lossless,
so this is the main speed lever on wide time windows. Defaults to
5000; pass `None` to draw every sample.
"""
import matplotlib.pyplot as plt
if style not in ('seismic', 'normal'):
raise ValueError(f"style must be 'seismic' or 'normal', got {style!r}")
sub = d if (ch_range is None and t_range is None) else d.truncate(
ch_range=ch_range, t_range=t_range,
)
if n_max_ch is not None and sub.nx > n_max_ch:
step = int(np.ceil(sub.nx / n_max_ch))
ch_indices = np.arange(0, sub.nx, step)
else:
step = 1
ch_indices = np.arange(sub.nx)
# Decimate the time axis for display. Wiggle cost scales with
# n_traces * n_samples; at screen resolution a few thousand points per
# trace is already pixel-lossless, so striding wide windows is the main
# speedup (analogous to n_max_ch for channels).
if max_nt is not None and sub.nt > max_nt:
t_step = int(np.ceil(sub.nt / max_nt))
else:
t_step = 1
# Rows select the traces; channel numbers place them on the axis.
ch_vals = sub.channels()[ch_indices]
data = sub.data[ch_indices][:, ::t_step]
if normalize:
peak = np.maximum(np.abs(data).max(axis=1, keepdims=True), 1e-30)
data = data / peak
# Scale by the plotted spacing (`step` traces apart, `dch` channels each)
# so a peak ±1 normalized trace just fills the gap between neighbours —
# keeps visual density constant regardless of decimation.
gap = step * sub.dch
data = data * scale * gap
from matplotlib.collections import LineCollection
# Numeric time axis for the line geometry. A LineCollection bypasses
# matplotlib's datetime unit conversion, so convert explicitly here and
# flip the relevant axis into date mode afterwards.
if usedatetime:
from matplotlib import dates as mdates
t_vals = mdates.date2num(
sub.times(type='datetime')[::t_step].astype('datetime64[us]').astype(object)
)
else:
t_vals = sub.times()[::t_step]
if ax is None:
if figsize is None:
n_ch = len(ch_indices)
figsize = (
(max(2, n_ch * 0.25), 10) if style == 'seismic'
else (10, max(2, n_ch * 0.25))
)
_, ax = plt.subplots(figsize=figsize)
# One LineCollection (single artist) instead of an `ax.plot` per channel —
# the per-trace Python loop is the wiggle bottleneck on wide windows.
offs = data + ch_vals[:, None] # (n_lines, nt) offset traces
tt = np.broadcast_to(t_vals, offs.shape) # (n_lines, nt)
if style == 'seismic':
segs = np.stack([offs, tt], axis=-1) # x = amplitude, y = time
else:
segs = np.stack([tt, offs], axis=-1) # x = time, y = amplitude
ax.add_collection(LineCollection(segs, colors=color, linewidths=lw, **kwargs))
if style == 'seismic':
ax.set_xlim(ch_vals[0] - gap, ch_vals[-1] + gap)
ax.set_ylim(t_vals[-1], t_vals[0]) # inverted (time grows down)
ax.set_xlabel('channel')
ax.set_ylabel('time' if usedatetime else 't (s)')
if usedatetime:
ax.yaxis_date()
else:
ax.set_xlim(t_vals[0], t_vals[-1])
ax.set_ylim(ch_vals[-1] + gap, ch_vals[0] - gap) # inverted (first ch at top)
ax.set_xlabel('time' if usedatetime else 't (s)')
ax.set_ylabel('channel')
if usedatetime:
ax.xaxis_date()
return ax
[docs]
def plot_xcorr(
cc,
lag_s=None,
fs: Optional[float] = None,
offsets=None,
ax=None,
perc: float = 99.0,
vmin: Optional[float] = None,
vmax: Optional[float] = None,
cmap: str = 'seismic',
cbar: bool = True,
figsize: Optional[Tuple[float, float]] = None,
**kwargs,
):
"""Image a cross-correlation gather — pair (x) vs lag time (y, downward).
Plots the ``(npair, nlag)`` array returned by `dasio.noise.xcorr_stack`
(the first tuple element) or `xcorr_dataset` (``result['cc']``): each column
is one channel pair's stacked correlation, zero lag at the centre row. Lag
time runs down the vertical axis (seismic-gather convention) and pairs along
the horizontal axis, so a common-shot gather shows the virtual-source moveout
fanning out symmetrically in +/- lag.
Parameters
----------
cc : ``(npair, nlag)`` array.
lag_s : lag axis in seconds (e.g. ``xcorr_dataset``'s ``result['lag_s']``).
If None it is built from `fs` (or falls back to sample indices), with
zero lag at the centre. Drawn on the vertical axis, increasing downward.
fs : sample rate, used to label the lag axis in seconds when `lag_s` is None.
offsets : optional per-pair x values (e.g. inter-channel offset in km for a
common-offset gather); defaults to the pair index.
perc, vmin, vmax, cmap, cbar, figsize : as in `imshow` — `perc` sets
symmetric colour limits from the `perc`-th percentile of ``|cc|`` unless
`vmin`/`vmax` are given.
**kwargs : forwarded to `ax.imshow`.
Returns ``(ax, im)``.
"""
import matplotlib.pyplot as plt
cc = np.asarray(cc)
npair, nlag = cc.shape
half = (nlag - 1) // 2
if lag_s is None:
lag_s = (np.arange(nlag) - half) / fs if fs else (np.arange(nlag) - half)
else:
lag_s = np.asarray(lag_s)
x = np.arange(npair) if offsets is None else np.asarray(offsets)
if vmin is None or vmax is None:
v = float(np.nanpercentile(np.abs(cc), perc))
if vmin is None:
vmin = -v
if vmax is None:
vmax = v
if ax is None:
_, ax = plt.subplots(figsize=figsize or (8, 10))
# lag on the vertical axis, increasing downward (top = first/most-negative lag)
extent = (float(x[0]), float(x[-1]), float(lag_s[-1]), float(lag_s[0]))
im = ax.imshow(cc.T, aspect='auto', cmap=cmap, vmin=vmin, vmax=vmax,
extent=extent, **kwargs)
ax.set_xlabel('offset' if offsets is not None else 'pair index')
ax.set_ylabel('lag (s)' if (lag_s is not None or fs) else 'lag (samples)')
if cbar:
ax.figure.colorbar(im, ax=ax, label='cross-correlation')
return ax, im
def _period_ticks(fmin, fmax, nmax=8):
"""Round periods (1-2-5 series, seconds) inside a frequency band."""
cands = [m * 10.0 ** e for e in range(-3, 7) for m in (1, 2, 5)]
ticks = [T for T in cands if 1.0 / fmax <= T <= 1.0 / fmin]
while len(ticks) > nmax:
ticks = ticks[::2]
return ticks
[docs]
def plot_spectrogram(
S,
freq=None,
times=None,
ich: int = 0,
t0=None,
usedatetime: Optional[bool] = None,
ax=None,
log_color: bool = True,
period_axis: Optional[bool] = None,
perc: float = 99.5,
vmin: Optional[float] = None,
vmax: Optional[float] = None,
cmap: str = 'turbo',
cbar: bool = True,
units: str = '',
figsize: Optional[Tuple[float, float]] = None,
**kwargs,
):
"""Image a time-resolved PSD from `psd.spectrogram` or `psd.spectrogram_cwt`.
Data first, like the other plotters here — either the whole result or the
array plus its axes:
plot_spectrogram(spectrogram_cwt(d.data, d.fs, 1/200, 0.5))
plot_spectrogram(S, freq, times, t0=d.start_time)
Parameters
----------
S : ``(nchan, nfreq, ntime)`` or ``(nfreq, ntime)`` array, or the whole
``(freq, times, S)`` tuple those functions return (accepted as-is
when `freq`/`times` are not given).
freq, times : the frequency (Hz) and time (s) axes; taken from the result
tuple, or bin/sample indices when absent.
ich : channel to draw when ``S`` is 3-D.
t0 : optional absolute start time — ``np.datetime64``, datetime, ISO
string, or unix seconds. Turns the x-axis into calendar time.
usedatetime : calendar-time x-axis, as in `imshow`. Default None means
"whenever `t0` is given". True without `t0` raises — ``times`` holds
seconds from the record start, so the epoch must come from somewhere.
False keeps the seconds axis even when `t0` is passed.
log_color : color by ``log10(S)`` (default) rather than linear PSD.
period_axis : log-frequency axis labeled in seconds of period. Default
None auto-enables it for a log-spaced grid (`spectrogram_cwt`) and
leaves a linear grid (`spectrogram`) on a linear frequency axis.
perc, vmin, vmax : color limits; unless given, the `perc`-th and
``(100-perc)``-th percentiles of the finite color values.
cmap, cbar, figsize : as in `imshow`.
units : physical unit of the data (e.g. ``'Pa'``) for the colorbar label,
shown as units^2/Hz.
**kwargs : forwarded to `ax.pcolormesh` (`rasterized` defaults to True —
thousands of quads otherwise swamp vector backends).
NaN cells — e.g. the cone of influence that `spectrogram_cwt` masks —
are left blank. Returns ``(ax, im)``.
"""
import matplotlib.pyplot as plt
if freq is None and times is None and isinstance(S, (tuple, list)) \
and len(S) == 3:
freq, times, S = S
S = np.asarray(S)
if S.ndim == 3:
S = S[ich]
if S.ndim != 2:
raise ValueError(f"S must be 2-D or 3-D, got shape {S.shape}")
freq = np.arange(S.shape[0]) if freq is None else np.asarray(freq, float)
times = np.arange(S.shape[1]) if times is None else np.asarray(times, float)
if freq.shape != (S.shape[0],) or times.shape != (S.shape[1],):
raise ValueError(
f"axis lengths {freq.shape}, {times.shape} do not match S {S.shape}"
)
if period_axis is None:
r = np.diff(np.log(freq)) if (freq > 0).all() and freq.size > 2 else None
period_axis = r is not None and np.allclose(r, r[0], rtol=1e-3)
with np.errstate(divide='ignore', invalid='ignore'):
C = np.log10(np.where(S > 0, S, np.nan)) if log_color else S
finite = C[np.isfinite(C)]
if finite.size == 0:
raise ValueError("spectrogram holds no positive finite value")
if vmax is None:
vmax = float(np.percentile(finite, perc))
if vmin is None:
vmin = float(np.percentile(finite, 100.0 - perc))
if usedatetime is None:
usedatetime = t0 is not None
if usedatetime and t0 is None:
raise ValueError(
"usedatetime=True needs t0 — `times` holds seconds from the "
"record start, so the absolute epoch must be supplied"
)
if usedatetime:
import matplotlib.dates as mdates
if isinstance(t0, (int, float)):
t0 = np.datetime64(int(round(t0)), 's')
x = mdates.date2num(np.datetime64(t0)) + times / 86400.0
else:
x = times
if ax is None:
_, ax = plt.subplots(figsize=figsize or (9, 4))
kwargs.setdefault('rasterized', True)
kwargs.setdefault('shading', 'auto')
im = ax.pcolormesh(x, freq, C, cmap=cmap, vmin=vmin, vmax=vmax, **kwargs)
if period_axis:
fpos = freq[freq > 0] # a linear grid's DC row has no period
ax.set_yscale('log')
ax.set_ylim(fpos.min(), fpos.max())
ticks = _period_ticks(fpos.min(), fpos.max())
if ticks:
ax.set_yticks([1.0 / T for T in ticks])
ax.set_yticklabels([f'{T:g}' for T in ticks])
ax.minorticks_off()
ax.set_ylabel('Period (s)')
else:
ax.set_ylabel('f (Hz)')
if usedatetime:
import matplotlib.dates as mdates
loc = mdates.AutoDateLocator()
ax.xaxis.set_major_locator(loc)
ax.xaxis.set_major_formatter(mdates.ConciseDateFormatter(loc))
ax.set_xlabel('time')
else:
ax.set_xlabel('t (s)')
if cbar:
ub = f'{units}$^2$/Hz' if units else 'data$^2$/Hz'
label = f'log$_{{10}}$ PSD ({ub})' if log_color else f'PSD ({ub})'
ax.figure.colorbar(im, ax=ax, label=label)
return ax, im
def _db_label(db, scaling):
"""Colour-bar label matching the scaling `_fk_amplitude` applied."""
if not db:
return 'amplitude'
return 'dB re max' if scaling is None else 'dB'
def _fk_amplitude(A, db, perc, vmin, vmax, scaling=None):
"""Amplitude in dB (or linear), plus the percentile colour limits.
Args:
A: complex or real panel.
db: convert to dB.
perc: percentile for `vmax` when it is not given.
vmin: lower colour limit, or None to derive it.
vmax: upper colour limit, or None to derive it.
scaling: what `A` is, which fixes both the dB law and the reference.
None treats it as a raw, uncalibrated amplitude and scales to the
panel maximum. 'density' treats it as the calibrated power spectral
density `fk(scaling='density')` returns — 10*log10, since it is a
power, against one unit rather than the panel peak.
Returns:
A: magnitude, in dB if requested.
vmin: resolved lower colour limit.
vmax: resolved upper colour limit.
Raises:
ValueError: if the panel holds no positive finite value, or `scaling`
is neither None nor 'density'.
The default normalises to the panel's own maximum, because a raw f-k
amplitude carries whatever scaling the gather and the transform size gave
it — only the contrast means anything, and `vmax` sits at 0 by
construction. The consequence is that two panels cannot be compared with
each other: each is measured against its own peak. Pass the `scaling`
that `fk` produced the panel with to keep absolute, comparable numbers.
"""
A = np.abs(np.asarray(A))
if db:
if scaling is None:
factor, denom = 20.0, (A.max() or 1.0)
elif scaling == 'density':
factor, denom = 10.0, 1.0 # a power, so 10*log10, absolute
else:
raise ValueError(
f"scaling must be None or 'density', got {scaling!r}")
with np.errstate(divide='ignore'):
A = factor * np.log10(np.where(A > 0, A, np.nan) / denom)
finite = A[np.isfinite(A)]
if finite.size == 0:
raise ValueError('the panel holds no positive finite value')
if vmax is None:
vmax = float(np.percentile(finite, perc))
if vmin is None:
vmin = vmax - 60.0 if db else float(np.percentile(finite, 100.0 - perc))
return A, vmin, vmax
[docs]
def plot_fk(F, f=None, k=None, ax=None, db: bool = True, perc: float = 100.0,
vmin: Optional[float] = None, vmax: Optional[float] = None,
scaling: Optional[str] = None, cmap: str = 'turbo',
cbar: bool = True,
velocities=None,
figsize: Optional[Tuple[float, float]] = None, **kwargs):
"""Image the f-k plane from `dasio.fk.fk`.
Args:
F: f-k amplitude from `fk`, or the whole `(F, f, k)` result.
f: frequency axis [Hz]; taken from `F` when it is the whole result.
k: wavenumber axis [rad/m]; likewise.
ax: axes to draw on; a new figure is made when None.
db: scale to dB re the panel maximum.
perc: percentile for the upper colour limit.
vmin: lower colour limit; defaults to 60 dB below `vmax`.
vmax: upper colour limit.
scaling: what the panel is, matching how `fk` produced it — the same
token in both places. None (default) is a raw complex amplitude:
dB is taken relative to the panel maximum, because the absolute
level of a raw transform has no meaning (it moves with the sample
count), and the consequence is that two panels cannot be compared
with each other. 'density' is what `fk(..., scaling='density')`
returns: a power, so 10*log10 against one unit, absolute and
comparable between panels.
cmap: colormap name.
cbar: draw the colour bar.
velocities: `c = 2*pi*f/k` guide lines [m/s] — positive values draw
both directions, signed values draw one.
figsize: figure size when `ax` is None.
**kwargs: forwarded to `pcolormesh`.
Returns:
ax: the axes drawn on.
im: the `QuadMesh`, for a colour bar or limit tweak.
Raises:
ValueError: if `F` is not 2-D or the axes do not match it.
Data first, like the other plotters here — either the whole result or the
array plus its axes:
plot_fk(*fk(d.data, d.dx, d.dt))
plot_fk(F, f, k, velocities=[1500, 3000])
Frequency runs along x and wavenumber along y, so the apparent-velocity
fans open sideways and a dispersive arrival reads left to right. `db`
scales relative to the panel maximum, since an f-k amplitude inherits the
gather's units and only its contrast means anything.
"""
import matplotlib.pyplot as plt
if f is None and k is None and isinstance(F, (tuple, list)) and len(F) == 3:
F, f, k = F
F = np.asarray(F)
if F.ndim != 2:
raise ValueError(f'F must be 2-D, got shape {F.shape}')
k = np.arange(F.shape[0]) if k is None else np.asarray(k, float)
f = np.arange(F.shape[1]) if f is None else np.asarray(f, float)
if k.shape != (F.shape[0],) or f.shape != (F.shape[1],):
raise ValueError(
f'axis lengths {k.shape}, {f.shape} do not match F {F.shape}')
A, vmin, vmax = _fk_amplitude(F, db, perc, vmin, vmax, scaling)
if ax is None:
_, ax = plt.subplots(figsize=figsize or (7, 5))
kwargs.setdefault('rasterized', True)
kwargs.setdefault('shading', 'auto')
im = ax.pcolormesh(f, k, A, cmap=cmap, vmin=vmin, vmax=vmax, **kwargs)
vs = np.atleast_1d(velocities if velocities is not None else [])
if vs.size and np.all(vs > 0):
vs = np.concatenate([vs, -vs]) # unsigned means "both fans"
for c in vs:
# k = 2*pi*f/c is a line through the origin; drawing it against the
# frequency axis keeps it right whatever the k limits are.
ax.plot(f, 2.0 * np.pi * f / c, color='w', lw=0.8, ls='--', alpha=0.7)
ax.set(
xlabel='f (Hz)', ylabel='k (rad/m)',
xlim=(f.min(), f.max()), ylim=(k.min(), k.max()),
)
if cbar:
ax.figure.colorbar(im, ax=ax, label=_db_label(db, scaling))
return ax, im
[docs]
def plot_fc(A, f=None, c=None, ax=None, db: bool = False, perc: float = 100.0,
vmin: Optional[float] = None, vmax: Optional[float] = None,
scaling: Optional[str] = None,
cmap: str = 'turbo', cbar: bool = True, picks=None,
figsize: Optional[Tuple[float, float]] = None, **kwargs):
"""Image the dispersion panel from `dasio.fk.fk2fc`.
Args:
A: dispersion panel from `fk2fc`, or the whole `(A, f, c)` result.
f: frequency axis [Hz]; taken from `A` when it is the whole result.
c: velocity axis [m/s]; likewise.
ax: axes to draw on; a new figure is made when None.
db: scale to dB re the panel maximum.
perc: percentile for the upper colour limit.
vmin: lower colour limit.
vmax: upper colour limit.
scaling: what the panel is, as in `plot_fk`; pass 'density' for one
built from `fk(..., scaling='density')`.
cmap: colormap name.
cbar: draw the colour bar.
picks: `(f, c)` sequences to overlay, e.g. a picked dispersion curve.
figsize: figure size when `ax` is None.
**kwargs: forwarded to `pcolormesh`.
Returns:
ax: the axes drawn on.
im: the `QuadMesh`, for a colour bar or limit tweak.
Raises:
ValueError: if `A` is not 2-D or the axes do not match it.
plot_fc(*fk2fc(F, f, k, 200, 5000))
Frequency on x, apparent velocity on y — the orientation dispersion curves
are picked in. Linear colour by default: a dispersion image is read by
where its ridge is, and dB flattens the contrast that makes it stand out.
"""
import matplotlib.pyplot as plt
if f is None and c is None and isinstance(A, (tuple, list)) and len(A) == 3:
A, f, c = A
A = np.asarray(A)
if A.ndim != 2:
raise ValueError(f'A must be 2-D, got shape {A.shape}')
c = np.arange(A.shape[0]) if c is None else np.asarray(c, float)
f = np.arange(A.shape[1]) if f is None else np.asarray(f, float)
if c.shape != (A.shape[0],) or f.shape != (A.shape[1],):
raise ValueError(
f'axis lengths {c.shape}, {f.shape} do not match A {A.shape}')
A, vmin, vmax = _fk_amplitude(A, db, perc, vmin, vmax, scaling)
if ax is None:
_, ax = plt.subplots(figsize=figsize or (7, 5))
kwargs.setdefault('rasterized', True)
kwargs.setdefault('shading', 'auto')
im = ax.pcolormesh(f, c, A, cmap=cmap, vmin=vmin, vmax=vmax, **kwargs)
if picks is not None:
pf, pc = np.asarray(picks[0]), np.asarray(picks[1])
ax.plot(pf, pc, 'w.-', lw=1.0, ms=4)
ax.set(
xlabel='f (Hz)', ylabel='c (m/s)',
xlim=(f.min(), f.max()), ylim=(c.min(), c.max()),
)
if cbar:
ax.figure.colorbar(im, ax=ax, label=_db_label(db, scaling))
return ax, im
[docs]
def plot_timeline(db, ax=None, by_month=None):
"""Gantt-style coverage timeline for a `DASdb`.
Green where data exists, red for the acquisition gaps between.
With `by_month=True` each calendar month becomes its own row against a
shared day-of-month axis, which keeps a multi-month archive legible —
a year on one continuous axis compresses every gap into a hairline.
`None` (default) picks that layout automatically when the catalog spans
more than one month.
Bar geometry is computed in matplotlib date numbers rather than handed
over as Timestamp/Timedelta objects. Mixing the two lets the unit
machinery interpret a width in the same frame as a position, and a
single bad conversion is enough to push the axes limits far enough that
rendering a `bbox_inches='tight'` figure asks for a canvas of millions
of pixels.
"""
try:
import matplotlib.pyplot as plt
from matplotlib import dates as mdates
except ImportError as e:
raise ImportError(
'plot_timeline needs matplotlib; already a core dep for '
'the monitor but was not found in this env'
) from e
# Only this plotter needs pandas; `imshow` and `wiggle` are numpy-only and
# this module is imported for them far more often.
import pandas as pd
segments = list(db.segments())
if not segments:
raise RuntimeError('DASdb is empty; nothing to plot')
spans = [(s['begin_time'].iloc[0], s['end_time'].iloc[-1])
for s in segments]
gaps = [(a[1], b[0]) for a, b in zip(spans[:-1], spans[1:]) if b[0] > a[1]]
months = pd.period_range(spans[0][0], spans[-1][1], freq='M')
if by_month is None:
by_month = len(months) > 1
if not by_month:
if ax is None:
_, ax = plt.subplots(figsize=(12, 2.0), layout='constrained')
for t0, t1 in spans:
x0, x1 = mdates.date2num(t0), mdates.date2num(t1)
ax.barh(0, x1 - x0, left=x0, height=0.6,
color='#66bb66', edgecolor='#2e7d32', alpha=0.9)
for t0, t1 in gaps:
x0, x1 = mdates.date2num(t0), mdates.date2num(t1)
ax.barh(0, x1 - x0, left=x0, height=0.6, color='#d55', alpha=0.4)
ax.set_xlim(mdates.date2num(spans[0][0]), mdates.date2num(spans[-1][1]))
ax.set_ylim(-0.5, 0.5)
ax.set_yticks([])
ax.xaxis.set_major_formatter(
mdates.ConciseDateFormatter(ax.xaxis.get_major_locator()))
ax.set_xlabel('time')
else:
if ax is None:
_, ax = plt.subplots(figsize=(12, 0.42 * len(months) + 1.1),
layout='constrained')
for row, month in enumerate(months):
m0 = month.start_time.tz_localize(spans[0][0].tz)
m1 = (month + 1).start_time.tz_localize(spans[0][0].tz)
span_days = (m1 - m0).total_seconds() / 86400.0
ax.barh(row, span_days, left=0, height=0.62,
color='#f2f2f2', edgecolor='#dddddd', zorder=0)
for (t0, t1), color, alpha in (
[(s, '#66bb66', 0.95) for s in spans]
+ [(g, '#d55', 0.55) for g in gaps]):
lo, hi = max(t0, m0), min(t1, m1)
if hi <= lo:
continue # this bar does not reach into this month
left = (lo - m0).total_seconds() / 86400.0
ax.barh(row, (hi - lo).total_seconds() / 86400.0, left=left,
height=0.62, color=color, alpha=alpha, zorder=2)
ax.set_ylim(len(months) - 0.5, -0.5) # first month at the top
ax.set_yticks(range(len(months)))
ax.set_yticklabels([str(m) for m in months], fontsize=9)
ax.set_xlim(0, 31)
ax.set_xticks(range(0, 31, 5))
ax.set_xticklabels([str(d + 1) for d in range(0, 31, 5)])
ax.set_xlabel('day of month')
ax.grid(axis='x', color='white', lw=0.8, zorder=1)
ax.set_title(
f'{db.format} | {db.n_files:,} files | '
f'{db.n_segments} segments | '
f'{spans[0][0]:%Y-%m-%d} to {spans[-1][1]:%Y-%m-%d}',
fontsize=10,
)
return ax
class _PlotAccessor:
"""Bound plot interface for a single DASdata.
Reachable as `d.plot`. `d.plot()` is shorthand for `d.plot.imshow()`;
`d.plot.imshow(...)` and `d.plot.wiggle(...)` reach the named
plotters. `pcolormesh` will be added when a non-uniform channel
axis use case shows up (truncate + uniform `dx` covers the rest).
"""
__slots__ = ('_d',)
def __init__(self, d):
self._d = d
def __call__(self, **kw):
return self.imshow(**kw)
def imshow(self, **kw):
return imshow(self._d, **kw)
def wiggle(self, **kw):
return wiggle(self._d, **kw)