diff --git a/CHANGES.rst b/CHANGES.rst
index b8baec1f..7c867a93 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -87,6 +87,10 @@ Bug Fixes
 - Fix ``background_deviation_box`` discarding the result of the functional
   array update, which left the output at the global standard deviation on
   immutable array-API backends such as JAX. [#963]
+- Add a ``median`` fallback written purely in terms of the array API
+  standard, built on the existing NaN-aware ``nanmedian``, and use it in
+  ``subtract_overscan`` when the selected array namespace has no ``median``,
+  instead of raising ``AttributeError``. [#989]
 
 2.5.1 (2025-07-05)
 ------------------
diff --git a/ccdproc/_nanfuncs.py b/ccdproc/_nanfuncs.py
index ad3e341d..75afb378 100644
--- a/ccdproc/_nanfuncs.py
+++ b/ccdproc/_nanfuncs.py
@@ -1,26 +1,29 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 """
-NaN-aware sum/mean/standard deviation/median written only in terms of the
-array API.
+NaN-aware sum/mean/standard deviation/median, and a NaN-propagating median,
+written only in terms of the array API.
 
 ``nansum``/``nanmean``/``nanstd``/``nanmedian`` are not part of the array
 API standard, so this module provides fallbacks that work on any conforming
 namespace (``array-api-strict``, ``jax``, ``dask``, ``numpy``, ...). They
 are used by `ccdproc.combiner.Combiner` when the selected namespace does
-not provide the native versions.
+not provide the native versions. ``median`` is not part of the standard
+either, and this module also provides a fallback for it, built on
+``nanmedian``, used by `ccdproc.core.subtract_overscan` when the selected
+namespace does not provide the native version.
 
-All four functions promote integer and boolean input to the namespace's
+All five functions promote integer and boolean input to the namespace's
 default real floating dtype, which is where they part company with
 ``numpy.nansum``: numpy preserves an integer dtype, these do not. Every
 caller in `ccdproc` combines floating point image data, and the promotion
-keeps the four functions consistent with each other.
+keeps the five functions consistent with each other.
 """
 
 import operator
 
 import array_api_compat
 
-__all__ = ["nanmean", "nanmedian", "nanstd", "nansum"]
+__all__ = ["median", "nanmean", "nanmedian", "nanstd", "nansum"]
 
 
 def _setup(x, axis, xp):
@@ -347,3 +350,41 @@ def nanmedian(x, /, *, axis=0, xp=None):
     # makes an all-NaN slice yield NaN. Do not remove it as redundant.
     nan = xp.asarray(xp.nan, dtype=s.dtype, device=device)
     return xp.where(xp.squeeze(n, axis=axis) == 0, nan, result)
+
+
+def median(x, /, *, axis=0, xp=None):
+    """
+    Median along an axis, using only array-API functions.
+
+    Parameters
+    ----------
+    x : array
+        Input array. Integer and boolean inputs are promoted to the
+        namespace's default real floating dtype.
+    axis : int, optional
+        Axis along which to compute the median. Default is 0. Booleans,
+        ``None`` and tuples of axes are not supported; numpy integer
+        scalars are accepted.
+    xp : array namespace, optional
+        Namespace to use. Defaults to ``array_api_compat.array_namespace(x)``.
+
+    Returns
+    -------
+    array
+        Median of ``x`` along ``axis``, with that axis removed. Slices that
+        contain any NaN yield NaN, matching `numpy.median`; this is the
+        difference from `nanmedian`, which ignores NaNs entirely.
+
+    Notes
+    -----
+    On input with no NaNs this is exactly `nanmedian` -- the same sorting
+    and index-picking algorithm is used, so the values agree bit for bit.
+    The two differ only in how a NaN in the reduced slice is handled: this
+    function propagates it to the result, matching `numpy.median`, while
+    `nanmedian` ignores it. That NaN-propagating behaviour is restored here
+    with a final `where` over whether any NaN is present along ``axis``,
+    since `nanmedian` alone would silently drop NaNs instead.
+    """
+    x, axis, xp, device = _setup(x, axis, xp)
+    nan = xp.asarray(xp.nan, dtype=x.dtype, device=device)
+    return xp.where(xp.any(xp.isnan(x), axis=axis), nan, nanmedian(x, axis=axis, xp=xp))
diff --git a/ccdproc/core.py b/ccdproc/core.py
index 81fbc796..72f28a7b 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -28,6 +28,7 @@
     _unwrap_ccddata_for_array_api,
     _wrap_ccddata_for_array_api,
 )
+from ._nanfuncs import median as _nanfuncs_median
 from .log_meta import log_to_metadata
 from .utils.slices import slice_from_string
 
@@ -152,6 +153,39 @@ def _percentile_fallback(array, percentiles, xp=None):
     return sorted_array[indexes]
 
 
+def _median_fallback(array, axis, xp=None):
+    """
+    Try calculating the median using the namespace, otherwise fall back to
+    `ccdproc._nanfuncs.median`. As of the 2023 version of the array API
+    there is no median function in the API.
+
+    Parameters
+    ----------
+    array : array_like
+        Array from which to calculate the median.
+
+    axis : int
+        Axis along which to calculate the median.
+
+    xp : array namespace, optional
+        Array namespace to use for calculations. If not provided, the
+        namespace will be determined from the array.
+
+    Returns
+    -------
+    median : array
+        Median of ``array`` along ``axis``.
+    """
+    xp = xp or array_api_compat.array_namespace(array)
+    try:
+        return xp.median(array, axis=axis)
+    except AttributeError:
+        # median is not part of the array API standard; fall back to an
+        # implementation built on nanmedian, which also matches
+        # numpy.median's NaN-propagating semantics.
+        return _nanfuncs_median(array, axis=axis, xp=xp)
+
+
 @log_to_metadata
 def ccd_process(
     ccd,
@@ -626,7 +660,7 @@ def subtract_overscan(
         overscan_axis = 0 if overscan.shape[1] > overscan.shape[0] else 1
 
     if median:
-        oscan = xp.median(overscan.data, axis=overscan_axis)
+        oscan = _median_fallback(overscan.data, overscan_axis, xp=xp)
     else:
         oscan = xp.mean(overscan.data, axis=overscan_axis)
 
diff --git a/ccdproc/tests/test_ccdproc.py b/ccdproc/tests/test_ccdproc.py
index 1d27b1cc..683ca05e 100644
--- a/ccdproc/tests/test_ccdproc.py
+++ b/ccdproc/tests/test_ccdproc.py
@@ -1,5 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 
+import types
 import warnings
 
 import array_api_compat
@@ -26,6 +27,7 @@
 from ccdproc.conftest import testing_array_library as xp
 from ccdproc.core import (
     Keyword,
+    _median_fallback,
     ccd_process,
     cosmicray_lacosmic,
     cosmicray_median,
@@ -365,6 +367,26 @@ def test_subtract_overscan_fails():
         subtract_overscan(xp.zeros((10, 10)), fits_section="[1:10]")
 
 
+def test_median_fallback_without_native_median():
+    # The except branch of _median_fallback only runs naturally on namespaces
+    # with no ``median`` (array-api-strict is the only such backend in CI, and
+    # it does not report coverage), so hide the native function behind a proxy
+    # namespace that otherwise delegates to the backend under test.
+    class _NamespaceWithoutMedian(types.ModuleType):
+        def __getattr__(self, name):
+            if name == "median":
+                raise AttributeError(name)
+            return getattr(xp, name)
+
+    proxy = _NamespaceWithoutMedian("xp_without_median")
+    data = xp.asarray(np_array([[1.0, 2.0, 4.0], [3.0, 5.0, 6.0]]), device=xp_device)
+
+    result = _median_fallback(data, 0, xp=proxy)
+
+    expected = xp.asarray(np_array([2.0, 3.5, 5.0]), device=xp_device)
+    assert xp.all(xpx.isclose(result, expected))
+
+
 def test_trim_image_fits_section_requires_string():
     ccd_data = ccd_data_func()
     with pytest.raises(TypeError):
diff --git a/ccdproc/tests/test_nanfuncs.py b/ccdproc/tests/test_nanfuncs.py
index c7237930..cd5f58c3 100644
--- a/ccdproc/tests/test_nanfuncs.py
+++ b/ccdproc/tests/test_nanfuncs.py
@@ -5,7 +5,7 @@
 import numpy as np
 import pytest
 
-from ccdproc._nanfuncs import nanmean, nanmedian, nanstd, nansum
+from ccdproc._nanfuncs import median, nanmean, nanmedian, nanstd, nansum
 from ccdproc.conftest import testing_array_device as xp_device
 from ccdproc.conftest import testing_array_library as xp
 
@@ -26,6 +26,7 @@
     # every multi-element float row below.
     pytest.param(nanstd, np.nanstd, id="nanstd"),
     pytest.param(nanmedian, np.nanmedian, id="nanmedian"),
+    pytest.param(median, np.median, id="median"),
 ]
 
 _DATA = [
@@ -74,14 +75,17 @@ def test_matches_numpy(func, reference, data, axis):
     assert xp.all(xpx.isclose(result, expected, equal_nan=True))
 
 
-@pytest.mark.parametrize("func", [nansum, nanmean, nanstd, nanmedian])
+@pytest.mark.parametrize("func", [nansum, nanmean, nanstd, nanmedian, median])
 def test_no_warning_on_all_nan_slice(func):
     """
     All-NaN slices are handled silently.
 
     Most of the numpy counterparts warn here, and ccdproc's pytest
     configuration turns warnings into errors, so a fallback that warned would
-    fail every ``Combiner`` test with a fully masked pixel.
+    fail every ``Combiner`` test with a fully masked pixel. ``numpy.median``
+    itself does not warn on NaN input, but ``median`` is included here too
+    since it shares the ``_setup``/``nanmedian`` machinery with the other
+    fallbacks.
     """
     data = xp.asarray(
         np.array([[1.0, np.nan], [2.0, np.nan], [3.0, np.nan]]), device=xp_device
@@ -101,7 +105,7 @@ def test_nansum_all_nan_slice_is_zero():
     assert xp.all(xpx.isclose(result, xp.asarray([3.0, 0.0], device=xp_device)))
 
 
-@pytest.mark.parametrize("func", [nansum, nanmean, nanstd, nanmedian])
+@pytest.mark.parametrize("func", [nansum, nanmean, nanstd, nanmedian, median])
 @pytest.mark.parametrize(
     ("axis", "error"),
     [