diff --git a/.github/workflows/ci_tests.yml b/.github/workflows/ci_tests.yml
index 2fbbd2ce..57f4346f 100644
--- a/.github/workflows/ci_tests.yml
+++ b/.github/workflows/ci_tests.yml
@@ -56,7 +56,7 @@ jobs:
           - name: 'ubuntu-py312-bottleneck'
             os: ubuntu-latest
             python: '3.12'
-            tox_env: 'py312-test-alldeps-numpy126-bottleneck'
+            tox_env: 'py312-test-alldeps-numpy200-bottleneck'
 
           - name: 'macos-py312-dask'
             os: macos-latest
diff --git a/CHANGES.rst b/CHANGES.rst
index 546cb592..205f78d2 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -22,6 +22,10 @@ Other Changes and Additions
   without having to name the interpreter. [#986]
 - ``combine`` no longer accepts an array as its ``array_package`` argument;
   pass the array namespace or module instead, as for ``Combiner``. [#997]
+- The minimum supported numpy is now 2.0 (``normalize_axis_tuple`` is
+  imported from its numpy 2 location with no 1.26 fallback), which raises
+  the minimum astropy to 6.1, reproject to 0.14 and astroscrappy to 1.2,
+  the oldest releases that work with numpy 2. [#1000]
 
 Bug Fixes
 ^^^^^^^^^
@@ -121,6 +125,15 @@ Bug Fixes
   dtype instead of passing it through unchanged. [#999]
 - Size the overscan model fit in ``subtract_overscan`` with ``shape``
   instead of ``len``, which the array API standard does not provide. [#999]
+- Compute ``sigma_func`` (the ``median_combine`` uncertainty) with a median
+  absolute deviation written purely in terms of the array API when the array
+  namespace is not numpy, instead of converting the data to numpy through
+  ``astropy.stats.median_absolute_deviation``. Numpy input still uses
+  astropy. [#1000]
+- ``sigma_func`` now always excludes the masked pixels of a masked
+  ``CCDData``. Previously the mask was only honored on numpy for small
+  arrays with a single integer ``axis`` and ``ignore_nan=True``, and, when
+  bottleneck is installed, not for float64 data at all. [#1000]
 
 2.5.1 (2025-07-05)
 ------------------
diff --git a/ccdproc/_nanfuncs.py b/ccdproc/_nanfuncs.py
index 75afb378..3f66b2da 100644
--- a/ccdproc/_nanfuncs.py
+++ b/ccdproc/_nanfuncs.py
@@ -26,6 +26,35 @@
 __all__ = ["median", "nanmean", "nanmedian", "nanstd", "nansum"]
 
 
+def _promote_to_real(x, xp, device):
+    """
+    Promote integer and boolean ``x`` to the namespace's default real
+    floating dtype; a real floating ``x``, including float32, passes
+    through unchanged.
+
+    Parameters
+    ----------
+    x : array
+        Input array.
+    xp : array namespace
+        Namespace to use.
+    device : device
+        Device on which to resolve the default real floating dtype.
+
+    Returns
+    -------
+    array
+        ``x``, promoted if necessary.
+    """
+    if xp.isdtype(x.dtype, "real floating"):
+        return x
+    # Promote to the namespace's default real dtype rather than hardcoding
+    # float64: jax without JAX_ENABLE_X64 has no float64 and warns when one
+    # is requested, which pytest's filterwarnings turns into an error.
+    info = xp.__array_namespace_info__()
+    return xp.astype(x, info.default_dtypes(device=device)["real floating"])
+
+
 def _setup(x, axis, xp):
     """
     Validate ``axis``, resolve the namespace and device, promote to float.
@@ -83,13 +112,7 @@ def _setup(x, axis, xp):
     axis = axis % ndim
 
     device = array_api_compat.device(x)
-
-    if not xp.isdtype(x.dtype, "real floating"):
-        # Promote to the namespace's default real dtype rather than hardcoding
-        # float64: jax without JAX_ENABLE_X64 has no float64 and warns when one
-        # is requested, which pytest's filterwarnings turns into an error.
-        info = xp.__array_namespace_info__()
-        x = xp.astype(x, info.default_dtypes(device=device)["real floating"])
+    x = _promote_to_real(x, xp, device)
 
     return x, axis, xp, device
 
diff --git a/ccdproc/core.py b/ccdproc/core.py
index 9b70e154..35f72c21 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -21,6 +21,7 @@
 from astropy.utils import deprecated, deprecated_renamed_argument
 from astropy.wcs.utils import proj_plane_pixel_area
 from numpy import mgrid as np_mgrid
+from numpy.lib.array_utils import normalize_axis_tuple
 from numpy.ma import nomask as np_ma_nomask
 from packaging import version as pkgversion
 from scipy import ndimage
@@ -29,7 +30,9 @@
     _unwrap_ccddata_for_array_api,
     _wrap_ccddata_for_array_api,
 )
+from ._nanfuncs import _promote_to_real
 from ._nanfuncs import median as _nanfuncs_median
+from ._nanfuncs import nanmedian as _nanfuncs_nanmedian
 from .log_meta import log_to_metadata
 from .utils.slices import slice_from_string
 
@@ -286,6 +289,92 @@ def _median_fallback(array, axis, xp=None):
         return _nanfuncs_median(array, axis=axis, xp=xp)
 
 
+def _mad_fallback(data, axis, ignore_nan, xp=None):
+    """
+    Median absolute deviation written purely in terms of the array API.
+
+    This is the non-numpy branch of `sigma_func`;
+    `astropy.stats.median_absolute_deviation` is numpy-only, so any other
+    namespace computes the deviation here, on the device of ``data``.
+
+    Parameters
+    ----------
+    data : array
+        Array whose median absolute deviation is to be calculated. Integer
+        and boolean input is promoted to the namespace's default real
+        floating dtype.
+
+    axis : int, tuple of int, list of int or None
+        Axis or axes along which the deviation is computed. ``None``
+        flattens ``data`` first; a tuple or list reduces over all the
+        listed axes. Negative values count from the last axis.
+
+    ignore_nan : bool
+        If `True`, NaNs are ignored; otherwise a NaN anywhere in a reduced
+        slice makes that slice's result NaN, as `numpy.median` does.
+
+    xp : array namespace, optional
+        Namespace to use. If not provided, it is determined from ``data``.
+
+    Returns
+    -------
+    mad : array
+        Median absolute deviation of ``data`` along ``axis``, with the
+        reduced axes removed (0-d when ``axis`` is `None`), in the namespace
+        and on the device of ``data``.
+
+    Raises
+    ------
+    ValueError
+        If a tuple or list ``axis`` is out of bounds for ``data`` or lists
+        an axis more than once (including via a negative alias).
+    TypeError
+        If a tuple or list ``axis`` contains a bool or a non-integer entry.
+
+    Notes
+    -----
+    The namespace's own ``nanmedian``/``median`` is used when it has one;
+    otherwise the sort-based fallbacks from `ccdproc._nanfuncs` are used
+    instead, which cost O(n log n) along the reduced axes rather than the
+    O(n) of a native median, and yield NaN for an entirely NaN slice
+    without the ``RuntimeWarning`` numpy emits.
+    """
+    xp = xp or array_api_compat.array_namespace(data)
+    device = array_api_compat.device(data)
+
+    # The medians below promote internally, but ``data - center`` below
+    # needs the promoted input too: subtracting a float from an integer array
+    # is not allowed by the standard and raises on array-api-strict.
+    data = _promote_to_real(data, xp, device)
+
+    if axis is None:
+        data = xp.reshape(data, (-1,))
+        axis = 0
+    elif isinstance(axis, (tuple, list)):
+        if any(isinstance(ax, bool) for ax in axis):
+            raise TypeError("axis entries must be integers, not bool")
+        ndim = data.ndim
+        axes = normalize_axis_tuple(axis, ndim)
+        # Move the reduced axes to the end and merge them into one, so that
+        # the single-axis medians below reduce over all of them at once.
+        kept = [ax for ax in range(ndim) if ax not in axes]
+        data = xp.permute_dims(data, tuple(kept + list(axes)))
+        data = xp.reshape(
+            data, tuple(data.shape[ax] for ax in range(len(kept))) + (-1,)
+        )
+        axis = -1
+
+    def med(d, ax):
+        try:
+            return (xp.nanmedian if ignore_nan else xp.median)(d, axis=ax)
+        except AttributeError:
+            fallback = _nanfuncs_nanmedian if ignore_nan else _nanfuncs_median
+            return fallback(d, axis=ax, xp=xp)
+
+    center = med(data, axis)
+    return med(xp.abs(data - xp.expand_dims(center, axis=axis)), axis)
+
+
 @log_to_metadata
 def ccd_process(
     ccd,
@@ -1385,7 +1474,7 @@ def sigma_func(arr, axis=None, ignore_nan=False):
 
     Parameters
     ----------
-    arr : `~astropy.nddata.CCDData` or `numpy.ndarray`
+    arr : `~astropy.nddata.CCDData` or array
         Array whose deviation is to be calculated.
 
     axis : int, tuple of ints or None, optional
@@ -1395,20 +1484,57 @@ def sigma_func(arr, axis=None, ignore_nan=False):
         it counts from the last to the first axis.
         Default is ``None``.
 
+    ignore_nan : bool, optional
+        If `True`, NaNs are ignored when computing the medians; otherwise a
+        NaN in a reduced slice makes that slice's result NaN.
+        Default is ``False``.
+
     Returns
     -------
-    uncertainty : float
-        uncertainty of array estimated from median absolute deviation.
+    uncertainty : array
+        Uncertainty of the array estimated from the median absolute
+        deviation, in the array namespace and on the device of the input.
+        It is 0-d when ``axis`` is ``None``.
+
+    Notes
+    -----
+    For numpy input the deviation is computed by
+    `astropy.stats.median_absolute_deviation`, exactly as before. For every
+    other array namespace it is computed by a fallback written purely in
+    terms of the array API standard, on the device of the input, using the
+    namespace's own ``median``/``nanmedian`` when it provides one and
+    otherwise a sort-based implementation that costs O(n log n) along the
+    reduced axes rather than O(n); it promotes integer and boolean input to
+    the namespace's default real floating dtype, and it yields NaN silently
+    for slices that are entirely NaN, where numpy would warn.
+
+    The masked pixels of a masked `~astropy.nddata.CCDData` are excluded
+    from the statistics on every backend, which implies ``ignore_nan``; a
+    slice that is entirely masked gives NaN.
     """
     if isinstance(arr, CCDData):
-        xp = array_api_compat.array_namespace(arr.data)
+        data = arr.data
+        mask = arr.mask
     else:
-        xp = array_api_compat.array_namespace(arr)
-
-    return xp.asarray(
-        stats.median_absolute_deviation(arr, axis=axis, ignore_nan=ignore_nan)
-        * 1.482602218505602
-    )
+        data = arr
+        mask = None
+    xp = array_api_compat.array_namespace(data)
+
+    if mask is not None:
+        device = array_api_compat.device(data)
+        # The promotion has to happen before the ``where`` below: on
+        # array-api-strict, ``xp.where(mask, float_nan, int_data)`` raises
+        # rather than silently promoting the way numpy does.
+        data = _promote_to_real(data, xp, device)
+        nan = xp.asarray(xp.nan, dtype=data.dtype, device=device)
+        data = xp.where(xp.asarray(mask, device=device), nan, data)
+        ignore_nan = True
+
+    if array_api_compat.is_numpy_namespace(xp):
+        result = stats.median_absolute_deviation(data, axis=axis, ignore_nan=ignore_nan)
+        return xp.asarray(result * 1.482602218505602)
+
+    return _mad_fallback(data, axis, ignore_nan, xp=xp) * 1.482602218505602
 
 
 def setbox(x, y, mbox, xmax, ymax):
diff --git a/ccdproc/tests/array_escape_baseline.txt b/ccdproc/tests/array_escape_baseline.txt
index 6b82c205..634a21d9 100644
--- a/ccdproc/tests/array_escape_baseline.txt
+++ b/ccdproc/tests/array_escape_baseline.txt
@@ -23,5 +23,4 @@ core.py      block_reduce                 numpy.asanyarray  BOUNDARY: astropy.nd
 core.py      block_replicate              numpy.asanyarray  BOUNDARY: astropy.nddata.block_replicate is numpy-only
 core.py      cosmicray_median             numpy.asarray     BOUNDARY: scipy.ndimage median_filter/maximum_filter are numpy-only
 core.py      create_deviation             numpy.asarray     BOUNDARY: astropy CCDData.copy()/StdDevUncertainty are numpy-backed
-core.py      sigma_func                   numpy.asanyarray  BOUNDARY: astropy.stats.median_absolute_deviation is numpy-only
 core.py      subtract_overscan            numpy.asanyarray  BOUNDARY: astropy.modeling fit/evaluation (model= path) is numpy-only
diff --git a/ccdproc/tests/test_ccdproc.py b/ccdproc/tests/test_ccdproc.py
index 8f51a54e..27ef1ab5 100644
--- a/ccdproc/tests/test_ccdproc.py
+++ b/ccdproc/tests/test_ccdproc.py
@@ -16,17 +16,24 @@
     StdDevUncertainty,
     VarianceUncertainty,
 )
+from astropy.stats import mad_std, median_absolute_deviation
 from astropy.units.quantity import Quantity
 from astropy.utils.exceptions import AstropyUserWarning
 from astropy.wcs import WCS
 from numpy import array as np_array
+from numpy import asarray as np_asarray
+from numpy import float32 as np_float32
 from numpy import mgrid as np_mgrid
+from numpy import nan as np_nan
+from numpy import nan_to_num as np_nan_to_num
+from numpy import ones as np_ones
 from numpy import random as np_random
 
 from ccdproc.conftest import testing_array_device as xp_device
 from ccdproc.conftest import testing_array_library as xp
 from ccdproc.core import (
     Keyword,
+    _mad_fallback,
     _median_fallback,
     ccd_process,
     cosmicray_lacosmic,
@@ -396,6 +403,167 @@ def __getattr__(self, name):
     assert xp.all(xpx.isclose(result, expected))
 
 
+_MAD_RNG = np_random.default_rng(929)
+_MAD_3D = _MAD_RNG.normal(size=(5, 4, 3))
+_MAD_3D[[0, 1, 2, 4], [1, 2, 0, 3], [0, 2, 1, 1]] = np_nan  # scattered NaNs
+
+_MAD_CASES = [
+    # axis=None (background_deviation_box), an int (median_combine), a tuple
+    # with a negative, unsorted entry, and every axis.
+    *[(_MAD_3D, axis) for axis in (None, 0, (-1, 0), (0, 1, 2))],
+    (np_array([[1, 4], [2, 3], [5, 6], [4, 1]]), 0),  # int -> default float
+    (np_array([[1.0, 2.0, 3.0], [4.0, 5.0, 7.0]], dtype=np_float32), 1),  # float32 kept
+]
+
+
+# _mad_fallback agrees with astropy.stats.median_absolute_deviation on every
+# backend for each axis form sigma_func's callers use, and handles dtype
+# like the _nanfuncs medians (int promoted, float32 kept). This is a
+# differential check against astropy over the axis forms, not a coverage
+# device.
+#
+# The reference gets a tuple axis with its negative entries normalised:
+# astropy's bottleneck dispatch (astropy.stats.nanfunctions.
+# _move_tuple_axes_last) transposes with the tuple as given and raises
+# on a negative entry. The fallback still receives the tuple as written.
+@pytest.mark.filterwarnings("ignore:All-NaN slice encountered:RuntimeWarning")
+@pytest.mark.parametrize("ignore_nan", [False, True])
+@pytest.mark.parametrize(("data", "axis"), _MAD_CASES)
+def test_mad_fallback_matches_astropy(data, axis, ignore_nan):
+    reference_axis = axis
+    if isinstance(axis, tuple):
+        reference_axis = tuple(ax % data.ndim for ax in axis)
+    expected_np = np_asarray(
+        median_absolute_deviation(data, axis=reference_axis, ignore_nan=ignore_nan)
+    )
+    expected = xp.asarray(expected_np, device=xp_device)
+
+    result = _mad_fallback(xp.asarray(data, device=xp_device), axis, ignore_nan)
+
+    assert result.shape == expected.shape
+    assert result.dtype == expected.dtype
+    assert xp.all(xpx.isclose(result, expected, equal_nan=True))
+
+
+# Duplicate (including a negative alias), out-of-bounds and bool entries in
+# a tuple axis raise instead of silently reducing the wrong axes.
+@pytest.mark.parametrize(
+    ("axis", "error", "match"),
+    [
+        pytest.param((0, 0), ValueError, "repeated axis", id="duplicate"),
+        # -3 is axis 0 of a 3-D array, so this is a duplicate too.
+        pytest.param(
+            (0, -3), ValueError, "repeated axis", id="duplicate-negative-alias"
+        ),
+        pytest.param((0, 3), ValueError, "out of bounds", id="out-of-bounds"),
+        # normalize_axis_tuple would silently treat True as 1.
+        pytest.param((0, True), TypeError, "not bool", id="bool-entry"),
+    ],
+)
+def test_mad_fallback_rejects_bad_axis_tuple(axis, error, match):
+    data = xp.asarray(_MAD_3D, device=xp_device)
+    with pytest.raises(error, match=match):
+        _mad_fallback(data, axis, True)
+
+
+# The except branch of _mad_fallback's med() only runs naturally on
+# namespaces without a native nanmedian/median (array-api-strict is the only
+# such backend in CI, and it does not report coverage), so hide the native
+# functions behind a proxy namespace that otherwise delegates to the backend
+# under test, as test_median_fallback_without_native_median does.
+@pytest.mark.parametrize("ignore_nan", [False, True])
+def test_mad_fallback_without_native_medians(ignore_nan):
+    class _NamespaceWithoutMedians(types.ModuleType):
+        def __getattr__(self, name):
+            if name in ("median", "nanmedian"):
+                raise AttributeError(name)
+            return getattr(xp, name)
+
+    proxy = _NamespaceWithoutMedians("xp_without_medians")
+    data = xp.asarray(_MAD_3D, device=xp_device)
+
+    result = _mad_fallback(data, 0, ignore_nan, xp=proxy)
+
+    expected = _mad_fallback(data, 0, ignore_nan)
+    assert xp.all(xpx.isclose(result, expected, equal_nan=True))
+
+
+# The public entry point matches astropy.stats.mad_std, stays in the
+# input's namespace and device, and for axis=None gives a 0-d result that
+# converts to float (what background_deviation_box relies on).
+@pytest.mark.parametrize("force_fallback", [False, True])
+@pytest.mark.parametrize(
+    ("data", "axis", "ignore_nan"),
+    [
+        pytest.param(np_nan_to_num(_MAD_3D), None, False, id="all-no_nan"),
+        pytest.param(_MAD_3D, 0, True, id="axis0-ignore_nan"),
+        pytest.param(_MAD_3D, (0, 1), True, id="axes01-ignore_nan"),
+    ],
+)
+def test_sigma_func_matches_mad_std(
+    data, axis, ignore_nan, force_fallback, monkeypatch
+):
+    if force_fallback:
+        # On numpy sigma_func hands the data to astropy directly; every
+        # other backend already takes the fallback branch, and none of them
+        # reports coverage, so make numpy take it too.
+        monkeypatch.setattr(array_api_compat, "is_numpy_namespace", lambda _xp: False)
+
+    expected_np = np_asarray(mad_std(data, axis=axis, ignore_nan=ignore_nan))
+    expected = xp.asarray(expected_np, device=xp_device)
+
+    arr = xp.asarray(data, device=xp_device)
+    result = sigma_func(arr, axis=axis, ignore_nan=ignore_nan)
+
+    assert result.shape == expected.shape
+    assert xp.all(xpx.isclose(result, expected, equal_nan=True))
+    assert array_api_compat.array_namespace(result) is array_api_compat.array_namespace(
+        arr
+    )
+    if xp_device is not None:
+        assert array_api_compat.device(result) == xp_device
+    if axis is None:
+        # background_deviation_box relies on the 0-d result converting to float.
+        assert result.shape == ()
+        assert float(result) == pytest.approx(float(expected_np))
+
+
+# Masked CCDData pixels are excluded on every backend by the same mask->NaN
+# substitution in sigma_func (which forces ignore_nan on); an all-masked
+# slice gives NaN, not 0.0. The ignore mark: the NaN-filled data reach
+# numpy.nanmedian, which warns on the entirely masked column and on the
+# all-masked input.
+@pytest.mark.filterwarnings("ignore:All-NaN slice encountered:RuntimeWarning")
+@pytest.mark.parametrize(
+    ("axis", "ignore_nan", "all_masked"),
+    [
+        # axis=0 with ignore_nan=True is the form median_combine uses.
+        pytest.param(0, True, False, id="axis0"),
+        pytest.param(None, False, False, id="axis-none"),
+        pytest.param(None, False, True, id="all-masked-axis-none"),
+    ],
+)
+def test_sigma_func_ccddata_mask_is_honored(axis, ignore_nan, all_masked):
+    ccd = ccd_data_func()
+    if all_masked:
+        mask_np = np_ones(ccd.shape, dtype=bool)
+    else:
+        mask_np = RNG(929).random(ccd.shape) > 0.7
+        mask_np[:, 3] = True  # an entirely masked column gives NaN along axis 0
+    mask = xp.asarray(mask_np, device=xp_device)
+    # TODO: Set .mask instead of ._mask when CCDData is array-api compliant
+    ccd._mask = mask
+    nan = xp.asarray(np_nan, dtype=ccd.data.dtype, device=xp_device)
+    nanned = xp.where(mask, nan, ccd.data)
+
+    expected = sigma_func(nanned, axis=axis, ignore_nan=True)
+    result = sigma_func(ccd, axis=axis, ignore_nan=ignore_nan)
+    assert result.shape == expected.shape
+    assert xp.all(xpx.isclose(result, expected, equal_nan=True))
+    if all_masked:
+        assert bool(xp.isnan(result))
+
+
 def test_trim_image_fits_section_requires_string():
     ccd_data = ccd_data_func()
     with pytest.raises(TypeError):
diff --git a/ccdproc/tests/test_combiner.py b/ccdproc/tests/test_combiner.py
index 5ff00b6b..002b20ce 100644
--- a/ccdproc/tests/test_combiner.py
+++ b/ccdproc/tests/test_combiner.py
@@ -720,7 +720,7 @@ def test_combiner_with_scaling():
     # Scale by a float
     avg_ccd = combiner.average_combine(scale_to=2.0)
     expected_avg = 2 * xp.mean(
-        xp.asarray((ccd_data.data, ccd_data_lower.data, ccd_data_higher.data))
+        xp.stack((ccd_data.data, ccd_data_lower.data, ccd_data_higher.data))
     )
     assert xp.all(xpx.isclose(xp.mean(avg_ccd.data), expected_avg))
     assert avg_ccd.shape == ccd_data.shape
diff --git a/docs/array_api.rst b/docs/array_api.rst
index 63e8bdac..43160157 100644
--- a/docs/array_api.rst
+++ b/docs/array_api.rst
@@ -149,14 +149,30 @@ effect only after they are merged.
 What limitations should I be aware of?
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-+ The ``median`` function is not part of the array API, but most array libraries
-  do provide a NaN-aware ``nanmedian``. When combining images, `ccdproc`_ uses
-  ``nanmedian`` from `bottleneck`_ for `numpy`_ arrays if `bottleneck`_ is
-  installed; otherwise it uses the ``nanmedian`` of the selected array library.
-  If that library has no ``nanmedian`` (for example ``array-api-strict``),
-  `ccdproc`_ falls back to a sort-based implementation written purely in terms
-  of the array API standard. That fallback is correct but slower
-  (O(n log n) along the combination axis rather than O(n)).
++ The NaN-aware reductions ``nansum``, ``nanmean``, ``nanstd`` and
+  ``nanmedian`` are not part of the array API, but most array libraries do
+  provide them. When combining images, `ccdproc`_ uses the versions from
+  `bottleneck`_ for `numpy`_ arrays if `bottleneck`_ is installed; otherwise
+  it uses the ones the selected array library provides. If that library has
+  none (for example ``array-api-strict``), `ccdproc`_ falls back to
+  implementations written purely in terms of the array API standard. They
+  are correct but slower: the sum, mean and standard deviation take a few
+  extra passes over the data, and the median is sort-based (O(n log n)
+  along the combination axis rather than O(n)). The fallbacks promote
+  integer and boolean input to the library's default real floating dtype
+  and yield NaN silently for slices that are entirely NaN.
++ ``median`` is not part of the array API either. ``subtract_overscan``
+  uses the ``median`` of the selected array library when there is one, and
+  otherwise the same sort-based fallback, propagating NaNs as
+  ``numpy.median`` does.
++ ``sigma_func``, the default uncertainty estimate of ``median_combine``,
+  is a median absolute deviation. For `numpy`_ arrays it calls
+  ``astropy.stats.median_absolute_deviation``, which is numpy-only; for
+  every other array library it uses a version written purely in terms of
+  the array API standard, preferring the library's own ``median``/
+  ``nanmedian`` and falling back to the sort-based medians above only when
+  the library has neither, so the extra sort cost and the silence on
+  all-NaN slices apply only there.
 
 Which array library should I use?
 ---------------------------------
diff --git a/pyproject.toml b/pyproject.toml
index 7f5a5ba4..189d7eb2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -17,10 +17,10 @@ authors = [
 dependencies = [
     "array_api_compat>=1.12.0",
     "array_api_extra>=0.7.0",
-    "astropy>=6.0.1",
-    "astroscrappy>=1.1.0",
-    "numpy>=1.26",
-    "reproject>=0.9.1",
+    "astropy>=6.1",
+    "astroscrappy>=1.2",
+    "numpy>=2.0",
+    "reproject>=0.14",
     "scipy",
 ]
 
diff --git a/tox.ini b/tox.ini
index 7b66af19..9fb76a3c 100644
--- a/tox.ini
+++ b/tox.ini
@@ -52,7 +52,6 @@ description =
     devdeps: with the latest developer version of key dependencies
     oldestdeps: with the oldest supported version of key dependencies
     cov: and test coverage
-    numpy126: with numpy 1.26.*
     numpy200: with numpy 2.0.*
     numpy210: with numpy 2.1.*
     bottleneck: with bottleneck
@@ -66,13 +65,9 @@ description =
 deps =
     cov: coverage
 
-    numpy126: numpy==1.26.*  # currently oldest support numpy version
-    numpy200: numpy==2.0.*
+    numpy200: numpy==2.0.*  # currently oldest supported numpy version
     numpy210: numpy==2.1.*
 
-    astroscrappy11: astroscrappy==1.1.*
-    astroscrappy11: numpy<2.0
-
     bottleneck: bottleneck>=1.3.2
 
     devdeps: astropy>=0.0.dev0
@@ -80,9 +75,9 @@ deps =
 
     # Remember to transfer any changes here to setup.cfg also. Only listing
     # packages which are constrained in the setup.cfg
-    oldestdeps: numpy==1.26.*
-    oldestdeps: astropy==6.0.*
-    oldestdeps: reproject==0.9.1
+    oldestdeps: numpy==2.0.*
+    oldestdeps: astropy==6.1.*
+    oldestdeps: reproject==0.14.0
 
     dask: dask
     # The strict env intentionally installs only the test extras plus
@@ -97,7 +92,7 @@ commands =
     cov: pytest --pyargs ccdproc {toxinidir}/docs --cov ccdproc --cov-config={toxinidir}/pyproject.toml {posargs}
     cov: coverage xml -o {toxinidir}/coverage.xml
     # install astroscrappy after numpy
-    oldestdeps: python -m pip install astroscrappy==1.1.0
+    oldestdeps: python -m pip install astroscrappy==1.2.0
     # Do not care about warnings on the oldest builds
     oldestdeps: pytest --pyargs ccdproc {toxinidir}/docs -W ignore {posargs}