def low_shear_regions(rho, shear, threshold, *, psi_n=None, q=None, pressure=None) -> tuple[LowShearRegion, ...]:
    """Connected radial intervals where ``|s| < threshold``.

    Parameters
    ----------
    rho : array_like
        Radial coordinate, increasing; the one the shear was taken in [-].
    shear : array_like
        Magnetic shear on ``rho`` [-].
    threshold : float
        The ``|s|`` bound, always explicit: no value of it is a stability boundary [-].
    psi_n : array_like, optional
        Normalized poloidal flux on the same grid, for the ``psi_n`` edges [-].
    q : array_like, optional
        Safety factor on the same grid, for the q range and q=1 containment [-].
    pressure : array_like, optional
        Pressure on the same grid, for the drive-context descriptors [Pa].

    Returns
    -------
    regions : tuple of LowShearRegion
        Every interval, innermost first, at sample resolution [-].

    Applicability
    -------------
    Machine-independent.

    Limitations
    -----------
    Intervals end at the last sample inside, without interpolating the
    threshold crossing, so a width is resolved to the grid spacing. An interval
    is closed only by a sample outside it: a non-finite shear sample breaks it.
    ``contains_q_min`` compares against the global ``|q|`` minimum on the grid.
    """
    if not threshold > 0:
        raise ValueError(f"threshold must be positive, not {threshold!r}")
    rho = np.asarray(rho, dtype=float).ravel()
    shear = np.asarray(shear, dtype=float).ravel()
    if rho.shape != shear.shape:
        raise ValueError(f"rho and shear differ in length: {rho.size} vs {shear.size}")
    psi = rho if psi_n is None else np.asarray(psi_n, dtype=float).ravel()
    q_abs = None if q is None else np.abs(np.asarray(q, dtype=float).ravel())
    p = None if pressure is None else np.asarray(pressure, dtype=float).ravel()
    inside = np.isfinite(shear) & (np.abs(shear) < threshold)
    i_min = None if q_abs is None or not np.isfinite(q_abs).any() else int(np.nanargmin(q_abs))
    regions = []
    for i, j in _interval_edges(rho, inside):
        sl = slice(i, j + 1)
        q_lo = q_hi = math.nan
        contains_q1 = False
        if q_abs is not None:
            q_lo, q_hi = float(np.nanmin(q_abs[sl])), float(np.nanmax(q_abs[sl]))
            contains_q1 = bool(q_lo <= 1.0 <= q_hi)
        pressure_fields: dict[str, Optional[float]] = {}
        if p is not None and np.isfinite(p[sl]).all():
            pressure_fields["pressure_drop"] = float(p[i] - p[j])
            if j > i:
                gradient = np.gradient(p[sl], rho[sl])
                pressure_fields["mean_abs_pressure_gradient"] = float(np.mean(np.abs(gradient)))
                pressure_fields["max_abs_pressure_gradient"] = float(np.max(np.abs(gradient)))
        regions.append(LowShearRegion(
            threshold=float(threshold),
            rho_start=float(rho[i]),
            rho_end=float(rho[j]),
            psi_n_start=float(psi[i]),
            psi_n_end=float(psi[j]),
            contains_axis=bool(i == 0),
            contains_q_min=bool(i_min is not None and i <= i_min <= j),
            contains_q1_surface=contains_q1,
            q_low=q_lo,
            q_high=q_hi,
            mean_shear=float(np.mean(shear[sl])),
            max_abs_shear=float(np.max(np.abs(shear[sl]))),
            **pressure_fields,
        ))
    return tuple(regions)