[FIX] Report plot polish — colorbar ticks, cyclic direction, T-S bounds, tilt panels

oceanarray/plotters/helpers.py

@@ -176,7 +176,7 @@ def _velocity_panel_style(
     div (diverging)  east/north/up/error_velocity  Spectral_r  ± div_abs_max
     seq (sequential) current_speed                 plasma       0 → 98th pctile
     seq              bin_pressure                  PuRd         2nd→98th pctile
-    cyc (cyclic)     current_direction             hsv          0–360°
+    cyc (cyclic)     current_direction             twilight     0–360°
 
     Parameters
     ----------
@@ -195,7 +195,9 @@ def _velocity_panel_style(
     tuple
         ``(bounds, norm, cmap, cb_label)`` where *bounds* is a 1-D array of
         colorbar tick/boundary values, *norm* is a ``BoundaryNorm``, *cmap*
-        is a colormap name string, and *cb_label* is the colorbar axis label.
+        is a colormap name string or a :class:`~matplotlib.colors.Colormap`
+        (cyclic panels return a resampled colormap), and *cb_label* is the
+        colorbar axis label.
 
     """
     import matplotlib.colors as mcolors
@@ -210,11 +212,17 @@ def _velocity_panel_style(
         bounds, norm = colorbar_norm(vmin=0.0, vmax=max(spd_max, 1e-4))
         return bounds, norm, "plasma", "m s⁻¹"
     if var == "current_direction":
-        bounds = np.linspace(0, 360, 21)
-        norm = mcolors.BoundaryNorm(bounds, ncolors=256)
+        import matplotlib.pyplot as plt
+
         # Cyclic colormap so 0° and 360° share a colour; twilight is perceptually
-        # uniform (hsv is not).  Hard-coded here, not in parameters.py.
-        return bounds, norm, "twilight", "°"
+        # uniform (hsv is not).  twilight has 510 native colours, so a plain
+        # BoundaryNorm(ncolors=256) maps 0–360° onto only the first half of the
+        # cycle (pale→dark, non-cyclic — Eleanor 2026-08-18).  Resample to one
+        # colour per bin and match ncolors so the full cycle shows and 0°==360°.
+        bounds = np.linspace(0, 360, 21)
+        cmap = plt.get_cmap("twilight", len(bounds) - 1)
+        norm = mcolors.BoundaryNorm(bounds, cmap.N)
+        return bounds, norm, cmap, "°"
     if var == "bin_pressure":
         p_lo = float(np.percentile(finite_vals, 2)) if len(finite_vals) else 0.0
         p_hi = float(np.percentile(finite_vals, 98)) if len(finite_vals) else 1000.0

oceanarray/plotters/timeseries.py

@@ -457,16 +457,25 @@ def draw_grid_n2(
     CT = gsw.CT_from_t(SA, T_pt, p_2d)
     N2, p_mid = gsw.Nsquared(SA, CT, p_2d, lat=lat)
     p_mid_1d = np.nanmean(p_mid, axis=1)
+    # Floor non-positive/unstable N² to 1e-12 so it renders as the lowest colour
+    # (weakly stratified / off-scale low) rather than blank.
     N2_log = np.log10(np.maximum(N2, 1e-12))
 
     fig, ax = plt.subplots(
         figsize=(width_in, params.GRID_PANEL_ROW_IN), layout="constrained"
     )
-    # Clip the colorbar to the 2.5-97.5 percentiles of log10(N²) so a few extreme
-    # cells don't wash out the stratification structure.
-    _finite = N2_log[np.isfinite(N2_log)]
-    _lo = float(np.nanpercentile(_finite, 2.5)) if _finite.size else -12.0
-    _hi = float(np.nanpercentile(_finite, 97.5)) if _finite.size else 0.0
+    # Colour limits from the 1st-99th percentile of log10(N²) over *positive*
+    # (statically stable) cells only.  Non-positive N² is floored to 1e-12 for the
+    # plot, but including that floor in the percentile would peg the low limit at
+    # log10(1e-12) = -12 and wash out the stratification structure whenever more
+    # than ~1% of cells are unstable (Eleanor 2026-08-18).
+    _pos = N2[np.isfinite(N2) & (N2 > 0)]
+    if _pos.size:
+        _logpos = np.log10(_pos)
+        _lo = float(np.nanpercentile(_logpos, 1))
+        _hi = float(np.nanpercentile(_logpos, 99))
+    else:
+        _lo, _hi = -12.0, 0.0
     bounds, norm = colorbar_norm(vmin=_lo, vmax=_hi)
     pc = ax.pcolormesh(
         time_vals, p_mid_1d, N2_log, shading="nearest", cmap="plasma_r", norm=norm

oceanarray/plotters/ts.py

@@ -283,6 +283,44 @@ def draw_ts_diagram(
     return fig
 
 
+def _nice_axis_limits(
+    x: np.ndarray, *, plow: float = 1.0, phigh: float = 99.0, pad_frac: float = 0.05
+) -> "tuple[float, float]":
+    """Return padded, outward-rounded axis limits from robust percentiles.
+
+    Takes the *plow*/*phigh* percentiles of *x* (default 1st/99th) so outliers do
+    not stretch the axis, pads by *pad_frac* of that percentile range on each side
+    (default 5%, i.e. 10% of the range total), then rounds the low bound down and
+    the high bound up to a clean step (one tenth of the range's order of
+    magnitude) so the axis ticks land on round values.
+
+    Parameters
+    ----------
+    x : numpy.ndarray
+        Finite data values for the axis.
+    plow, phigh : float, optional
+        Low/high percentiles bounding the data core.  Default 1 and 99.
+    pad_frac : float, optional
+        Fraction of the percentile range added as padding on each side.
+        Default 0.05.
+
+    Returns
+    -------
+    tuple of float
+        ``(lo, hi)`` axis limits.
+
+    """
+    p_lo = float(np.nanpercentile(x, plow))
+    p_hi = float(np.nanpercentile(x, phigh))
+    rng = p_hi - p_lo
+    if not np.isfinite(rng) or rng <= 0:
+        return p_lo - 0.5, p_hi + 0.5
+    lo = p_lo - pad_frac * rng
+    hi = p_hi + pad_frac * rng
+    step = 10.0 ** np.floor(np.log10(rng)) / 10.0
+    return float(np.floor(lo / step) * step), float(np.ceil(hi / step) * step)
+
+
 def draw_stack_ts_diagram(
     ds: "xr.Dataset", *, width_in: float = report_tokens.W_FULL
 ) -> "Optional[plt.Figure]":
@@ -392,13 +430,11 @@ def draw_stack_ts_diagram(
     ax_scatter.set_ylabel(params.vlabel("temperature"))
     ax_scatter.set_title("T-S (colour = pressure)")
 
-    # Shared limits (scatter's data range) so the heatmap matches (boxes are
-    # already square from square_axes_grid).
-    _sf, _tf = S_flat[finite], T_flat[finite]
-    _sp = 0.02 * (float(np.nanmax(_sf) - np.nanmin(_sf)) + 1e-9)
-    _tp = 0.02 * (float(np.nanmax(_tf) - np.nanmin(_tf)) + 1e-9)
-    _s_lim = (float(np.nanmin(_sf)) - _sp, float(np.nanmax(_sf)) + _sp)
-    _t_lim = (float(np.nanmin(_tf)) - _tp, float(np.nanmax(_tf)) + _tp)
+    # Shared limits from robust percentiles (1/99, padded 5% each side and rounded
+    # outward) so a few outliers don't stretch the box; the heatmap matches (boxes
+    # are already square from square_axes_grid).
+    _s_lim = _nice_axis_limits(S_flat[finite])
+    _t_lim = _nice_axis_limits(T_flat[finite])
     ax_scatter.set_xlim(*_s_lim)
     ax_scatter.set_ylim(*_t_lim)
 
@@ -496,11 +532,11 @@ def draw_grid_ts_diagram(
     if finite.sum() < 10:
         return None
 
-    # Axis limits from the data — returned to caller so hydro panels share scales.
-    s_lo = float(np.nanpercentile(S[finite], 0.01))
-    s_hi = float(np.nanpercentile(S[finite], 99.99))
-    t_lo = float(np.nanpercentile(T[finite], 0.01))
-    t_hi = float(np.nanpercentile(T[finite], 99.99))
+    # Axis limits from robust percentiles (1/99, padded and rounded outward) so
+    # outliers don't stretch the box — returned to caller so hydro panels share
+    # the same scales.
+    s_lo, s_hi = _nice_axis_limits(S[finite])
+    t_lo, t_hi = _nice_axis_limits(T[finite])
     ts_bounds: dict = {"t_lim": (t_lo, t_hi), "s_lim": (s_lo, s_hi), "o2_lim": None}
 
     has_o2 = "oxygen_saturation_pct" in ds.data_vars

oceanarray/reports/_stack.py

@@ -56,9 +56,11 @@
 def _make_aquadopp_tilt_panels(ds: Any, step: int = 1) -> Optional[str]:
     """One subplot per Aquadopp showing pitch, roll, and tilt_from_pressure.
 
-    All three curves share the same y-axis so they can be compared directly.
-    Horizontal reference lines are drawn at the suspect and fail thresholds
-    read from ds.attrs (falling back to 20° / 30° if absent).
+    Within each row the time-series panel and its paired scatter panel share one
+    0–90° y-axis, so the pitch/roll time series and the pressure-tilt scatter read
+    on the same scale.  Gridlines are drawn on both panels.  Horizontal reference
+    lines are drawn at the suspect and fail thresholds read from ds.attrs (falling
+    back to 20° / 30° if absent).
     Returns None if no Aquadopp levels are found or none of the relevant
     variables exist.
     """
@@ -87,7 +89,7 @@ def _make_aquadopp_tilt_panels(ds: Any, step: int = 1) -> Optional[str]:
         return None
 
     # Cap the number of stacked panels so the figure stays a sane height for PDF
-    # pagination (one panel is ~2.8 in; 5 → ~14 in).  Deep-first order is kept, so
+    # pagination (one panel is ~2.24 in; 5 → ~11 in).  Deep-first order is kept, so
     # the deepest Aquadopps are shown; a note flags any that were dropped.
     _MAX_TILT_ROWS = 5
     _n_dropped = max(0, len(aq_indices) - _MAX_TILT_ROWS)
@@ -104,8 +106,11 @@ def _make_aquadopp_tilt_panels(ds: Any, step: int = 1) -> Optional[str]:
 
     def _draw() -> "plt.Figure":
         fig = plt.figure(
-            figsize=(report_tokens.W_FULL, 2.8 * n_panels), constrained_layout=True
+            figsize=(report_tokens.W_FULL, 2.24 * n_panels), constrained_layout=True
         )
+        # Tighten the inter-row gap; upper rows drop their x tick labels/label so
+        # the rows can pack closer.
+        fig.get_layout_engine().set(h_pad=0.03, hspace=0.0)
         gs = fig.add_gridspec(n_panels, 3, width_ratios=[2, 2, 1])
 
         ax_ts_first = None
@@ -116,7 +121,9 @@ def _draw() -> "plt.Figure":
             ax_ts = fig.add_subplot(gs[row, :2], sharex=ax_ts_first)
             if ax_ts_first is None:
                 ax_ts_first = ax_ts
-            ax_sc = fig.add_subplot(gs[row, 2])
+            # Scatter shares this row's y-axis so pitch/roll read on the same
+            # 0–90° scale as the time series.
+            ax_sc = fig.add_subplot(gs[row, 2], sharey=ax_ts)
 
             p_data = r_data = tp_data = None
             if has_pitch:
@@ -143,8 +150,9 @@ def _draw() -> "plt.Figure":
 
             ax_ts.axhline(tilt_suspect, color="tab:orange", lw=0.8, ls="--", zorder=0)
             ax_ts.axhline(tilt_fail, color="tab:red", lw=0.8, ls=":", zorder=0)
-            ax_ts.set_ylim(bottom=0.0)
+            ax_ts.set_ylim(0.0, 90.0)
             ax_ts.set_ylabel("Degrees (°)")
+            grid_despine(ax_ts)
 
             _ref_note = ""
             if ref_habs is not None and np.isfinite(ref_habs[i]):
@@ -185,10 +193,9 @@ def _draw() -> "plt.Figure":
                         label="|roll|",
                         **sc_kw,
                     )
-                _lim = max(ax_sc.get_xlim()[1], ax_sc.get_ylim()[1], 35.0)
                 ax_sc.plot(
-                    [0, _lim],
-                    [0, _lim],
+                    [0, 90],
+                    [0, 90],
                     color="0.4",
                     lw=0.8,
                     ls="--",
@@ -203,10 +210,18 @@ def _draw() -> "plt.Figure":
                     tilt_suspect, color="tab:orange", lw=0.7, ls="--", zorder=0
                 )
                 ax_sc.axhline(tilt_fail, color="tab:red", lw=0.7, ls=":", zorder=0)
-                ax_sc.set_xlim(left=0.0)
-                ax_sc.set_ylim(bottom=0.0)
-                ax_sc.set_xlabel("tilt (pressure) [°]")
+                ax_sc.set_xlim(0.0, 90.0)
+                ax_sc.set_ylim(0.0, 90.0)  # re-assert after scatter autoscale
+                # Match x/y ticks (both are 0–90° tilt); the shared y also sets the
+                # time-series panel's y-ticks.
+                ax_sc.set_xticks([0, 30, 60, 90])
+                ax_sc.set_yticks([0, 30, 60, 90])
                 ax_sc.set_ylabel("|pitch|, |roll| [°]")
+                if row == n_panels - 1:
+                    ax_sc.set_xlabel("tilt (pressure) [°]")
+                else:
+                    ax_sc.tick_params(labelbottom=False)
+                grid_despine(ax_sc)
                 # No scatter legend — the time-series panel legend already names
                 # pitch/roll; a second legend here is a redundant repeat.
             else:

oceanarray/utilities.py

@@ -449,31 +449,31 @@ def _safe_rel(path: Path, root: Path) -> str:
 def _nice_colorbar_bounds(vmin: float, vmax: float, n: int = 20) -> np.ndarray:
     """Return a boundary array for a discrete colorbar with approximately *n* levels.
 
-    The step is rounded to 1 significant figure so tick labels land on clean values.
-    The range is centered on the midpoint of [vmin, vmax].
+    The step is snapped to the "nice" 1 / 2 / 2.5 / 5 family (the same set
+    :class:`matplotlib.ticker.MaxNLocator` uses), and the boundaries are aligned
+    to integer multiples of that step spanning ``[vmin, vmax]``.  This guarantees
+    that round labelled ticks (0.2, 100, …) land exactly on colour-step
+    boundaries instead of between them — otherwise a boundary step of, say, 0.03
+    or 30 leaves the labels sitting between the discrete colour changes.
 
     Examples
     --------
-    Temperature  vmin=0.87, vmax=7.23 → step=0.3,  bounds 1.2–7.2  (20 levels)
-    Salinity     vmin=34.782, vmax=35.139 → step=0.02, bounds 34.76–35.16 (20 levels)
+    Temperature  vmin=0.87, vmax=7.23 → step=0.25, bounds 0.75–7.25
+    Salinity     vmin=34.782, vmax=35.139 → step=0.02, bounds 34.78–35.14
 
     """
     span = vmax - vmin
     if span <= 0:
         return np.linspace(vmin - 1, vmin + 1, n + 1)
     raw_step = span / n
     mag = 10.0 ** math.floor(math.log10(raw_step))
-    rounded = round(raw_step / mag)
-    if rounded == 0:
-        rounded = 1
-    elif rounded >= 10:
-        mag *= 10
-        rounded = 1
-    nice_step = rounded * mag
-    mid = (vmin + vmax) / 2
-    mid_aligned = round(mid / nice_step) * nice_step
-    lo = mid_aligned - (n / 2) * nice_step
-    return np.array([lo + i * nice_step for i in range(n + 1)])
+    frac = raw_step / mag  # in [1, 10)
+    nice_frac = min((1.0, 2.0, 2.5, 5.0, 10.0), key=lambda c: abs(c - frac))
+    nice_step = nice_frac * mag
+    lo = math.floor(vmin / nice_step) * nice_step
+    hi = math.ceil(vmax / nice_step) * nice_step
+    n_steps = int(round((hi - lo) / nice_step))
+    return lo + nice_step * np.arange(n_steps + 1)
 
 
 def nice_colorbar_ticks(vmin: float, vmax: float, *, max_ticks: int = 6) -> np.ndarray:

tests/unit/test_plotters.py

@@ -165,3 +165,30 @@ def test_plot_trajectory_lc_array_length():
         f"expected {len(x) - 1} (one per segment)"
     )
     plt.close(fig)
+
+
+def test_nice_axis_limits_pads_and_rounds():
+    """1/99 percentiles, +5% pad each side, rounded outward to a clean step."""
+    from oceanarray.plotters.ts import _nice_axis_limits
+
+    # linspace 0..100 has exact 1st/99th percentiles of 1 and 99 (range 98):
+    # 1 - 0.05*98 = -3.9 floors to -4; 99 + 4.9 = 103.9 ceils to 104 (step 1).
+    x = np.linspace(0.0, 100.0, 101)
+    assert _nice_axis_limits(x) == (-4.0, 104.0)
+
+
+def test_nice_axis_limits_excludes_outliers():
+    """Outliers past the 1/99 percentiles must not stretch the limits."""
+    from oceanarray.plotters.ts import _nice_axis_limits
+
+    x = np.concatenate([np.linspace(34.66, 35.18, 999), [30.0, 40.0]])
+    lo, hi = _nice_axis_limits(x)
+    assert 34.5 < lo < 34.66 and 35.18 < hi < 35.4  # outliers 30/40 excluded
+
+
+def test_nice_axis_limits_degenerate_range():
+    """Zero-width data falls back to a symmetric ±0.5 window (no crash)."""
+    from oceanarray.plotters.ts import _nice_axis_limits
+
+    lo, hi = _nice_axis_limits(np.full(50, 5.0))
+    assert lo == 4.5 and hi == 5.5

tests/unit/test_utilities.py

@@ -264,10 +264,31 @@ def test_get_dims_2d_pressure():
 
 def test_nice_colorbar_bounds_standard():
     bounds = utilities._nice_colorbar_bounds(0.5, 7.5, n=20)
-    assert len(bounds) == 21  # n+1 edges
     assert np.all(np.diff(bounds) > 0), "bounds must be monotonically increasing"
-    # All values should be finite
     assert np.all(np.isfinite(bounds))
+    # Boundaries span the data range.
+    assert bounds[0] <= 0.5 and bounds[-1] >= 7.5
+    # Step snapped to the nice 1/2/2.5/5 family, and roughly n levels.
+    step = bounds[1] - bounds[0]
+    frac = step / 10.0 ** np.floor(np.log10(step))
+    assert np.isclose(frac, min((1.0, 2.0, 2.5, 5.0), key=lambda c: abs(c - frac)))
+    assert 10 <= len(bounds) <= 40  # count is adaptive, no longer fixed at n+1
+
+
+def test_nice_colorbar_ticks_align_to_bounds():
+    """Labelled ticks must land on discrete colour-step boundaries.
+
+    Regression guard: the boundary step used to be rounded to any 1-sig-fig value
+    (0.03, 30), so round ticks (0.2, 100) fell between the colour changes rather
+    than on them (Eleanor 2026-08-18).
+    """
+    for vmin, vmax in [(27.5, 27.95), (82.0, 808.0), (0.5, 7.5)]:
+        bounds = utilities._nice_colorbar_bounds(vmin, vmax)
+        ticks = utilities.nice_colorbar_ticks(float(bounds[0]), float(bounds[-1]))
+        for t in ticks:
+            assert np.any(np.isclose(t, bounds)), (
+                f"tick {t} off the boundary grid for ({vmin}, {vmax})"
+            )
 
 
 def test_nice_colorbar_bounds_zero_span():
@@ -278,11 +299,10 @@ def test_nice_colorbar_bounds_zero_span():
 
 
 def test_nice_colorbar_bounds_salinity_range():
-    """Typical salinity range gives clean step and correct count."""
+    """Typical salinity range gives a clean 0.02 step covering the data."""
     bounds = utilities._nice_colorbar_bounds(34.78, 35.14, n=20)
-    assert len(bounds) == 21
-    assert bounds[0] < 34.78 or np.isclose(bounds[0], 34.78, atol=0.1)
-    assert bounds[-1] > 35.14 or np.isclose(bounds[-1], 35.14, atol=0.1)
+    assert bounds[0] <= 34.78 and bounds[-1] >= 35.14
+    assert np.isclose(bounds[1] - bounds[0], 0.02)
 
 
 def _touch(path, mtime):