Repository: astropy/astropy

Issue:
Consider changing the default behavior of `Table.pformat()` to be independent of the terminal
### What is the problem this feature will solve?

Currently the `Table.pformat()` default is to respect the size of the interactive terminal when generating lines of output representing the formatted table. This matches what would be seen with the corresponding `Table.pprint()` call (more-or-less).

The problem is that in most cases `pformat()` is used in a context where an interactive terminal is not relevant or does not exist (in which case 24x80 is assumed). This leads to unexpected results or the need for workarounds. A lot of test code uses `pformat_all()` as a work-around. 

Issue #14010 is a direct result of this unexpected behavior because `str(table)` also respects the terminal size.

### Describe the desired outcome

By default `Table.pformat()` would be equivalent to `Table.pformat_all()`, i.e. use `max_lines=-1` and `max_width=-1`. This would make `pformat()` independent of the terminal and lead to fewer suprises.

### Additional context

The big concern is whether this change will cause unexpected code breakage. This would need to wait until 7.0 but if merged into `main` there would be plenty of time for testing.

Discussion:
github-actions[bot]: Thank you for your contribution to Astropy! 🌌 This checklist is meant to remind the package maintainers who will review this pull request of some common things to look for.

  - [ ] Do the proposed changes actually accomplish desired goals?
  - [ ] Do the proposed changes follow the [Astropy coding guidelines](https://docs.astropy.org/en/latest/development/codeguide.html)?
  - [ ] Are tests added/updated as required? If so, do they follow the [Astropy testing guidelines](https://docs.astropy.org/en/latest/development/testguide.html)?
  - [ ] Are docs added/updated as required? If so, do they follow the [Astropy documentation guidelines](https://docs.astropy.org/en/latest/development/docguide.html)?
  - [ ] Is rebase and/or squash necessary? If so, please provide the author with appropriate instructions. Also see instructions for [rebase](https://docs.astropy.org/en/latest/development/development_details.html#rebase-if-necessary) and [squash](https://docs.astropy.org/en/latest/development/development_details.html#squash-if-necessary).
  - [ ] Did the CI pass? If no, are the failures related? If you need to run daily and weekly cron jobs as part of the PR, please apply the "Extra CI" label. Codestyle issues can be fixed by the [bot](https://docs.astropy.org/en/latest/development/development_details.html#pre-commit).
  - [ ] Is a change log needed? If yes, did the change log check pass? If no, add the "no-changelog-entry-needed" label. If this is a manual backport, use the "skip-changelog-checks" label unless special changelog handling is necessary.
  - [ ] Is this a big PR that makes a "What's new?" entry worthwhile and if so, is (1) a "what's new" entry included in this PR and (2) the "whatsnew-needed" label applied?
  - [ ] At the time of adding the milestone, if the milestone set requires a backport to release branch(es), apply the appropriate "backport-X.Y.x" label(s) *before* merge.
github-actions[bot]: 👋 Thank you for your draft pull request! Do you know that you can use `[ci skip]` or `[skip ci]` in your commit messages to skip running continuous integration tests until you are ready?
taldcroft: Needs a rebase. :-(
neutrinoceros: suggestions integrated and branch rebased !
taldcroft: Sorry, this needs a rebase as well. But now I have auto-merge enabled so it should go through!

Patch:
diff --git a/astropy/stats/nanfunctions.py b/astropy/stats/nanfunctions.py
index 098aba060c..e18bb10fac 100644
--- a/astropy/stats/nanfunctions.py
+++ b/astropy/stats/nanfunctions.py
@@ -77,7 +77,11 @@ if HAS_BOTTLENECK:
 
         result = function(array, axis=axis, **kwargs)
         if isinstance(array, Quantity):
-            return array.__array_wrap__(result)
+            if function is bottleneck.nanvar:
+                result = array._result_as_quantity(result, array.unit**2, None)
+            else:
+                result = array.__array_wrap__(result)
+            return result
         elif isinstance(result, float):
             # For compatibility with numpy, always return a numpy scalar.
             return np.float64(result)
@@ -86,16 +90,22 @@ if HAS_BOTTLENECK:
 
     bn_funcs = dict(
         nansum=functools.partial(_apply_bottleneck, bottleneck.nansum),
+        nanmin=functools.partial(_apply_bottleneck, bottleneck.nanmin),
+        nanmax=functools.partial(_apply_bottleneck, bottleneck.nanmax),
         nanmean=functools.partial(_apply_bottleneck, bottleneck.nanmean),
         nanmedian=functools.partial(_apply_bottleneck, bottleneck.nanmedian),
         nanstd=functools.partial(_apply_bottleneck, bottleneck.nanstd),
+        nanvar=functools.partial(_apply_bottleneck, bottleneck.nanvar),
     )
 
     np_funcs = dict(
         nansum=np.nansum,
+        nanmin=np.nanmin,
+        nanmax=np.nanmax,
         nanmean=np.nanmean,
         nanmedian=np.nanmedian,
         nanstd=np.nanstd,
+        nanvar=np.nanvar,
     )
 
     def _dtype_dispatch(func_name):
@@ -115,15 +125,21 @@ if HAS_BOTTLENECK:
         return wrapped
 
     nansum = _dtype_dispatch("nansum")
+    nanmin = _dtype_dispatch("nanmin")
+    nanmax = _dtype_dispatch("nanmax")
     nanmean = _dtype_dispatch("nanmean")
     nanmedian = _dtype_dispatch("nanmedian")
     nanstd = _dtype_dispatch("nanstd")
+    nanvar = _dtype_dispatch("nanvar")
 
 else:
     nansum = np.nansum
+    nanmin = np.nanmin
+    nanmax = np.nanmax
     nanmean = np.nanmean
     nanmedian = np.nanmedian
     nanstd = np.nanstd
+    nanvar = np.nanvar
 
 
 def nanmadstd(
diff --git a/astropy/stats/sigma_clipping.py b/astropy/stats/sigma_clipping.py
index 7e5885e86f..334c16a8f7 100644
--- a/astropy/stats/sigma_clipping.py
+++ b/astropy/stats/sigma_clipping.py
@@ -8,7 +8,18 @@ from typing import TYPE_CHECKING
 import numpy as np
 
 from astropy.stats._fast_sigma_clip import _sigma_clip_fast
-from astropy.stats.nanfunctions import nanmadstd, nanmean, nanmedian, nanstd
+from astropy.stats.biweight import biweight_location, biweight_scale
+from astropy.stats.funcs import mad_std
+from astropy.stats.nanfunctions import (
+    nanmadstd,
+    nanmax,
+    nanmean,
+    nanmedian,
+    nanmin,
+    nanstd,
+    nansum,
+    nanvar,
+)
 from astropy.units import Quantity
 from astropy.utils import isiterable
 from astropy.utils.compat.numpycompat import NUMPY_LT_2_0
@@ -25,7 +36,7 @@ if TYPE_CHECKING:
 
     from numpy.typing import ArrayLike, NDArray
 
-__all__ = ["SigmaClip", "sigma_clip", "sigma_clipped_stats"]
+__all__ = ["SigmaClip", "sigma_clip", "SigmaClippedStats", "sigma_clipped_stats"]
 
 
 class SigmaClip:
@@ -112,7 +123,7 @@ class SigmaClip:
 
     See Also
     --------
-    sigma_clip, sigma_clipped_stats
+    sigma_clip, sigma_clipped_stats, SigmaClippedStats
 
     Notes
     -----
@@ -817,7 +828,7 @@ def sigma_clip(
 
     See Also
     --------
-    SigmaClip, sigma_clipped_stats
+    SigmaClip, sigma_clipped_stats, SigmaClippedStats
 
     Notes
     -----
@@ -879,6 +890,330 @@ def sigma_clip(
     )
 
 
+class SigmaClippedStats:
+    """
+    Class to calculate sigma-clipped statistics on the provided data.
+
+    Parameters
+    ----------
+    data : array-like or `~numpy.ma.MaskedArray`
+        Data array or object that can be converted to an array.
+
+    mask : `numpy.ndarray` (bool), optional
+        A boolean mask with the same shape as ``data``, where a `True`
+        value indicates the corresponding element of ``data`` is masked.
+        Masked pixels are excluded when computing the statistics.
+
+    mask_value : float, optional
+        A data value (e.g., ``0.0``) that is ignored when computing the
+        statistics. ``mask_value`` will be masked in addition to any
+        input ``mask``.
+
+    sigma : float, optional
+        The number of standard deviations to use for both the lower
+        and upper clipping limit. These limits are overridden by
+        ``sigma_lower`` and ``sigma_upper``, if input. The default is 3.
+
+    sigma_lower : float or None, optional
+        The number of standard deviations to use as the lower bound for
+        the clipping limit. If `None` then the value of ``sigma`` is
+        used. The default is `None`.
+
+    sigma_upper : float or None, optional
+        The number of standard deviations to use as the upper bound for
+        the clipping limit. If `None` then the value of ``sigma`` is
+        used. The default is `None`.
+
+    maxiters : int or None, optional
+        The maximum number of sigma-clipping iterations to perform or
+        `None` to clip until convergence is achieved (i.e., iterate
+        until the last iteration clips nothing). If convergence is
+        achieved prior to ``maxiters`` iterations, the clipping
+        iterations will stop. The default is 5.
+
+    cenfunc : {'median', 'mean'} or callable, optional
+        The statistic or callable function/object used to compute
+        the center value for the clipping. If using a callable
+        function/object and the ``axis`` keyword is used, then it must
+        be able to ignore NaNs (e.g., `numpy.nanmean`) and it must have
+        an ``axis`` keyword to return an array with axis dimension(s)
+        removed. The default is ``'median'``.
+
+    stdfunc : {'std', 'mad_std'} or callable, optional
+        The statistic or callable function/object used to compute the
+        standard deviation about the center value. If using a callable
+        function/object and the ``axis`` keyword is used, then it must
+        be able to ignore NaNs (e.g., `numpy.nanstd`) and it must have
+        an ``axis`` keyword to return an array with axis dimension(s)
+        removed. The default is ``'std'``.
+
+    axis : None or int or tuple of int, optional
+        The axis or axes along which to sigma clip the data. If `None`,
+        then the flattened data will be used. ``axis`` is passed to the
+        ``cenfunc`` and ``stdfunc``. The default is `None`.
+
+    grow : float or `False`, optional
+        Radius within which to mask the neighbouring pixels of those
+        that fall outwith the clipping limits (only applied along
+        ``axis``, if specified). As an example, for a 2D image a value
+        of 1 will mask the nearest pixels in a cross pattern around each
+        deviant pixel, while 1.5 will also reject the nearest diagonal
+        neighbours and so on.
+
+    Notes
+    -----
+    The best performance will typically be obtained by setting
+    ``cenfunc`` and ``stdfunc`` to one of the built-in functions
+    specified as a string. If one of the options is set to a string
+    while the other has a custom callable, you may in some cases
+    see better performance if you have the `bottleneck`_ package
+    installed. To preserve accuracy, bottleneck is only used for float64
+    computations.
+
+    .. _bottleneck:  https://github.com/pydata/bottleneck
+
+    See Also
+    --------
+    sigma_clipped_stats, SigmaClip, sigma_clip
+    """
+
+    def __init__(
+        self,
+        data: ArrayLike,
+        *,
+        mask: NDArray | None = None,
+        mask_value: float | None = None,
+        sigma: float = 3.0,
+        sigma_lower: float | None = None,
+        sigma_upper: float | None = None,
+        maxiters: int = 5,
+        cenfunc: Literal["median", "mean"] | Callable = "median",
+        stdfunc: Literal["std", "mad_std"] | Callable = "std",
+        axis: int | tuple[int, ...] | None = None,
+        grow: float | Literal[False] | None = False,
+    ) -> None:
+        sigclip = SigmaClip(
+            sigma=sigma,
+            sigma_lower=sigma_lower,
+            sigma_upper=sigma_upper,
+            maxiters=maxiters,
+            cenfunc=cenfunc,
+            stdfunc=stdfunc,
+            grow=grow,
+        )
+
+        if mask is not None:
+            data = np.ma.MaskedArray(data, mask)
+        if mask_value is not None:
+            data = np.ma.masked_values(data, mask_value)
+
+        if isinstance(data, np.ma.MaskedArray) and data.mask.all():
+            raise ValueError("input data is all masked")
+
+        self.data = sigclip(
+            data, axis=axis, masked=False, return_bounds=False, copy=True
+        )
+        self.axis = axis
+
+    def min(self) -> float | NDArray:
+        """
+        Calculate the minimum of the data.
+
+        NaN values are ignored.
+
+        Returns
+        -------
+        min : float or `~numpy.ndarray`
+            The minimum of the data.
+        """
+        return nanmin(self.data, axis=self.axis)
+
+    def max(self) -> float | NDArray:
+        """
+        Calculate the maximum of the data.
+
+        NaN values are ignored.
+
+        Returns
+        -------
+        max : float or `~numpy.ndarray`
+            The maximum of the data.
+        """
+        return nanmax(self.data, axis=self.axis)
+
+    def sum(self) -> float | NDArray:
+        """
+        Calculate the sum of the data.
+
+        NaN values are ignored.
+
+        Returns
+        -------
+        sum : float or `~numpy.ndarray`
+            The sum of the data.
+        """
+        return nansum(self.data, axis=self.axis)
+
+    def mean(self) -> float | NDArray:
+        """
+        Calculate the mean of the data.
+
+        NaN values are ignored.
+
+        Returns
+        -------
+        mean : float or `~numpy.ndarray`
+            The mean of the data.
+        """
+        return nanmean(self.data, axis=self.axis)
+
+    def median(self) -> float | NDArray:
+        """
+        Calculate the median of the data.
+
+        NaN values are ignored.
+
+        Returns
+        -------
+        median : float or `~numpy.ndarray`
+            The median of the data.
+        """
+        return nanmedian(self.data, axis=self.axis)
+
+    def mode(
+        self, median_factor: float = 3.0, mean_factor: float = 2.0
+    ) -> float | NDArray:
+        """
+        Calculate the mode of the data using a estimator of the form
+        ``(median_factor * median) - (mean_factor * mean)``.
+
+        NaN values are ignored.
+
+        Parameters
+        ----------
+        median_factor : float, optional
+            The multiplicative factor for the data median. Defaults to 3.
+
+        mean_factor : float, optional
+            The multiplicative factor for the data mean. Defaults to 2.
+
+        Returns
+        -------
+        mode : float or `~numpy.ndarray`
+            The estimated mode of the data.
+        """
+        return (median_factor * self.median()) - (mean_factor * self.mean())
+
+    def std(self, ddof: int = 0) -> float | NDArray:
+        """
+        Calculate the standard deviation of the data.
+
+        NaN values are ignored.
+
+        Parameters
+        ----------
+        ddof : int, optional
+            The delta degrees of freedom for the standard deviation
+            calculation. The divisor used in the calculation is ``N -
+            ddof``, where ``N`` represents the number of elements. For
+            a population standard deviation where you have data for the
+            entire population, use ``ddof=0``. For a sample standard
+            deviation where you have a sample of the population, use
+            ``ddof=1``. The default is 0.
+
+        Returns
+        -------
+        std : float or `~numpy.ndarray`
+            The standard deviation of the data.
+        """
+        return nanstd(self.data, axis=self.axis, ddof=ddof)
+
+    def var(self, ddof: int = 0) -> float | NDArray:
+        """
+        Calculate the variance of the data.
+
+        NaN values are ignored.
+
+        Parameters
+        ----------
+        ddof : int, optional
+            The delta degrees of freedom. The divisor used in the
+            calculation is ``N - ddof``, where ``N`` represents the
+            number of elements. For a population variance where you have
+            data for the entire population, use ``ddof=0``. For a sample
+            variance where you have a sample of the population, use
+            ``ddof=1``. The default is 0.
+
+        Returns
+        -------
+        var : float or `~numpy.ndarray`
+            The variance of the data.
+        """
+        return nanvar(self.data, axis=self.axis, ddof=ddof)
+
+    def biweight_location(
+        self, c: float = 6.0, M: float | None = None
+    ) -> float | NDArray:
+        """
+        Calculate the biweight location of the data.
+
+        NaN values are ignored.
+
+        Parameters
+        ----------
+        c : float, optional
+            Tuning constant for the biweight estimator. Default value is
+            6.0.
+
+        M : float or None, optional
+            Initial guess for the biweight location. Default value is
+            `None`.
+
+        Returns
+        -------
+        biweight_location : float or `~numpy.ndarray`
+            The biweight location of the data.
+        """
+        return biweight_location(self.data, c=c, M=M, axis=self.axis, ignore_nan=True)
+
+    def biweight_scale(self, c: float = 6.0, M: float | None = None) -> float | NDArray:
+        """
+        Calculate the biweight scale of the data.
+
+        NaN values are ignored.
+
+        Parameters
+        ----------
+        c : float, optional
+            Tuning constant for the biweight estimator. Default value is
+            6.0.
+
+        M : float or None, optional
+            Initial guess for the biweight location. Default value is
+            `None`.
+
+        Returns
+        -------
+        biweight_scale : float or `~numpy.ndarray`
+            The biweight scale of the data.
+        """
+        return biweight_scale(self.data, c=c, M=M, axis=self.axis, ignore_nan=True)
+
+    def mad_std(self) -> float | NDArray:
+        """
+        Calculate the median absolute deviation (MAD) based standard
+        deviation of the data.
+
+        NaN values are ignored.
+
+        Returns
+        -------
+        mad_std : float or `~numpy.ndarray`
+            The MAD-based standard deviation of the data.
+        """
+        return mad_std(self.data, axis=self.axis, ignore_nan=True)
+
+
 def sigma_clipped_stats(
     data: ArrayLike,
     mask: NDArray | None = None,
@@ -952,8 +1287,11 @@ def sigma_clipped_stats(
     std_ddof : int, optional
         The delta degrees of freedom for the standard deviation
         calculation. The divisor used in the calculation is ``N -
-        std_ddof``, where ``N`` represents the number of elements. The
-        default is 0.
+        std_ddof``, where ``N`` represents the number of elements. For a
+        population standard deviation where you have data for the entire
+        population, use ``std_ddof=0``. For a sample standard deviation
+        where you have a sample of the population, use ``std_ddof=1``.
+        The default is 0.
 
     axis : None or int or tuple of int, optional
         The axis or axes along which to sigma clip the data. If `None`,
@@ -987,7 +1325,7 @@ def sigma_clipped_stats(
 
     See Also
     --------
-    SigmaClip, sigma_clip
+    SigmaClippedStats, SigmaClip, sigma_clip
     """
     if mask is not None:
         data = np.ma.MaskedArray(data, mask)
@@ -997,21 +1335,16 @@ def sigma_clipped_stats(
     if isinstance(data, np.ma.MaskedArray) and data.mask.all():
         return np.ma.masked, np.ma.masked, np.ma.masked
 
-    sigclip = SigmaClip(
+    stats = SigmaClippedStats(
+        data,
         sigma=sigma,
         sigma_lower=sigma_lower,
         sigma_upper=sigma_upper,
         maxiters=maxiters,
         cenfunc=cenfunc,
         stdfunc=stdfunc,
+        axis=axis,
         grow=grow,
     )
-    data_clipped = sigclip(
-        data, axis=axis, masked=False, return_bounds=False, copy=True
-    )
-
-    mean = nanmean(data_clipped, axis=axis)
-    median = nanmedian(data_clipped, axis=axis)
-    std = nanstd(data_clipped, ddof=std_ddof, axis=axis)
 
-    return mean, median, std
+    return stats.mean(), stats.median(), stats.std(ddof=std_ddof)
diff --git a/astropy/table/column.py b/astropy/table/column.py
index e8812ee0a5..1cf4a9b1a4 100644
--- a/astropy/table/column.py
+++ b/astropy/table/column.py
@@ -852,7 +852,7 @@ class BaseColumn(_ColumnGetitemShim, np.ndarray):
 
     def pformat(
         self,
-        max_lines=None,
+        max_lines=-1,
         show_name=True,
         show_unit=False,
         show_dtype=False,
@@ -860,17 +860,19 @@ class BaseColumn(_ColumnGetitemShim, np.ndarray):
     ):
         """Return a list of formatted string representation of column values.
 
-        If no value of ``max_lines`` is supplied then the height of the
+        If ``max_lines=None`` is supplied then the height of the
         screen terminal is used to set ``max_lines``.  If the terminal
         height cannot be determined then the default will be
         determined using the ``astropy.conf.max_lines`` configuration
         item. If a negative value of ``max_lines`` is supplied then
-        there is no line limit applied.
+        there is no line limit applied (default).
 
         Parameters
         ----------
-        max_lines : int
-            Maximum lines of output (header + data rows)
+        max_lines : int or None
+            Maximum lines of output (header + data rows).
+            -1 (default) implies no limit, ``None`` implies using the
+            height of the current terminal.
 
         show_name : bool
             Include column name. Default is True.
@@ -904,7 +906,7 @@ class BaseColumn(_ColumnGetitemShim, np.ndarray):
     def pprint(self, max_lines=None, show_name=True, show_unit=False, show_dtype=False):
         """Print a formatted string representation of column values.
 
-        If no value of ``max_lines`` is supplied then the height of the
+        If ``max_lines=None`` (default) then the height of the
         screen terminal is used to set ``max_lines``.  If the terminal
         height cannot be determined then the default will be
         determined using the ``astropy.conf.max_lines`` configuration
diff --git a/astropy/table/pprint.py b/astropy/table/pprint.py
index a4512fcca7..c782b79eac 100644
--- a/astropy/table/pprint.py
+++ b/astropy/table/pprint.py
@@ -557,8 +557,8 @@ class TableFormatter:
     def _pformat_table(
         self,
         table,
-        max_lines=None,
-        max_width=None,
+        max_lines=-1,
+        max_width=-1,
         show_name=True,
         show_unit=None,
         show_dtype=False,
@@ -574,9 +574,13 @@ class TableFormatter:
         ----------
         max_lines : int or None
             Maximum number of rows to output
+            -1 (default) implies no limit, ``None`` implies using the
+            height of the current terminal.
 
         max_width : int or None
             Maximum character width of output
+            -1 (default) implies no limit, ``None`` implies using the
+            width of the current terminal.
 
         show_name : bool
             Include a header row for column names. Default is True.
diff --git a/astropy/table/table.py b/astropy/table/table.py
index a5338f1024..0a29c28400 100644
--- a/astropy/table/table.py
+++ b/astropy/table/table.py
@@ -1649,7 +1649,7 @@ class Table:
         return self._base_repr_(html=False, max_width=None)
 
     def __str__(self):
-        return "\n".join(self.pformat())
+        return "\n".join(self.pformat(max_lines=None, max_width=None))
 
     def __bytes__(self):
         return str(self).encode("utf-8")
@@ -1919,8 +1919,8 @@ class Table:
     @format_doc(_pformat_docs, id="{id}")
     def pformat(
         self,
-        max_lines=None,
-        max_width=None,
+        max_lines=-1,
+        max_width=-1,
         show_name=True,
         show_unit=None,
         show_dtype=False,
@@ -1932,12 +1932,12 @@ class Table:
         """Return a list of lines for the formatted string representation of
         the table.
 
-        If no value of ``max_lines`` is supplied then the height of the
+        If ``max_lines=None`` is supplied then the height of the
         screen terminal is used to set ``max_lines``.  If the terminal
-        height cannot be determined then the default is taken from the
-        configuration item ``astropy.conf.max_lines``.  If a negative
-        value of ``max_lines`` is supplied then there is no line limit
-        applied.
+        height cannot be determined then the default will be
+        determined using the ``astropy.conf.max_lines`` configuration
+        item. If a negative value of ``max_lines`` is supplied then
+        there is no line limit applied (default).
 
         The same applies for ``max_width`` except the configuration item  is
         ``astropy.conf.max_width``.
@@ -1961,6 +1961,7 @@ class Table:
 
         return lines
 
+    @deprecated(since="7.0", alternative="Table.pformat")
     @format_doc(_pformat_docs, id="{id}")
     def pformat_all(
         self,
@@ -1977,12 +1978,12 @@ class Table:
         """Return a list of lines for the formatted string representation of
         the entire table.
 
-        If no value of ``max_lines`` is supplied then the height of the
+        If ``max_lines=None`` is supplied then the height of the
         screen terminal is used to set ``max_lines``.  If the terminal
-        height cannot be determined then the default is taken from the
-        configuration item ``astropy.conf.max_lines``.  If a negative
-        value of ``max_lines`` is supplied then there is no line limit
-        applied.
+        height cannot be determined then the default will be
+        determined using the ``astropy.conf.max_lines`` configuration
+        item. If a negative value of ``max_lines`` is supplied then
+        there is no line limit applied (default).
 
         The same applies for ``max_width`` except the configuration item  is
         ``astropy.conf.max_width``.
diff --git a/astropy/units/astrophys.py b/astropy/units/astrophys.py
index 3e16d5ba34..b942884f9f 100644
--- a/astropy/units/astrophys.py
+++ b/astropy/units/astrophys.py
@@ -114,18 +114,6 @@ def_unit(
 ##########################################################################
 # ENERGY
 
-# Here, explicitly convert the planck constant to 'eV s' since the constant
-# can override that to give a more precise value that takes into account
-# covariances between e and h.  Eventually, this may also be replaced with
-# just `_si.Ryd.to(eV)`.
-def_unit(
-    ["Ry", "rydberg"],
-    (_si.Ryd * _si.c * _si.h.to(si.eV * si.s)).to(si.eV),
-    namespace=_ns,
-    prefixes=True,
-    doc="Rydberg: Energy of a photon whose wavenumber is the Rydberg constant",
-    format={"latex": r"R_{\infty}", "unicode": "R∞"},
-)
 def_unit(
     ["foe", "Bethe", "bethe"],
     1e51 * si.g * si.cm**2 / si.s**2,
diff --git a/astropy/units/equivalencies.py b/astropy/units/equivalencies.py
index 4e226ecc2e..f89a1942a5 100644
--- a/astropy/units/equivalencies.py
+++ b/astropy/units/equivalencies.py
@@ -386,11 +386,11 @@ def doppler_radio(rest):
         return restwav * ckms / (ckms - x)
 
     def to_vel_en(x):
-        resten = rest.to_value(si.eV, equivalencies=spectral())
+        resten = rest.to_value(misc.eV, equivalencies=spectral())
         return (resten - x) / (resten) * ckms
 
     def from_vel_en(x):
-        resten = rest.to_value(si.eV, equivalencies=spectral())
+        resten = rest.to_value(misc.eV, equivalencies=spectral())
         voverc = x / ckms
         return resten * (1 - voverc)
 
@@ -398,7 +398,7 @@ def doppler_radio(rest):
         [
             (si.Hz, si.km / si.s, to_vel_freq, from_vel_freq),
             (si.AA, si.km / si.s, to_vel_wav, from_vel_wav),
-            (si.eV, si.km / si.s, to_vel_en, from_vel_en),
+            (misc.eV, si.km / si.s, to_vel_en, from_vel_en),
         ],
         "doppler_radio",
         {"rest": rest},
@@ -456,11 +456,11 @@ def doppler_optical(rest):
         return restwav * (1 + voverc)
 
     def to_vel_en(x):
-        resten = rest.to_value(si.eV, equivalencies=spectral())
+        resten = rest.to_value(misc.eV, equivalencies=spectral())
         return ckms * (resten - x) / x
 
     def from_vel_en(x):
-        resten = rest.to_value(si.eV, equivalencies=spectral())
+        resten = rest.to_value(misc.eV, equivalencies=spectral())
         voverc = x / ckms
         return resten / (1 + voverc)
 
@@ -468,7 +468,7 @@ def doppler_optical(rest):
         [
             (si.Hz, si.km / si.s, to_vel_freq, from_vel_freq),
             (si.AA, si.km / si.s, to_vel_wav, from_vel_wav),
-            (si.eV, si.km / si.s, to_vel_en, from_vel_en),
+            (misc.eV, si.km / si.s, to_vel_en, from_vel_en),
         ],
         "doppler_optical",
         {"rest": rest},
@@ -533,11 +533,11 @@ def doppler_relativistic(rest):
         return restwav * ((1 + voverc) / (1 - voverc)) ** 0.5
 
     def to_vel_en(x):
-        resten = rest.to_value(si.eV, spectral())
+        resten = rest.to_value(misc.eV, spectral())
         return (resten**2 - x**2) / (resten**2 + x**2) * ckms
 
     def from_vel_en(x):
-        resten = rest.to_value(si.eV, spectral())
+        resten = rest.to_value(misc.eV, spectral())
         voverc = x / ckms
         return resten * ((1 - voverc) / (1 + (voverc))) ** 0.5
 
@@ -545,7 +545,7 @@ def doppler_relativistic(rest):
         [
             (si.Hz, si.km / si.s, to_vel_freq, from_vel_freq),
             (si.AA, si.km / si.s, to_vel_wav, from_vel_wav),
-            (si.eV, si.km / si.s, to_vel_en, from_vel_en),
+            (misc.eV, si.km / si.s, to_vel_en, from_vel_en),
         ],
         "doppler_relativistic",
         {"rest": rest},
@@ -817,7 +817,7 @@ def temperature_energy():
     e = _si.e.value
     k_B = _si.k_B.value
     return Equivalency(
-        [(si.K, si.eV, lambda x: x / (e / k_B), lambda x: x * (e / k_B))],
+        [(si.K, misc.eV, lambda x: x / (e / k_B), lambda x: x * (e / k_B))],
         "temperature_energy",
     )
 
diff --git a/astropy/units/misc.py b/astropy/units/misc.py
index ff3440bdfc..2fa55c6c8b 100644
--- a/astropy/units/misc.py
+++ b/astropy/units/misc.py
@@ -104,6 +104,30 @@ def_unit(
     doc="Unified atomic mass unit",
 )
 
+##########################################################################
+# ENERGY
+
+def_unit(
+    ["eV", "electronvolt"],
+    _si.e.value * si.J,
+    namespace=_ns,
+    prefixes=True,
+    doc="Electron Volt",
+)
+
+# Here, explicitly convert the planck constant to 'eV s' since the constant
+# can override that to give a more precise value that takes into account
+# covariances between e and h.  Eventually, this may also be replaced with
+# just `_si.Ryd.to(eV)`.
+def_unit(
+    ["Ry", "rydberg"],
+    (_si.Ryd * _si.c * _si.h.to(eV * si.s)).to(eV),
+    namespace=_ns,
+    prefixes=True,
+    doc="Rydberg: Energy of a photon whose wavenumber is the Rydberg constant",
+    format={"latex": r"R_{\infty}", "unicode": "R∞"},
+)
+
 
 ###########################################################################
 # COMPUTER
diff --git a/astropy/units/quantity_helper/function_helpers.py b/astropy/units/quantity_helper/function_helpers.py
index c515d28fc9..be464c573b 100644
--- a/astropy/units/quantity_helper/function_helpers.py
+++ b/astropy/units/quantity_helper/function_helpers.py
@@ -97,8 +97,7 @@ SUBCLASS_SAFE_FUNCTIONS |= {
     np.isposinf, np.isneginf, np.isreal, np.iscomplex,
     np.average, np.mean, np.std, np.var, np.trace,
     np.nanmax, np.nanmin, np.nanargmin, np.nanargmax, np.nanmean,
-    np.nansum, np.nancumsum, np.nanstd, np.nanvar,
-    np.nanprod, np.nancumprod,
+    np.nansum, np.nancumsum, np.nanprod, np.nancumprod,
     np.einsum_path, np.linspace,
     np.sort, np.partition, np.meshgrid,
     np.common_type, np.result_type, np.can_cast, np.min_scalar_type,
@@ -219,6 +218,7 @@ dispatched_function = FunctionAssigner(DISPATCHED_FUNCTIONS)
         np.fft.fft2, np.fft.ifft2, np.fft.rfft2, np.fft.irfft2,
         np.fft.fftn, np.fft.ifftn, np.fft.rfftn, np.fft.irfftn,
         np.fft.hfft, np.fft.ihfft,
+        np.nanstd,  # See comment on nanvar helper.
         np.linalg.eigvals, np.linalg.eigvalsh,
     } | ({np.asfarray} if NUMPY_LT_2_0 else set())  # noqa: NPY201
 )  # fmt: skip
@@ -252,6 +252,17 @@ def like_helper(a, *args, **kwargs):
     return (a.view(np.ndarray),) + args, kwargs, unit, None
 
 
+# nanvar is safe for Quantity and was previously in SUBCLASS_FUNCTIONS, but it
+# is not safe for Angle, since the resulting unit is inconsistent with being
+# an Angle. By using FUNCTION_HELPERS, the unit gets passed through
+# _result_as_quantity, which will correctly drop to Quantity.
+# A side effect would be that np.nanstd then also produces Quantity; this
+# is avoided by it being helped by invariant_a_helpers above.
+@function_helper
+def nanvar(a, *args, **kwargs):
+    return (a.view(np.ndarray),) + args, kwargs, a.unit**2, None
+
+
 @function_helper
 def sinc(x):
     from astropy.units.si import radian
diff --git a/astropy/units/si.py b/astropy/units/si.py
index 08aaefdd41..7a1201044e 100644
--- a/astropy/units/si.py
+++ b/astropy/units/si.py
@@ -8,8 +8,6 @@ This package defines the SI units.  They are also available in
 
 import numpy as np
 
-from astropy.constants import si as _si
-
 from .core import Unit, UnitBase, def_unit
 
 __all__: list[str] = []  #  Units are added at the end
@@ -61,7 +59,12 @@ def_unit(
     namespace=_ns,
     doc="ångström: 10 ** -10 m",
     prefixes=[(["m", "milli"], ["milli", "m"], 1.0e-3)],
-    format={"latex": r"\mathring{A}", "unicode": "Å", "vounit": "Angstrom"},
+    format={
+        "latex": r"\mathring{A}",
+        "ogip": "angstrom",
+        "unicode": "Å",
+        "vounit": "Angstrom",
+    },
 )
 
 
@@ -310,13 +313,6 @@ def_unit(
     prefixes=True,
     doc="Joule: energy",
 )
-def_unit(
-    ["eV", "electronvolt"],
-    _si.e.value * J,
-    namespace=_ns,
-    prefixes=True,
-    doc="Electron Volt",
-)
 
 
 ##########################################################################
diff --git a/astropy/visualization/interval.py b/astropy/visualization/interval.py
index 66fac6a117..4f4d6253de 100644
--- a/astropy/visualization/interval.py
+++ b/astropy/visualization/interval.py
@@ -146,19 +146,25 @@ class AsymmetricPercentileInterval(BaseInterval):
 
     Parameters
     ----------
-    lower_percentile : float
-        The lower percentile below which to ignore pixels.
-    upper_percentile : float
-        The upper percentile above which to ignore pixels.
+    lower_percentile : float or None
+        The lower percentile below which to ignore pixels. If None, then
+        defaults to 0.
+    upper_percentile : float or None
+        The upper percentile above which to ignore pixels. If None, then
+        defaults to 100.
     n_samples : int, optional
         Maximum number of values to use. If this is specified, and there
         are more values in the dataset as this, then values are randomly
         sampled from the array (with replacement).
     """
 
-    def __init__(self, lower_percentile, upper_percentile, n_samples=None):
-        self.lower_percentile = lower_percentile
-        self.upper_percentile = upper_percentile
+    def __init__(self, lower_percentile=None, upper_percentile=None, n_samples=None):
+        self.lower_percentile = (
+            lower_percentile if lower_percentile is not None else 0.0
+        )
+        self.upper_percentile = (
+            upper_percentile if upper_percentile is not None else 100.0
+        )
         self.n_samples = n_samples
 
     def get_limits(self, values):
diff --git a/astropy/visualization/mpl_normalize.py b/astropy/visualization/mpl_normalize.py
index 78a10a35a2..356a9b186f 100644
--- a/astropy/visualization/mpl_normalize.py
+++ b/astropy/visualization/mpl_normalize.py
@@ -37,7 +37,7 @@ else:
             raise ImportError("matplotlib is required in order to use this class.")
 
 
-__all__ = ["ImageNormalize", "simple_norm", "imshow_norm"]
+__all__ = ["ImageNormalize", "SimpleNorm", "simple_norm", "imshow_norm"]
 
 __doctest_requires__ = {"*": ["matplotlib"]}
 
@@ -212,6 +212,229 @@ class ImageNormalize(Normalize):
         return values_norm * (self.vmax - self.vmin) + self.vmin
 
 
+class SimpleNorm:
+    """
+    Class to create a normalization object that can be used for
+    displaying images with Matplotlib.
+
+    This convenience class provides the most common image stretching
+    functions. Additional stretch functions are available in
+    `~astropy.visualization.mpl_normalize.ImageNormalize`.
+
+    Parameters
+    ----------
+    stretch : {'linear', 'sqrt', 'power', log', 'asinh', 'sinh'}, optional
+        The stretch function to apply to the image. The default is
+        'linear'.
+
+    percent : float, optional
+        The percentage of the image values used to determine the pixel
+        values of the minimum and maximum cut levels. The lower cut
+        level will set at the ``(100 - percent) / 2`` percentile, while
+        the upper cut level will be set at the ``(100 + percent) / 2``
+        percentile. The default is 100.0. ``percent`` is ignored if
+        either ``min_percent`` or ``max_percent`` is input.
+
+    min_percent : float, optional
+        The percentile value used to determine the pixel value of
+        minimum cut level. The default is 0.0. ``min_percent`` overrides
+        ``percent``.
+
+    max_percent : float, optional
+        The percentile value used to determine the pixel value of
+        maximum cut level. The default is 100.0. ``max_percent``
+        overrides ``percent``.
+
+    vmin : float, optional
+        The pixel value of the minimum cut level. Data values less
+        than ``vmin`` will set to ``vmin`` before stretching the
+        image. The default is the image minimum. ``vmin`` overrides
+        ``min_percent``.
+
+    vmax : float, optional
+        The pixel value of the maximum cut level. Data values greater
+        than ``vmax`` will set to ``vmax`` before stretching the
+        image. The default is the image maximum. ``vmax`` overrides
+        ``max_percent``.
+
+    power : float, optional
+        The power index for ``stretch='power'``. The default is 1.0.
+
+    log_a : float, optional
+        The log index for ``stretch='log'``. The default is 1000.
+
+    asinh_a : float, optional
+        For ``stretch='asinh'``, the value where the asinh curve
+        transitions from linear to logarithmic behavior, expressed as a
+        fraction of the normalized image. Must be in the range between 0
+        and 1. The default is 0.1.
+
+    sinh_a : float, optional
+        The scaling parameter for ``stretch='sinh'``. The default is
+        0.3.
+
+    clip : bool, optional
+        If `True`, data values outside the [0:1] range are clipped to
+        the [0:1] range.
+
+    invalid : None or float, optional
+        Value to assign NaN values generated by the normalization. NaNs
+        in the input ``data`` array are not changed. For matplotlib
+        normalization, the ``invalid`` value should map to the
+        matplotlib colormap "under" value (i.e., any finite value < 0).
+        If `None`, then NaN values are not replaced. This keyword has no
+        effect if ``clip=True``.
+
+    See Also
+    --------
+    simple_norm
+
+    Examples
+    --------
+    .. plot::
+        :include-source:
+
+        import numpy as np
+        import matplotlib.pyplot as plt
+        from astropy.visualization import SimpleNorm
+
+        image = np.arange(65536).reshape((256, 256))
+        snorm = SimpleNorm('sqrt', percent=98)
+        norm = snorm(image)
+        fig, ax = plt.subplots()
+        axim = ax.imshow(image, norm=norm, origin='lower')
+        fig.colorbar(axim)
+    """
+
+    def __init__(
+        self,
+        stretch="linear",
+        percent=None,
+        *,
+        min_percent=None,
+        max_percent=None,
+        vmin=None,
+        vmax=None,
+        power=1.0,
+        log_a=1000,
+        asinh_a=0.1,
+        sinh_a=0.3,
+        clip=False,
+        invalid=-1.0,
+    ):
+        if percent is not None:
+            interval = PercentileInterval(percent)
+        elif min_percent is not None or max_percent is not None:
+            interval = AsymmetricPercentileInterval(
+                lower_percentile=min_percent, upper_percentile=max_percent
+            )
+        elif vmin is not None or vmax is not None:
+            interval = ManualInterval(vmin, vmax)
+        else:
+            interval = MinMaxInterval()
+        self.interval = interval
+
+        if stretch == "linear":
+            stretch = LinearStretch()
+        elif stretch == "sqrt":
+            stretch = SqrtStretch()
+        elif stretch == "power":
+            stretch = PowerStretch(power)
+        elif stretch == "log":
+            stretch = LogStretch(log_a)
+        elif stretch == "asinh":
+            stretch = AsinhStretch(asinh_a)
+        elif stretch == "sinh":
+            stretch = SinhStretch(sinh_a)
+        else:
+            raise ValueError(f"Unknown stretch: {stretch}.")
+        self.stretch = stretch
+
+        self.clip = clip
+        self.invalid = invalid
+
+    def __call__(self, data):
+        """
+        Return an `ImageNormalize` instance that can be used for displaying
+        images with Matplotlib.
+
+        Parameters
+        ----------
+        data : ndarray
+            The image array.
+
+        Returns
+        -------
+        result : `ImageNormalize` instance
+            An `ImageNormalize` instance that can be used for
+            displaying images with Matplotlib.
+        """
+        vmin, vmax = self.interval.get_limits(data)
+        return ImageNormalize(
+            vmin=vmin,
+            vmax=vmax,
+            stretch=self.stretch,
+            clip=self.clip,
+            invalid=self.invalid,
+        )
+
+    def imshow(self, data, ax=None, **kwargs):
+        """
+        A convenience function to display an image using matplotlib's
+        `matplotlib.pyplot.imshow` function with the normalization
+        defined by this class.
+
+        Parameters
+        ----------
+        data : 2D or 3D array-like
+            The data to display. Can be whatever
+            `~matplotlib.pyplot.imshow` and `ImageNormalize` both
+            accept.
+
+        ax : None or `~matplotlib.axes.Axes`, optional
+            The matplotlib axes on which to plot. If `None`, then the
+            current `~matplotlib.axes.Axes` instance is used.
+
+        **kwargs : dict, optional
+            Keywords arguments passed to `~matplotlib.pyplot.imshow`.
+            Cannot include the ``norm`` or ``X`` keyword.
+
+        Returns
+        -------
+        result : `~matplotlib.image.AxesImage`
+            The `~matplotlib.image.AxesImage` generated by
+            `~matplotlib.pyplot.imshow`.
+
+        Examples
+        --------
+        .. plot::
+            :include-source:
+
+            import numpy as np
+            import matplotlib.pyplot as plt
+            from astropy.visualization import SimpleNorm
+
+            image = np.arange(65536).reshape((256, 256))
+            snorm = SimpleNorm('sqrt', percent=98)
+            fig, ax = plt.subplots()
+            axim = snorm.imshow(image, ax=ax, origin='lower')
+            fig.colorbar(axim)
+        """
+        import matplotlib.pyplot as plt
+
+        if ax is None:
+            ax = plt.gca()
+
+        if "norm" in kwargs:
+            raise ValueError(
+                "This class already defines the norm. Use "
+                "matplotlib.pyplot.imshow directly to use your norm."
+            )
+
+        axim = ax.imshow(data, norm=self(data), **kwargs)
+        return axim
+
+
 @deprecated_renamed_argument(["min_cut", "max_cut"], ["vmin", "vmax"], ["6.1", "6.1"])
 def simple_norm(
     data,
@@ -309,38 +532,26 @@ def simple_norm(
     result : `ImageNormalize` instance
         An `ImageNormalize` instance that can be used for displaying
         images with Matplotlib.
+
+    See Also
+    --------
+    SimpleNorm
     """
-    if percent is not None:
-        interval = PercentileInterval(percent)
-    elif min_percent is not None or max_percent is not None:
-        interval = AsymmetricPercentileInterval(
-            min_percent or 0.0, max_percent or 100.0
-        )
-    elif vmin is not None or vmax is not None:
-        interval = ManualInterval(vmin, vmax)
-    else:
-        interval = MinMaxInterval()
-
-    if stretch == "linear":
-        stretch = LinearStretch()
-    elif stretch == "sqrt":
-        stretch = SqrtStretch()
-    elif stretch == "power":
-        stretch = PowerStretch(power)
-    elif stretch == "log":
-        stretch = LogStretch(log_a)
-    elif stretch == "asinh":
-        stretch = AsinhStretch(asinh_a)
-    elif stretch == "sinh":
-        stretch = SinhStretch(sinh_a)
-    else:
-        raise ValueError(f"Unknown stretch: {stretch}.")
-
-    vmin, vmax = interval.get_limits(data)
-
-    return ImageNormalize(
-        vmin=vmin, vmax=vmax, stretch=stretch, clip=clip, invalid=invalid
+    simple_norm = SimpleNorm(
+        stretch=stretch,
+        percent=percent,
+        min_percent=min_percent,
+        max_percent=max_percent,
+        vmin=vmin,
+        vmax=vmax,
+        power=power,
+        log_a=log_a,
+        asinh_a=asinh_a,
+        sinh_a=sinh_a,
+        clip=clip,
+        invalid=invalid,
     )
+    return simple_norm(data)
 
 
 # used in imshow_norm
@@ -402,9 +613,6 @@ def imshow_norm(data, ax=None, **kwargs):
 
         ax = plt.gca()
 
-    if "X" in kwargs:
-        raise ValueError("Cannot give both ``X`` and ``data``")
-
     if "norm" in kwargs:
         raise ValueError(
             "There is no point in using imshow_norm if you give "
diff --git a/docs/changes/coordinates/17239.bugfix.rst b/docs/changes/coordinates/17239.bugfix.rst
new file mode 100644
index 0000000000..34f071f1f9
--- /dev/null
+++ b/docs/changes/coordinates/17239.bugfix.rst
@@ -0,0 +1,2 @@
+``np.nanvar(angle)`` now produces a ``Quantity`` with the correct
+unit, rather than raising an exception.
diff --git a/docs/changes/stats/17221.feature.rst b/docs/changes/stats/17221.feature.rst
new file mode 100644
index 0000000000..665fd7fd20
--- /dev/null
+++ b/docs/changes/stats/17221.feature.rst
@@ -0,0 +1,2 @@
+Added a SigmaClippedStats convenience class for computing sigma-clipped
+statistics.
diff --git a/docs/changes/table/17184.api.rst b/docs/changes/table/17184.api.rst
new file mode 100644
index 0000000000..4c85e1a862
--- /dev/null
+++ b/docs/changes/table/17184.api.rst
@@ -0,0 +1,6 @@
+The default behavior of ``Table.pformat`` was changed to include all rows and columns
+instead of truncating the outputs to fit the current terminal.  The new default
+keyword arguments ``max_width=-1`` and ``max_lines=-1`` now match those in
+``Table.pformat_all``. Since the ``Table.pformat_all`` method is now redundant, it is
+pending deprecation. Similarly, the default behavior of ``Column.pformat`` was changed
+to include all rows instead of truncating the outputs to fit the current terminal.
diff --git a/docs/changes/units/17241.bugfix.rst b/docs/changes/units/17241.bugfix.rst
new file mode 100644
index 0000000000..82b8fc3f1b
--- /dev/null
+++ b/docs/changes/units/17241.bugfix.rst
@@ -0,0 +1,2 @@
+The ``"ogip"`` unit format now represents the unit angstrom as ``"angstrom"``
+instead of ``"0.1 nm"``.
diff --git a/docs/changes/units/17246.api.rst b/docs/changes/units/17246.api.rst
new file mode 100644
index 0000000000..40288788a8
--- /dev/null
+++ b/docs/changes/units/17246.api.rst
@@ -0,0 +1,4 @@
+The ``eV`` and ``rydberg`` units were moved to ``astropy.units.misc`` (from
+``astropy.units.si`` and ``astropy.units.astrophys``, respectively). They remain
+available under the root namespace ``astropy.units`` as ``astropy.units.eV``
+and ``astropy.units.rydberg``.
diff --git a/docs/changes/visualization/17217.feature.rst b/docs/changes/visualization/17217.feature.rst
new file mode 100644
index 0000000000..8d58074a9c
--- /dev/null
+++ b/docs/changes/visualization/17217.feature.rst
@@ -0,0 +1,1 @@
+Added a SimpleNorm class to create a matplotlib normalization object.
diff --git a/docs/stats/robust.rst b/docs/stats/robust.rst
index 172ac46e05..4d47bf6a1f 100644
--- a/docs/stats/robust.rst
+++ b/docs/stats/robust.rst
@@ -116,6 +116,22 @@ calculation:
      >>> sigma_clipped_stats(y, sigma=3, maxiters=10)  # doctest: +FLOAT_CMP
      (np.float64(-0.0228473012826993), np.float64(-0.02356858871405204), np.float64(0.2079616996908159))
 
+
+:class:`~astropy.stats.SigmaClippedStats` is a
+convenience class that extends the functionality of
+:func:`~astropy.stats.sigma_clipped_stats`:
+
+.. doctest-requires:: scipy
+
+     >>> from astropy.stats import SigmaClippedStats
+     >>> stats = SigmaClippedStats(y, sigma=3, maxiters=10)
+     >>> stats.mean(), stats.median(), stats.std()  # doctest: +FLOAT_CMP
+     (np.float64(-0.0228473012826993), np.float64(-0.02356858871405204), np.float64(0.2079616996908159))
+     >>> stats.mode(), stats.var(), stats.mad_std()  # doctest: +FLOAT_CMP
+     (np.float64(-0.025011163576757534), np.float64(0.043248068538293126), np.float64(0.21277510956855722))
+     >>> stats.biweight_location(), stats.biweight_scale()  # doctest: +FLOAT_CMP
+     (np.float64(-0.0183718864859565), np.float64(0.21730062377965248))
+
 :func:`~astropy.stats.sigma_clip` and
 :class:`~astropy.stats.SigmaClip` can be combined with other robust
 statistics to provide improved outlier rejection as well.
diff --git a/docs/visualization/normalization.rst b/docs/visualization/normalization.rst
index 5b6c5e5faf..ec2bae986c 100644
--- a/docs/visualization/normalization.rst
+++ b/docs/visualization/normalization.rst
@@ -39,8 +39,15 @@ The Quick Way
 -------------
 
 ``astropy`` provides a convenience
-:func:`~astropy.visualization.mpl_normalize.simple_norm` function that can be
-useful for quick interactive analysis:
+:func:`~astropy.visualization.mpl_normalize.simple_norm` function and
+a :class:`~astropy.visualization.mpl_normalize.SimpleNorm` class that
+can be useful for quick interactive analysis. These convenience tools
+create a :class:`~astropy.visualization.mpl_normalize.ImageNormalize`
+normalization object that can be used with Matplotlib's
+:meth:`~matplotlib.axes.Axes.imshow` method
+
+Here's an example using the
+:func:`~astropy.visualization.mpl_normalize.simple_norm` function:
 
 .. plot::
     :include-source:
@@ -62,12 +69,29 @@ useful for quick interactive analysis:
     im = ax.imshow(image, origin='lower', norm=norm)
     fig.colorbar(im)
 
-This convenience function combines a :class:`Stretch
-<astropy.visualization.stretch.BaseStretch>` object with an :class:`Interval
-<astropy.visualization.interval.BaseInterval>` object.
-We recommend using
-:class:`~astropy.visualization.mpl_normalize.ImageNormalize` directly
-in scripted programs instead of this convenience function.
+Here's an example using the
+:class:`~astropy.visualization.mpl_normalize.SimpleNorm` class with its
+:meth:`~astropy.visualization.mpl_normalize.SimpleNorm.imshow`
+method:
+
+.. plot::
+    :include-source:
+    :align: center
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from astropy.visualization import SimpleNorm
+
+    # Generate a test image
+    image = np.arange(65536).reshape((256, 256))
+
+    # Create an ImageNormalize object
+    snorm = SimpleNorm('sqrt', percent=98)
+
+    # Display the image
+    fig, ax = plt.subplots()
+    axim = snorm.imshow(image, ax=ax, origin='lower')
+    fig.colorbar(axim)
 
 
 The detailed way
diff --git a/docs/whatsnew/7.0.rst b/docs/whatsnew/7.0.rst
index 5348992eaf..b770ee3c1e 100644
--- a/docs/whatsnew/7.0.rst
+++ b/docs/whatsnew/7.0.rst
@@ -16,6 +16,7 @@ In particular, this release includes:
 * :ref:`whatsnew-7.0-table-frames`
 * :ref:`whatsnew-7.0-table-notebook-backends`
 * :ref:`whatsnew-7.0-table-init-columns-order`
+* :ref:`whatsnew_7_0_table_pformat_default_changes`
 * :ref:`whatsnew_7_0_quantity_to_string_formatter`
 * :ref:`whatsnew_7_0_numpy_constructors_like_quantity`
 * :ref:`whatsnew_7_0_ecsv_meta_default_dict`
@@ -26,6 +27,8 @@ In particular, this release includes:
 * :ref:`whatsnew_7_0_rgb_image_visualization_enhancement`
 * :ref:`whatsnew_7_0_lorentz2d_model`
 * :ref:`whatsnew_7_0_votable_1_5`
+* :ref:`whatsnew_7_0_simplenorm`
+* :ref:`whatsnew_7_0_sigmaclippedstats`
 
 In addition to these major changes, Astropy v7.0 includes a large number of
 smaller improvements and bug fixes, which are described in the :ref:`changelog`.
@@ -115,6 +118,21 @@ Starting with ``7.0`` the table would instead look like this::
      10   7  --
      20  35  15
 
+.. _whatsnew_7_0_table_pformat_default_changes:
+
+``Table.pformat`` is now independent of terminal dimensions
+===========================================================
+
+``Table.pformat`` and ``Column.pformat`` do not truncate their outputs according
+to terminal height and width by default any more. The new default behavior is
+intended to be less surprising.
+
+Truncating representations to fit the current terminal is still supported but
+now requires explicitly passing ``max_lines=None`` and/or ``max_width=None``.
+
+``Table.pformat_all`` is deprecated as it is now fully redundant.
+
+
 .. _whatsnew_7_0_quantity_to_string_formatter:
 
 ``Quantity.to_string`` supports ``formatter`` for formatting
@@ -355,6 +373,68 @@ now supports version 1.5 of the VOTable standard.  The main new feature is that
 At this writing, version 1.5 is a proposed standard, but it is expected to be approved as an
 official recommendation soon.
 
+.. _whatsnew_7_0_simplenorm:
+
+New ``SimpleNorm`` class
+========================
+
+A new convenience class, :class:`~astropy.visualization.SimpleNorm`,
+has been added to the ``astropy.visualization`` module.
+This class provides a simple interface to create
+a :class:`~astropy.visualization.ImageNormalize`
+normalization object that can be used with Matplotlib's
+:meth:`~matplotlib.axes.Axes.imshow` method. It also provides
+a :meth:`~astropy.visualization.SimpleNorm.imshow` method that
+wraps Matplotlib's :meth:`matplotlib.axes.Axes.imshow` method and
+automatically sets the normalization.
+
+Here's an example using the :func:`~astropy.visualization.SimpleNorm`
+function with its :meth:`~astropy.visualization.SimpleNorm.imshow`
+method:
+
+.. plot::
+    :include-source:
+
+    import numpy as np
+    import matplotlib.pyplot as plt
+    from astropy.visualization import SimpleNorm
+
+    # Generate a test image
+    image = np.arange(65536).reshape((256, 256))
+
+    # Create an ImageNormalize object
+    snorm = SimpleNorm('sqrt', percent=98)
+
+    # Display the image
+    fig, ax = plt.subplots()
+    axim = snorm.imshow(image, ax=ax, origin='lower')
+    fig.colorbar(axim)
+
+.. _whatsnew_7_0_sigmaclippedstats:
+
+New ``SigmaClippedStats`` class
+===============================
+
+A new convenience class, :class:`~astropy.stats.SigmaClippedStats`, has
+been added to the :mod:`~astropy.stats` module. This class provides a
+convenient way to compute statistics of an array with sigma-clipping. A
+simple example might be:
+
+.. doctest-skip::
+
+    >>> import numpy as np
+    >>> from astropy.stats import SigmaClippedStats
+    >>> rng = np.random.default_rng(seed=42)
+    >>> data = rng.exponential(scale=100, size=500)
+    >>> stats = SigmaClippedStats(data, sigma=3, maxiters=10)
+    >>> stats.min(), stats.max(), stats.sum()  # doctest: +FLOAT_CMP
+    (np.float64(0.8129422833034009), np.float64(255.34193997940474), np.float64(36783.895498717866))
+    >>> stats.mean(), stats.median(), stats.std()  # doctest: +FLOAT_CMP
+    (np.float64(79.61882142579624), np.float64(60.01103363578014), np.float64(65.4457063794851))
+    >>> stats.mode(), stats.var(), stats.mad_std()  # doctest: +FLOAT_CMP
+    (np.float64(20.79545805574793), np.float64(4283.1404835097765), np.float64(62.608350894722484))
+    >>> stats.biweight_location(), stats.biweight_scale()  # doctest: +FLOAT_CMP
+    (np.float64(67.98055399699436), np.float64(64.82889460022386))
 
 Full change log
 ===============


Test patch:
diff --git a/astropy/coordinates/tests/test_angles.py b/astropy/coordinates/tests/test_angles.py
index 0ee290877b..0c5f17290d 100644
--- a/astropy/coordinates/tests/test_angles.py
+++ b/astropy/coordinates/tests/test_angles.py
@@ -216,6 +216,28 @@ def test_angle_methods():
     assert a_min == 0.0 * u.degree
 
 
+def test_angle_nan_functions():
+    # Most numpy functions tested as part of the Quantity tests.
+    # But check that we drop to Quantity when appropriate; see
+    # https://github.com/astropy/astropy/pull/17221#discussion_r1813060768
+    a = Angle([0.0, 2.0, np.nan], "deg")
+    a_mean = np.nanmean(a)
+    assert type(a_mean) is Angle
+    assert a_mean == 1.0 * u.degree
+    a_std = np.nanstd(a)
+    assert type(a_std) is Angle
+    assert a_std == 1.0 * u.degree
+    a_var = np.nanvar(a)
+    assert type(a_var) is u.Quantity
+    assert a_var == 1.0 * u.degree**2
+    a_max = np.nanmax(a)
+    assert type(a_max) is Angle
+    assert a_max == 2.0 * u.degree
+    a_min = np.nanmin(a)
+    assert type(a_min) is Angle
+    assert a_min == 0.0 * u.degree
+
+
 def test_angle_convert():
     """
     Test unit conversion of Angle objects
diff --git a/astropy/io/ascii/tests/test_read.py b/astropy/io/ascii/tests/test_read.py
index 1cdff46e02..c8223085cd 100644
--- a/astropy/io/ascii/tests/test_read.py
+++ b/astropy/io/ascii/tests/test_read.py
@@ -2061,7 +2061,7 @@ def test_read_converters_simplified():
 
     converters = {"a": float, "*": [np.int64, float, bool, str]}
     t2 = Table.read(out.getvalue(), format="ascii.basic", converters=converters)
-    assert t2.pformat_all(show_dtype=True) == [
+    assert t2.pformat(show_dtype=True) == [
         "   a       b      c     d     e  ",
         "float64 float64  bool  str5 int64",
         "------- ------- ----- ----- -----",
diff --git a/astropy/stats/tests/test_sigma_clipping.py b/astropy/stats/tests/test_sigma_clipping.py
index dcc52beab6..886e0aa140 100644
--- a/astropy/stats/tests/test_sigma_clipping.py
+++ b/astropy/stats/tests/test_sigma_clipping.py
@@ -5,8 +5,14 @@ import pytest
 from numpy.testing import assert_allclose, assert_equal
 
 from astropy import units as u
+from astropy.coordinates import Angle
 from astropy.stats import mad_std
-from astropy.stats.sigma_clipping import SigmaClip, sigma_clip, sigma_clipped_stats
+from astropy.stats.sigma_clipping import (
+    SigmaClip,
+    SigmaClippedStats,
+    sigma_clip,
+    sigma_clipped_stats,
+)
 from astropy.table import MaskedColumn
 from astropy.utils.compat import COPY_IF_NEEDED
 from astropy.utils.compat.optional_deps import HAS_BOTTLENECK, HAS_SCIPY
@@ -174,6 +180,91 @@ def test_sigma_clipped_stats_masked_col():
     sigma_clipped_stats(col)
 
 
+@pytest.mark.parametrize("units", [False, True])
+@pytest.mark.parametrize("sigma", np.linspace(2, 5, 10))
+def test_sigmaclippedstats_stats(units, sigma):
+    """Test SigmaClippedStats class."""
+    data = np.ones((101, 205))
+    rng = np.random.default_rng(0)
+    idx = rng.integers(0, 100, 10)
+    data[idx, idx] = 1000
+
+    if units:
+        unit = u.m
+        data <<= unit
+    else:
+        unit = 1
+
+    stats = SigmaClippedStats(data, sigma=sigma)
+    assert stats.min() == 1.0 * unit
+    assert stats.max() == 1.0 * unit
+    assert stats.sum() == np.sum(data) - (1000 * 10) * unit
+    assert stats.mean() == 1.0 * unit
+    assert stats.median() == 1.0 * unit
+    assert stats.mode() == 1.0 * unit
+    assert stats.mode(median_factor=3.5, mean_factor=2.5) == 1.0 * unit
+    assert stats.std() == 0.0 * unit
+    assert stats.var() == 0.0 * unit**2
+    assert stats.mad_std() == 0.0 * unit
+    assert stats.biweight_location() == 1.0 * unit
+    assert stats.biweight_scale() == 0.0 * unit
+
+    # test nanvar with Angle
+    data = Angle([10, 20, 30, 40], "deg")
+    stats = SigmaClippedStats(data, sigma=5)
+    assert isinstance(stats.mean(), u.Quantity)
+    assert isinstance(stats.var(), u.Quantity)
+    assert not isinstance(stats.var(), Angle)
+    assert stats.mean() == 25 * u.deg
+    assert stats.var() == 125 * u.deg**2
+
+
+def test_sigmaclippedstats_mask():
+    data = np.ones((101, 205))
+    rng = np.random.default_rng(0)
+    idx = rng.integers(0, 100, 10)
+    data[idx, idx] = 1000
+    stats = SigmaClippedStats(data, sigma=3.2)
+    assert stats.min() == 1.0
+    assert stats.max() == 1.0
+    assert stats.sum() == np.sum(data) - 1000 * 10
+    assert stats.mean() == 1.0
+    assert stats.median() == 1.0
+    assert stats.mode() == 1.0
+    assert stats.mode(median_factor=3.5, mean_factor=2.5) == 1.0
+    assert stats.std() == 0.0
+    assert stats.var() == 0.0
+    assert stats.mad_std() == 0.0
+    assert stats.biweight_location() == 1.0
+    assert stats.biweight_scale() == 0.0
+
+    mask = np.zeros_like(data, dtype=bool)
+    mask[idx, idx] = True
+    stats3 = SigmaClippedStats(data, sigma=3.2, mask=mask)
+    assert stats3.mode() == stats.mode()
+
+    mask = np.zeros_like(data, dtype=bool)
+    stats4 = SigmaClippedStats(data, sigma=3.2, mask_value=1000)
+    assert stats4.mode() == stats3.mode()
+
+    data = np.ones((101, 205))
+    mask = np.ones_like(data, dtype=bool)
+    match = "input data is all masked"
+    with pytest.raises(ValueError, match=match):
+        SigmaClippedStats(data, sigma=3, mask=mask)
+
+
+def test_sigmaclippedstats_axis():
+    data = np.ones((101, 205))
+    stats = SigmaClippedStats(data, sigma=2.8, axis=1)
+    shape = (data.shape[0],)
+    assert stats.mean().shape == shape
+    assert stats.mode().shape == shape
+    assert stats.var().shape == shape
+    assert stats.mad_std().shape == shape
+    assert stats.biweight_location().shape == shape
+
+
 @pytest.mark.slow
 @pytest.mark.skipif(
     not HAS_BOTTLENECK,
diff --git a/astropy/table/mixins/tests/test_dask.py b/astropy/table/mixins/tests/test_dask.py
index 6e88137dcf..fda88c4e9d 100644
--- a/astropy/table/mixins/tests/test_dask.py
+++ b/astropy/table/mixins/tests/test_dask.py
@@ -39,7 +39,7 @@ class TestDaskHandler:
         assert_equal(sub["a"].compute(), np.arange(10))
 
     def test_pformat(self):
-        assert self.t.pformat_all() == [
+        assert self.t.pformat() == [
             " a ",
             "---",
             "  0",
diff --git a/astropy/table/tests/test_pprint.py b/astropy/table/tests/test_pprint.py
index ce76ab3881..9f9819345c 100644
--- a/astropy/table/tests/test_pprint.py
+++ b/astropy/table/tests/test_pprint.py
@@ -12,6 +12,7 @@ from astropy import units as u
 from astropy.io import ascii
 from astropy.table import Column, QTable, Table
 from astropy.table.table_helpers import simple_table
+from astropy.utils.exceptions import AstropyDeprecationWarning
 
 BIG_WIDE_ARR = np.arange(2000, dtype=np.float64).reshape(100, 20)
 SMALL_ARR = np.arange(18, dtype=np.int64).reshape(6, 3)
@@ -162,7 +163,7 @@ class TestPprint:
         self._setup(table_type)
         arr = np.arange(4000, dtype=np.float64).reshape(100, 40)
         with conf.set_temp("max_width", None), conf.set_temp("max_lines", None):
-            lines = table_type(arr).pformat()
+            lines = table_type(arr).pformat(max_lines=None, max_width=None)
         width, nlines = get_terminal_size()
         assert len(lines) == nlines
         for line in lines[:-1]:  # skip last "Length = .. rows" line
@@ -303,7 +304,14 @@ class TestPprint:
     def test_pformat_all(self, table_type):
         """Test that all rows are printed by default"""
         self._setup(table_type)
-        lines = self.tb.pformat_all()
+        with pytest.warns(
+            AstropyDeprecationWarning,
+            match=(
+                r"The pformat_all function is deprecated "
+                r"and may be removed in a future version\."
+            ),
+        ):
+            lines = self.tb.pformat_all()
         # +3 accounts for the three header lines in this  table
         assert len(lines) == BIG_WIDE_ARR.shape[0] + 3
 
@@ -693,7 +701,7 @@ def test_pprint_structured_with_format():
         "  1    1.23 -20.0 003 bar ",
         "  2   12.35   4.6 033 foo ",
     ]
-    assert t.pformat_all() == exp
+    assert t.pformat() == exp
 
 
 def test_pprint_nameless_col():
@@ -1032,26 +1040,26 @@ class TestColumnsShowHide:
         ]
 
         with t.pprint_exclude_names.set(["a", "c"]):
-            out = t.pformat_all()
+            out = t.pformat()
         assert out == exp
 
         with t.pprint_include_names.set(["b", "d"]):
-            out = t.pformat_all()
+            out = t.pformat()
         assert out == exp
 
         with t.pprint_exclude_names.set(["a", "c"]):
-            out = t.pformat_all()
+            out = t.pformat()
         assert out == exp
 
         with t.pprint_include_names.set(["b", "d"]):
-            out = t.pformat_all()
+            out = t.pformat()
         assert out == exp
 
         with (
             t.pprint_include_names.set(["b", "c", "d"]),
             t.pprint_exclude_names.set(["c"]),
         ):
-            out = t.pformat_all()
+            out = t.pformat()
         assert out == exp
 
     def test_output_globs(self):
@@ -1067,7 +1075,7 @@ class TestColumnsShowHide:
             "  1   1   2",
         ]
         with t.pprint_include_names.set("a*"):
-            out = t.pformat_all()
+            out = t.pformat()
         assert out == exp
 
         # Show a* but exclude a??
@@ -1077,7 +1085,7 @@ class TestColumnsShowHide:
             "  1   1",
         ]
         with t.pprint_include_names.set("a*"), t.pprint_exclude_names.set("a??"):
-            out = t.pformat_all()
+            out = t.pformat()
         assert out == exp
 
         # Exclude a??
@@ -1087,7 +1095,7 @@ class TestColumnsShowHide:
             "  1   2   3   4   1",
         ]
         with t.pprint_exclude_names.set("a??"):
-            out = t.pformat_all()
+            out = t.pformat()
         assert out == exp
 
 
@@ -1105,7 +1113,7 @@ def test_embedded_newline_tab():
         r"   a b \n c \t \n d",
         r"   x            y\n",
     ]
-    assert t.pformat_all() == exp
+    assert t.pformat() == exp
 
 
 def test_multidims_with_zero_dim():
@@ -1120,7 +1128,7 @@ def test_multidims_with_zero_dim():
         "   a             ",
         "   b             ",
     ]
-    assert t.pformat_all(show_dtype=True) == exp
+    assert t.pformat(show_dtype=True) == exp
 
 
 def test_zero_length_string():
@@ -1132,4 +1140,4 @@ def test_zero_length_string():
         "------ -----",
         "          12",
     ]
-    assert t.pformat_all(show_dtype=True) == exp
+    assert t.pformat(show_dtype=True) == exp
diff --git a/astropy/units/tests/test_format.py b/astropy/units/tests/test_format.py
index 01c8a73ae6..4b5386a355 100644
--- a/astropy/units/tests/test_format.py
+++ b/astropy/units/tests/test_format.py
@@ -529,7 +529,7 @@ def test_latex_inline_scale():
         ("unicode", "erg Å⁻¹ s⁻¹ cm⁻²", "10000000 kg m⁻¹ s⁻³"),
         (">25s", "   erg / (Angstrom s cm2)", "        1e+07 kg / (m s3)"),
         ("cds", "erg.Angstrom-1.s-1.cm-2", "10000000kg.m-1.s-3"),
-        ("ogip", "10 erg / (nm s cm**2)", "1e+07 kg / (m s**3)"),
+        ("ogip", "erg / (angstrom s cm**2)", "1e+07 kg / (m s**3)"),
         ("fits", "erg Angstrom-1 s-1 cm-2", "10**7 kg m-1 s-3"),
         ("vounit", "erg.Angstrom**-1.s**-1.cm**-2", "10000000kg.m**-1.s**-3"),
         # TODO: make fits and vounit less awful!
diff --git a/astropy/visualization/tests/test_norm.py b/astropy/visualization/tests/test_norm.py
index 0e703680a5..093540aaf1 100644
--- a/astropy/visualization/tests/test_norm.py
+++ b/astropy/visualization/tests/test_norm.py
@@ -7,7 +7,12 @@ from numpy.testing import assert_allclose, assert_array_equal, assert_equal
 
 from astropy.utils.compat.optional_deps import HAS_MATPLOTLIB, HAS_PLT
 from astropy.visualization.interval import ManualInterval, PercentileInterval
-from astropy.visualization.mpl_normalize import ImageNormalize, imshow_norm, simple_norm
+from astropy.visualization.mpl_normalize import (
+    ImageNormalize,
+    SimpleNorm,
+    imshow_norm,
+    simple_norm,
+)
 from astropy.visualization.stretch import LogStretch, PowerStretch, SqrtStretch
 
 DATA = np.linspace(0.0, 15.0, 6)
@@ -281,31 +286,55 @@ class TestImageScaling:
             simple_norm(DATA2, stretch="invalid")
 
 
+@pytest.mark.skipif(not HAS_PLT, reason="requires matplotlib.pyplot")
+@pytest.mark.parametrize("stretch", ["linear", "sqrt", "power", "log", "asinh", "sinh"])
+def test_simplenorm(stretch):
+    data = np.arange(25).reshape((5, 5))
+    snorm = SimpleNorm(stretch, percent=99)
+    norm = snorm(data)
+    assert isinstance(norm, ImageNormalize)
+    assert_allclose(norm(data), simple_norm(data, stretch, percent=99)(data))
+
+
+@pytest.mark.skipif(not HAS_PLT, reason="requires matplotlib.pyplot")
+def test_simplenorm_imshow():
+    import matplotlib.pyplot as plt
+    from matplotlib.image import AxesImage
+
+    data = np.arange(25).reshape((5, 5))
+    fig, ax = plt.subplots()
+    snorm = SimpleNorm("sqrt", percent=99)
+    axim = snorm.imshow(data, ax=ax)
+    assert isinstance(axim, AxesImage)
+    keys = ("vmin", "vmax", "stretch", "clip", "invalid")
+    for key in keys:
+        assert getattr(axim.norm, key) == getattr(snorm(data), key)
+
+    fig.clear()
+    axim = snorm.imshow(data, ax=None)
+
+    with pytest.raises(ValueError):
+        snorm.imshow(data, ax=ax, norm=ImageNormalize())
+
+
 @pytest.mark.skipif(not HAS_PLT, reason="requires matplotlib.pyplot")
 def test_imshow_norm():
     import matplotlib.pyplot as plt
 
     image = np.random.randn(10, 10)
 
-    plt.clf()
-    ax = plt.subplot(label="test_imshow_norm")
+    fig, ax = plt.subplots(label="test_imshow_norm")
     imshow_norm(image, ax=ax)
 
-    with pytest.raises(ValueError):
-        # X and data are the same, can't give both
-        imshow_norm(image, X=image, ax=ax)
-
     with pytest.raises(ValueError):
         # illegal to manually pass in normalization since that defeats the point
         imshow_norm(image, ax=ax, norm=ImageNormalize())
 
-    plt.clf()
+    fig.clear()
     imshow_norm(image, ax=ax, vmin=0, vmax=1)
 
     # make sure the pyplot version works
-    plt.clf()
+    fig.clear()
     imres, norm = imshow_norm(image, ax=None)
 
     assert isinstance(norm, ImageNormalize)
-
-    plt.close("all")