From b07879161617986b63876d7ef1b02d7b1711f480 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 27 Sep 2026 14:03:59 -0500
Subject: [PATCH 01/13] Return host results in the floating dtype of the input

reproject and astropy.modeling always return float64 and astroscrappy
always returns float32, and _from_numpy passed that dtype through, so
wcs_project turned float32 data into float64 and cosmicray_lacosmic
turned float64 data into float32. _from_numpy now casts a floating
result to the dtype of like when like is floating too. Boolean masks
and results for integer input keep their own dtype, so integer data
still gets a floating result.

This flips test_from_numpy_restores_namespace_device_and_dtype, which
pinned the old pass-through, and adds a wcs_project dtype test. Part
of #1023.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RgzdGvWt8d9CYfT2ayUpxz
---
 ccdproc/core.py                | 33 ++++++++++--
 ccdproc/tests/test_hostcopy.py | 94 ++++++++++++++++++++++++++++++----
 2 files changed, 113 insertions(+), 14 deletions(-)

diff --git a/ccdproc/core.py b/ccdproc/core.py
index d094cde7..03abae4d 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -215,9 +215,11 @@ def _from_numpy(arr, like=None, *, xp=None):
     array or None
         ``arr`` as an array of ``xp`` on the device of ``like`` (on the
         default device of ``xp`` if ``like`` is `None`), or `None` if
-        ``arr`` is `None`. The result has the dtype of ``arr``, not of
-        ``like``, unless ``xp`` cannot represent it: JAX without 64-bit
-        mode silently converts float64 to float32.
+        ``arr`` is `None`. If ``arr`` and ``like`` both have a real
+        floating dtype, the result has the dtype of ``like``; otherwise,
+        including when there is no ``like``, it has the dtype of ``arr``,
+        unless ``xp`` cannot represent it: JAX without 64-bit mode
+        silently converts float64 to float32.
 
     Notes
     -----
@@ -225,12 +227,28 @@ def _from_numpy(arr, like=None, *, xp=None):
     NumPy-only operation to the caller's array namespace and device, so
     that the caller never receives a NumPy array in place of what it
     passed in.
+
+    The NumPy-only libraries decide the dtype of their own output:
+    ``reproject`` and ``astropy.modeling`` always return float64 and
+    astroscrappy always returns float32. Floating input keeps its dtype
+    instead (#1023), so a floating result is cast to the floating dtype of
+    ``like``. A boolean result, such as a mask, and the result for integer
+    input keep their own dtype: casting a floating result back to an
+    integer dtype would wrap negative values and could not hold the NaN
+    that ``wcs_project`` needs.
     """
     if arr is None:
         return None
     xp = xp or array_api_compat.array_namespace(like)
-    device = None if like is None else array_api_compat.device(like)
-    return xp.asarray(arr, device=device)
+    device = None
+    dtype = None
+    if like is not None:
+        device = array_api_compat.device(like)
+        if np.isdtype(arr.dtype, "real floating") and xp.isdtype(
+            like.dtype, "real floating"
+        ):
+            dtype = like.dtype
+    return xp.asarray(arr, dtype=dtype, device=device)
 
 
 def _ccddata_from_numpy(ccd, like=None, *, xp=None):
@@ -1873,6 +1891,11 @@ def wcs_project(ccd, target_wcs, target_shape=None, order="bilinear"):
     data and mask are copied to the host, and the reprojected data and mask
     are copied back to the array namespace and device of the input; the
     copy is announced with a `HostCopyWarning`.
+
+    Floating data keeps its dtype, even though ``reproject`` works in
+    float64. Integer data gives reproject's float64 result, which is not
+    cast back to an integer dtype because pixels that fall outside the
+    input image are NaN.
     """
     from astropy.nddata.ccddata import _generate_wcs_and_update_header
     from reproject import reproject_interp
diff --git a/ccdproc/tests/test_hostcopy.py b/ccdproc/tests/test_hostcopy.py
index 909f1e73..86b0b068 100644
--- a/ccdproc/tests/test_hostcopy.py
+++ b/ccdproc/tests/test_hostcopy.py
@@ -138,17 +138,48 @@ def test_host_copy_warning_is_public_and_an_astropy_warning():
     assert issubclass(HostCopyWarning, AstropyUserWarning)
 
 
-@pytest.mark.parametrize("dtype", ["float64", "bool"])
-def test_from_numpy_restores_namespace_device_and_dtype(dtype):
+# Each case is a result one of the NumPy-only libraries hands back, for an
+# input of the dtype of ``like``: reproject and astropy.modeling return
+# float64, astroscrappy returns float32, and a mask comes back as bool.
+@pytest.mark.parametrize(
+    ("result_dtype", "like_dtype", "expected_dtype"),
+    [
+        # Floating input keeps its dtype (#1023), in both directions.
+        ("float64", "float32", "float32"),
+        ("float32", "float64", "float64"),
+        ("float64", "float64", "float64"),
+        # A mask stays a mask.
+        ("bool", "float32", "bool"),
+        # Integer input keeps the library's floating result.
+        ("float64", "uint16", "float64"),
+        ("float32", "int32", "float32"),
+    ],
+)
+def test_from_numpy_restores_namespace_device_and_dtype(
+    result_dtype, like_dtype, expected_dtype
+):
     """
     ``_from_numpy`` is the inverse of ``_to_numpy``: it puts a NumPy result
     back in the namespace *and* on the device of the array it came from,
-    keeping its dtype. The device half only bites on a backend with more
-    than one device (array-api-strict on its non-default device), which is
-    exactly the case that motivated the helper.
-    """
-    like = xp.asarray(np.zeros((3, 3)), device=xp_device)
-    np_array = np.ones((3, 3), dtype=dtype)
+    and a floating result back in the floating dtype of that array.
+
+    Notes
+    -----
+    The device half only bites on a backend with more than one device
+    (array-api-strict on its non-default device), which is exactly the case
+    that motivated the helper.
+
+    The dtype half is #1023: without it, the dtype of the output of
+    ``wcs_project``, ``subtract_overscan(model=...)`` and
+    ``cosmicray_lacosmic`` was set by the NumPy-only library doing the
+    work, so ``float32`` input came back as ``float64`` from the first two
+    and ``float64`` input as ``float32`` from the third. Only floating
+    input is matched: a floating result for integer input must not be cast
+    back to an integer dtype, where negative values would wrap and NaN
+    cannot be stored, and a boolean mask must stay boolean.
+    """
+    like = xp.asarray(np.zeros((3, 3), dtype=like_dtype), device=xp_device)
+    np_array = np.ones((3, 3), dtype=result_dtype)
 
     result = _from_numpy(np_array, like=like)
 
@@ -156,7 +187,7 @@ def test_from_numpy_restores_namespace_device_and_dtype(dtype):
         like
     )
     assert array_api_compat.device(result) == array_api_compat.device(like)
-    assert result.dtype == getattr(xp, dtype)
+    assert result.dtype == getattr(xp, expected_dtype)
 
 
 def test_from_numpy_passes_none_through():
@@ -437,3 +468,48 @@ def test_scalar_quantity_gain_keeps_namespace_and_device():
     cleaned = cosmicray_lacosmic(ccd, gain=gain, gain_apply=True)
     assert array_api_compat.array_namespace(cleaned.data) is input_namespace
     assert array_api_compat.device(cleaned.data) == array_api_compat.device(ccd.data)
+
+
+def _ccd_of_dtype(dtype):
+    """
+    A ``DATA_SIZE`` square image around 1000 ADU, in ``dtype``.
+
+    The level is that of a raw frame, so that the same values make sense
+    as floating and as unsigned integer data.
+    """
+    ccd = ccd_data_func(data_size=DATA_SIZE, data_mean=1000.0, data_scale=10.0)
+    ccd.data = xp.astype(ccd.data, getattr(xp, dtype))
+    return ccd
+
+
+@pytest.mark.parametrize(
+    ("dtype", "expected_dtype"),
+    [
+        ("float32", "float32"),
+        ("float64", "float64"),
+        ("uint16", "float64"),
+        ("int32", "float64"),
+    ],
+)
+def test_wcs_project_keeps_floating_dtype(dtype, expected_dtype):
+    """
+    ``wcs_project`` returns floating data in its own dtype, and integer data
+    as reproject's float64 (#1023).
+
+    Notes
+    -----
+    ``reproject`` always works in float64, and its result used to be
+    returned as it was, so float32 data came back as float64. Integer data
+    is not cast back to its own dtype: pixels that fall outside the input
+    image are NaN, which an integer dtype cannot hold. The shifted target
+    WCS makes sure there are such pixels.
+    """
+    ccd = _ccd_of_dtype(dtype)
+    ccd.wcs = wcs_for_testing(ccd.shape)
+    target_wcs = wcs_for_testing(ccd.shape)
+    target_wcs.wcs.crpix += [1, 1]
+
+    result = wcs_project(ccd, target_wcs)
+
+    assert result.data.dtype == getattr(xp, expected_dtype)
+    assert bool(xp.any(xp.isnan(result.data)))

From 69f7edaae7a98213905111f73d509ef966a88131 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 27 Sep 2026 14:51:26 -0500
Subject: [PATCH 02/13] Give subtract_overscan's fitted overscan the plain
 path's dtype

The fitted overscan is now converted like the reduced overscan instead
of like ccd.data. The reduced overscan is always floating, so a model
gives the same dtype as the plain mean or median path for every input,
including integer data, where that dtype depends on the array library
(for uint16 and int32, float64 on NumPy but float32 on JAX). Converting
like ccd.data would have kept astropy.modeling's float64 for integer
data on JAX.

The integer cases of the new test are a non-strict backend_xfail on
array-api-strict: subtract_overscan cannot reduce an integer overscan
there because of ccdproc's own arithmetic, not a limit of the backend,
so they will XPASS once that is fixed rather than be skipped where
nobody would notice. The reason points at the array API tracking issue
#971. Part of #1023.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RgzdGvWt8d9CYfT2ayUpxz
---
 ccdproc/core.py                | 12 +++++++++-
 ccdproc/tests/test_hostcopy.py | 44 ++++++++++++++++++++++++++++++++++
 2 files changed, 55 insertions(+), 1 deletion(-)

diff --git a/ccdproc/core.py b/ccdproc/core.py
index 03abae4d..9f812881 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -1250,6 +1250,11 @@ def subtract_overscan(
     it is subtracted. The median and mean paths stay in the input's array
     namespace.
 
+    The result has the same dtype with or without a ``model``, even though
+    ``astropy.modeling`` fits in float64: the fitted overscan is converted
+    to the dtype of the mean or median of the overscan. Floating data
+    therefore keeps its dtype.
+
     Examples
     --------
     Creating a 100x100 array containing ones just for demonstration purposes::
@@ -1309,7 +1314,12 @@ def subtract_overscan(
         of = fitting.LinearLSQFitter()
         yarr = np.arange(oscan_np.shape[0])
         fitted = of(model, yarr, oscan_np)
-        oscan = _from_numpy(fitted(yarr), like=ccd.data, xp=xp)
+        # Convert like the reduced overscan rather than like ccd.data: it is
+        # always floating, so the fit comes back in the dtype the median and
+        # mean paths subtract, including for integer data, where that dtype
+        # depends on the array library (for example float64 for uint16 on
+        # NumPy, float32 on JAX).
+        oscan = _from_numpy(fitted(yarr), like=oscan, xp=xp)
         if overscan_axis == 1:
             oscan = xp.reshape(oscan, (oscan.size, 1))
         else:
diff --git a/ccdproc/tests/test_hostcopy.py b/ccdproc/tests/test_hostcopy.py
index 86b0b068..b495ea85 100644
--- a/ccdproc/tests/test_hostcopy.py
+++ b/ccdproc/tests/test_hostcopy.py
@@ -482,6 +482,50 @@ def _ccd_of_dtype(dtype):
     return ccd
 
 
+_NO_INTEGER_OVERSCAN_ON_STRICT = pytest.mark.backend_xfail(
+    "array-api-strict",
+    reason="array-api-strict allows the mean only of floating data, so "
+    "subtract_overscan cannot reduce an integer overscan there yet "
+    "(https://github.com/astropy/ccdproc/issues/971)",
+)
+
+
+@pytest.mark.parametrize(
+    "dtype",
+    [
+        "float32",
+        "float64",
+        pytest.param("uint16", marks=_NO_INTEGER_OVERSCAN_ON_STRICT),
+        pytest.param("int32", marks=_NO_INTEGER_OVERSCAN_ON_STRICT),
+    ],
+)
+def test_subtract_overscan_model_keeps_dtype_of_plain_path(dtype):
+    """
+    Fitting a ``model`` to the overscan gives the same dtype as subtracting
+    its plain mean, so floating data keeps its dtype (#1023).
+
+    Notes
+    -----
+    ``astropy.modeling`` always fits in float64, and the fitted overscan
+    used to come back as float64, so float32 data became float64 only when
+    a ``model`` was given. Integer data gets whatever floating dtype the
+    mean of the overscan has, which for uint16 and int32 is float64 on NumPy
+    but float32 on JAX, so the two paths are compared with each other rather
+    than with a fixed dtype.
+    """
+    ccd = _ccd_of_dtype(dtype)
+    overscan = ccd[:, :5]
+
+    plain = subtract_overscan(ccd, overscan=overscan, overscan_axis=1)
+    fitted = subtract_overscan(
+        ccd, overscan=overscan, overscan_axis=1, model=models.Polynomial1D(1)
+    )
+
+    assert fitted.data.dtype == plain.data.dtype
+    if xp.isdtype(ccd.data.dtype, "real floating"):
+        assert fitted.data.dtype == ccd.data.dtype
+
+
 @pytest.mark.parametrize(
     ("dtype", "expected_dtype"),
     [

From 4dd276802371fd61caeaff347eb23b9a93af6be3 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 27 Sep 2026 14:51:28 -0500
Subject: [PATCH 03/13] Keep the pixels astroscrappy does not flag exactly

astroscrappy scales every pixel by the gain, rounds it to float32 and
scales it back, so its output differs from the input even where it
finds no cosmic ray: by up to 6e-5 at 1000 ADU for float64 data, and
for float32 data when the gain is not a power of two. _lacosmic_on_host
now takes only the flagged pixels from astroscrappy and keeps the rest
of the input as it was.

The merge, the removal of the pssl offset and the choice of dtype all
happen on the host copy that astroscrappy needs anyway, and the result
is converted back once. Taking the unflagged pixels from the caller's
array instead would put a lazy (dask) input's whole graph into the
result a second time, so computing the cleaned image would re-run every
earlier step of the pipeline.

The dtype of the result is named rather than left to type promotion,
which differs between integer types (NumPy gives float32 for uint16 but
float64 for int32): floating input keeps its dtype and integer input
gets astroscrappy's float32, as decided in #1023.

Tests pin both dtypes, that the unflagged pixels come back bit for bit
with and without a pssl, and, on dask only, that the input's blocks are
evaluated once. The docstring Notes say the unflagged pixels are exact
apart from the gain when gain_apply is set, and why integer input gives
float32. Part of #1023.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RgzdGvWt8d9CYfT2ayUpxz
---
 ccdproc/core.py                 |  39 +++++++--
 ccdproc/tests/test_cosmicray.py | 141 ++++++++++++++++++++++++++++++++
 2 files changed, 173 insertions(+), 7 deletions(-)

diff --git a/ccdproc/core.py b/ccdproc/core.py
index 9f812881..471e8313 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -2726,6 +2726,13 @@ def cosmicray_lacosmic(
     `HostCopyWarning`. The cleaned data and cosmic-ray mask are returned in
     the array namespace and on the device of the input.
 
+    astroscrappy works in float32 and rounds every pixel to float32, so
+    only the pixels it flags as cosmic rays are taken from its output; the
+    rest keep their input values exactly, apart from the multiplication by
+    the gain when ``gain_apply`` is set. Floating data keeps its dtype,
+    with float32 precision in the replaced pixels only. Integer data gives
+    a float32 result, the dtype astroscrappy provides, whatever the gain.
+
     Returns
     -------
     nccd : `~astropy.nddata.CCDData` or array
@@ -2970,12 +2977,16 @@ def _lacosmic_on_host(
     cleanarr, crmask : array
         The cleaned data, not yet gain-corrected, and the boolean
         cosmic-ray mask, in the namespace of ``data`` and on its device.
+        ``cleanarr`` has the dtype of ``data`` if that is floating, and is
+        float32 otherwise; only the pixels in ``crmask`` differ from
+        ``data`` converted to that dtype.
 
     Notes
     -----
     This is the part of `cosmicray_lacosmic` shared by its ``CCDData`` and
     bare-array input: the `HostCopyWarning`, the copy of every input array
-    to the host, the astroscrappy call, and the copy of the results back.
+    to the host, the astroscrappy call, the merge of its cleaned pixels with
+    the input, and the copy of the results back.
     Validation that should fail before any copy stays in the caller, and
     so does everything that depends on the astroscrappy version: the old
     interface's ``pssl`` goes in ``detect_kwargs``, and its ``inbkg`` and
@@ -3005,21 +3016,35 @@ def _lacosmic_on_host(
     # to add a zero offset just to copy the data: astroscrappy only reads
     # its input and works on a copy of its own.
     host_data = _to_numpy(data)
-    if data_offset != 0:
-        host_data = host_data + data_offset
     crmask, cleanarr = detect_cosmics(
-        host_data,
+        host_data + data_offset if data_offset != 0 else host_data,
         inmask=None if mask is None else _to_numpy(mask),
         gain=gain,
         **detect_kwargs,
         **background_kwargs,
     )
 
-    # Back to the caller's namespace and device before any arithmetic, so
-    # that everything below runs natively.
-    cleanarr = _from_numpy(cleanarr, like=data, xp=xp)
+    # astroscrappy rounds every pixel to float32, and scales by the gain and
+    # back, so take only the pixels it replaced from its output and keep the
+    # rest of the input exactly. The dtype is named rather than left to type
+    # promotion, which for integer input depends on the integer type (NumPy
+    # gives float32 for uint16 but float64 for int32): floating input keeps
+    # its own dtype, and integer input gets astroscrappy's float32. This is
+    # done on the host copy of the input, which already exists: using
+    # ``data`` itself would put its graph into the result a second time for
+    # a lazy namespace such as dask, so computing the result would compute
+    # everything upstream of this call again.
+    if np.isdtype(host_data.dtype, "real floating"):
+        dtype = host_data.dtype
+    else:
+        dtype = np.float32
+    cleanarr = cleanarr.astype(dtype, copy=False)
     if data_offset != 0:
         cleanarr = cleanarr - data_offset
+    cleanarr = np.where(crmask, cleanarr, host_data.astype(dtype, copy=False))
+
+    # Back to the caller's namespace and device, once, at the end.
+    cleanarr = _from_numpy(cleanarr, like=data, xp=xp)
     crmask = _from_numpy(crmask, like=data, xp=xp)
     return cleanarr, crmask
 
diff --git a/ccdproc/tests/test_cosmicray.py b/ccdproc/tests/test_cosmicray.py
index 1cbc5a59..38d5deca 100644
--- a/ccdproc/tests/test_cosmicray.py
+++ b/ccdproc/tests/test_cosmicray.py
@@ -929,6 +929,147 @@ def no_cosmics(data, **_kwargs):
     assert_allclose(_to_numpy(cleaned), frame, rtol=1e-6)
 
 
+def _lacosmic_frame(dtype):
+    """
+    A raw-looking frame around 1000 ADU in ``dtype``, with cosmic rays.
+
+    Returns the frame as a ``CCDData`` on the test backend.
+    """
+    ccd_data = ccd_data_func(data_scale=DATA_SCALE, data_mean=1000.0)
+    add_cosmicrays(ccd_data, DATA_SCALE, threshold=10, ncrays=NCRAYS)
+    ccd_data.data = xp.astype(ccd_data.data, getattr(xp, dtype))
+    return ccd_data
+
+
+@pytest.mark.parametrize("array_input", [True, False])
+@pytest.mark.parametrize("dtype", ["float32", "float64"])
+def test_cosmicray_lacosmic_keeps_floating_dtype(dtype, array_input):
+    """
+    Floating data comes back from ``cosmicray_lacosmic`` in its own dtype,
+    for both a ``CCDData`` and a bare array (#1023).
+
+    Notes
+    -----
+    astroscrappy always returns float32, and its result used to be returned
+    as it was, so float64 data came back as float32.
+    """
+    ccd_data = _lacosmic_frame(dtype)
+    ccd = ccd_data.data if array_input else ccd_data
+
+    result = cosmicray_lacosmic(ccd, sigclip=5.9, gain=2.0)
+
+    cleaned = result[0] if array_input else result.data
+    assert cleaned.dtype == getattr(xp, dtype)
+
+
+@pytest.mark.parametrize("gain", [1.0, 2.0])
+@pytest.mark.parametrize("dtype", ["uint16", "int32"])
+def test_cosmicray_lacosmic_integer_input_gives_float32(dtype, gain):
+    """
+    Integer data comes back from ``cosmicray_lacosmic`` as float32, whatever
+    the integer type and whatever the gain.
+
+    Notes
+    -----
+    float32 is the only dtype astroscrappy provides. This was settled in
+    #1023: floating input keeps its dtype, but integer input is not
+    promoted to float64 just to match other ccdproc functions. The integer
+    types are the ones NumPy promotes differently with float32 (uint16
+    gives float32, int32 gives float64), so a result left to type promotion
+    would differ between them. In 2.5.1 a gain other than 1 promoted the
+    result to float64 by accident, through a NumPy float64 scalar; the
+    ``gain`` parametrization pins that it no longer does.
+    """
+    ccd_data = _lacosmic_frame(dtype)
+
+    result = cosmicray_lacosmic(ccd_data, sigclip=5.9, gain=gain)
+
+    assert result.data.dtype == xp.float32
+
+
+@pytest.mark.parametrize("pssl", [0.0, 100.0])
+@pytest.mark.parametrize("dtype", ["float32", "float64"])
+def test_cosmicray_lacosmic_keeps_unflagged_pixels_exactly(dtype, pssl):
+    """
+    Only the pixels flagged as cosmic rays differ from the input; every
+    other pixel comes back bit for bit (#1023).
+
+    Notes
+    -----
+    astroscrappy multiplies the data by the gain, rounds every pixel to
+    float32, and divides by the gain again at the end, so taking its whole
+    output changed every pixel of float64 data, by up to 6e-5 at 1000 ADU,
+    and even pixels of float32 data when the gain is not a power of two, as
+    here. ``gain_apply=False`` keeps ccdproc's own gain multiplication out
+    of the comparison.
+
+    A ``pssl`` (a sky level of 100 ADU here) is added to the data before
+    astroscrappy runs and removed from its output afterwards. The unflagged
+    pixels must come straight from the input, not from the offset data with
+    the offset taken off again, which is not exact in floating point; the
+    ``pssl`` case pins that. ``pssl`` is deprecated, so passing it warns
+    even when it is zero, which is the case without an offset.
+    """
+    ccd_data = _lacosmic_frame(dtype)
+
+    with pytest.warns(AstropyDeprecationWarning):
+        cleaned, crmask = cosmicray_lacosmic(
+            ccd_data.data, sigclip=5.9, gain=1.7, gain_apply=False, pssl=pssl
+        )
+
+    np_cleaned = _to_numpy(cleaned)
+    np_data = _to_numpy(ccd_data.data)
+    np_crmask = _to_numpy(crmask)
+    # Every planted cosmic ray is found (with some of its neighbours at
+    # this gain), and the flagged pixels really were replaced.
+    assert np_crmask.sum() >= NCRAYS
+    assert not np.array_equal(np_cleaned[np_crmask], np_data[np_crmask])
+    np.testing.assert_array_equal(np_cleaned[~np_crmask], np_data[~np_crmask])
+
+
+@pytest.mark.backend_skip(
+    "numpy",
+    "jax",
+    "array-api-strict",
+    "cupy",
+    reason="only dask builds its arrays lazily",
+)
+def test_cosmicray_lacosmic_computes_lazy_input_once():
+    """
+    A lazy (dask) input is computed once, for the host copy astroscrappy
+    needs, and computing the result does not compute it again.
+
+    Notes
+    -----
+    Keeping the unflagged pixels exactly (#1023) needs the input as well as
+    astroscrappy's output. Taking them from the caller's array rather than
+    from the host copy put the whole graph upstream of the call into the
+    result a second time, so computing the cleaned image re-ran every
+    earlier step of the pipeline, such as reading the files or subtracting
+    the bias, which ``earlier_step`` stands in for here.
+    """
+    evaluated = []
+
+    def earlier_step(block):
+        evaluated.append(block.shape)
+        return block
+
+    ccd_data = _lacosmic_frame("float64")
+    # Several chunks, as a real pipeline would have. Giving the meta stops
+    # dask from calling earlier_step on small dummy arrays to work out what
+    # it returns, which would count as evaluations.
+    data = ccd_data.data.rechunk(50).map_blocks(
+        earlier_step, meta=np.empty((0, 0), dtype=ccd_data.data.dtype)
+    )
+
+    cleaned, crmask = cosmicray_lacosmic(data, sigclip=5.9, gain=2.0)
+    assert len(evaluated) == data.npartitions
+
+    _to_numpy(cleaned)
+    _to_numpy(crmask)
+    assert len(evaluated) == data.npartitions
+
+
 def test_cosmicray_median_mask_shape_mismatch():
     # CCDData and MaskedArray validate the mask shape themselves, so exercise
     # the shared helper directly.

From 5e0245c9f1d907503185ddac3d8163fdd789f07d Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 27 Sep 2026 14:51:28 -0500
Subject: [PATCH 04/13] Record the host round trips' dtype rule in CHANGES and
 the docs

Adds a Bug Fixes entry for #1023 and says in the array API guide that
floating data keeps its dtype through the host round trip. The #1011
entry's sentence about cosmicray_lacosmic returning float32 data is
narrowed to integer and float32 input, since float64 input now stays
float64, and #1023 is added to its references.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RgzdGvWt8d9CYfT2ayUpxz
---
 CHANGES.rst        | 18 +++++++++++++++---
 docs/array_api.rst |  3 +++
 2 files changed, 18 insertions(+), 3 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 3b1fa16b..3fab4dcd 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -80,9 +80,10 @@ New Features
   cleaned by ``cosmicray_lacosmic`` with a NumPy ``float64`` ``pssl`` (such
   as ``np.median`` of ``float64`` data) now stays ``float32``, as it already
   did for every other type of ``pssl``. Likewise, ``cosmicray_lacosmic`` with
-  ``gain_apply=True`` now returns astroscrappy's ``float32`` data for any
-  gain, as it already did for a gain of 1; previously any other gain
-  promoted the data to ``float64``. [#930, #933, #935, #1011]
+  ``gain_apply=True`` now returns ``float32`` data for integer and
+  ``float32`` input whatever the gain, as it already did for a gain of 1;
+  previously any other gain promoted the data to ``float64``.
+  [#930, #933, #935, #1011, #1023]
 
 Other Changes and Additions
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -131,6 +132,17 @@ Bug Fixes
 - ``wcs_project`` no longer returns a NumPy-backed ``CCDData`` for non-NumPy
   input: the reprojected data and mask are converted back to the array
   namespace and device of the input image. [#930, #1011]
+- ``wcs_project``, ``subtract_overscan`` with a ``model`` and
+  ``cosmicray_lacosmic`` now return floating data in its own dtype rather
+  than in the dtype of the NumPy-only library that does the work:
+  ``float32`` data no longer becomes ``float64`` in the first two, and
+  ``float64`` data no longer becomes ``float32`` in ``cosmicray_lacosmic``.
+  ``cosmicray_lacosmic`` also returns the pixels it does not flag as cosmic
+  rays exactly as they were (multiplied by the gain when ``gain_apply`` is
+  set); previously they came back rounded to ``float32`` by astroscrappy.
+  Integer data still gives a floating result: ``float32`` from
+  ``cosmicray_lacosmic``, ``float64`` from ``wcs_project``, and from
+  ``subtract_overscan`` the same dtype with or without a ``model``. [#1023]
 - ``flat_correct`` now checks the flat's mask with ``xp.any`` instead of the
   ``mask.any()`` method, which arrays of spec-only array-API namespaces do
   not have. [#1005]
diff --git a/docs/array_api.rst b/docs/array_api.rst
index 56189a57..27e7401f 100644
--- a/docs/array_api.rst
+++ b/docs/array_api.rst
@@ -241,6 +241,9 @@ These functions copy their input to host memory, run there, and copy every
 array they return -- the data and the mask -- back to the array namespace
 and the device of the array you passed in. You never get a `numpy`_ array
 back in place of what you handed over, whatever array library you use.
+Floating data also keeps its dtype, whatever dtype the library works in
+(``float64`` for `reproject`_ and `astropy.modeling`, ``float32`` for
+`astroscrappy`_); integer data comes back in a floating dtype.
 
 This copy does not work for `CuPy`_ arrays yet: `CuPy`_ refuses the implicit
 conversion to `numpy`_ that it relies on, so these three functions, and

From bbeb7e90cdbe1f9b82ee63fa4c86c078b07a6119 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 27 Sep 2026 15:19:34 -0500
Subject: [PATCH 05/13] Keep the dtype of files that combine reads next to
 images in memory

After the rebase onto #1028, _ccddata_from_numpy converted a file's data
and uncertainty with _from_numpy, whose #1023 rule casts a floating
result to the floating dtype of like. A file is not a result computed
from like, so combine cast a float64 file to the float32 of an image in
memory, losing precision before it cast to its dtype, and a float32
file uncertainty to float64, which the template relies on keeping.

_ccddata_from_numpy now converts with xp.asarray on the device of like,
keeping each array's dtype as its docstring says. This makes
test_ccddata_from_numpy_converts_data_and_uncertainty pass again, and a
new combine test pins the float64 precision in both list orders.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RgzdGvWt8d9CYfT2ayUpxz
---
 ccdproc/core.py                | 16 ++++++++++++----
 ccdproc/tests/test_combiner.py | 26 ++++++++++++++++++++++++++
 2 files changed, 38 insertions(+), 4 deletions(-)

diff --git a/ccdproc/core.py b/ccdproc/core.py
index 471e8313..d401bfcf 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -276,16 +276,24 @@ def _ccddata_from_numpy(ccd, like=None, *, xp=None):
     -----
     Arrays read from a FITS file may be big-endian, which only NumPy
     accepts, so each is put in native byte order with `_native_numpy`
-    before `_from_numpy` converts it.
+    before it is converted.
+
+    The conversion does not go through `_from_numpy`, whose floating dtype
+    rule is for results computed from ``like`` (#1023). A file is not such
+    a result, so its arrays keep the dtype they were read with: a
+    ``float64`` file combined with a ``float32`` image in memory must not
+    lose precision before ``combine`` casts it to its ``dtype``.
 
     The mask is left in NumPy: astropy's `~astropy.nddata.CCDData` mask
     setter converts any mask to NumPy, so a caller that needs the mask in
     ``xp`` moves it there itself.
     """
-    ccd.data = _from_numpy(_native_numpy(ccd.data), like=like, xp=xp)
+    xp = xp or array_api_compat.array_namespace(like)
+    device = None if like is None else array_api_compat.device(like)
+    ccd.data = xp.asarray(_native_numpy(ccd.data), device=device)
     if ccd.uncertainty is not None:
-        ccd.uncertainty.array = _from_numpy(
-            _native_numpy(ccd.uncertainty.array), like=like, xp=xp
+        ccd.uncertainty.array = xp.asarray(
+            _native_numpy(ccd.uncertainty.array), device=device
         )
     return ccd
 
diff --git a/ccdproc/tests/test_combiner.py b/ccdproc/tests/test_combiner.py
index d3736f8c..b5213d33 100644
--- a/ccdproc/tests/test_combiner.py
+++ b/ccdproc/tests/test_combiner.py
@@ -1803,6 +1803,32 @@ def test_combine_mixed_list_with_masked_file(tmp_path, file_first):
     assert not np.any(_to_numpy(result.mask))
 
 
+@pytest.mark.parametrize("file_first", [True, False], ids=["file-first", "ccd-first"])
+def test_combine_mixed_list_keeps_float64_file_precision(tmp_path, file_first):
+    """
+    A ``float64`` file combined with a ``float32`` image in memory keeps
+    its ``float64`` values until ``combine`` casts it to its ``dtype``.
+
+    Notes
+    -----
+    Guards the interaction of #1023 with #1025. ``combine`` reads files onto
+    the device of the image in memory with ``_ccddata_from_numpy``. If that
+    applied ``_from_numpy``'s floating dtype rule, meant for results
+    computed from an input, the file would first be cast to the ``float32``
+    of the image in memory, rounding ``1 + 2**-30`` to ``1``.
+    """
+    value = 1 + 2**-30
+    path = tmp_path / "float64.fits"
+    CCDData(np.full((3, 3), value), unit=u.adu).write(path)
+    in_memory = CCDData(to_xp(np.zeros((3, 3), dtype=np.float32)), unit=u.adu)
+    img_list = [str(path), in_memory] if file_first else [in_memory, str(path)]
+
+    result = combine(img_list)
+
+    assert result.data.dtype == xp.float64
+    np.testing.assert_array_equal(_to_numpy(result.data), np.full((3, 3), value / 2))
+
+
 def test_combine_files_keeps_uncertainty_dtype(tmp_path):
     """
     ``combine`` of files with a float32 uncertainty into ``array_package``

From ff0fb2f290d6c283dfe7eed0db0eee8b4626bf8d Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 27 Sep 2026 19:51:03 -0500
Subject: [PATCH 06/13] Take cosmicray_lacosmic's whole result from
 astroscrappy again

Reverts the merge from 4dd2768 that kept the pixels astroscrappy does
not flag exactly: rounding them to float32 is far below the noise in any
real frame, and anyone running a float32 tool should expect float32
precision. _lacosmic_on_host is back to what main does, converting
astroscrappy's output with _from_numpy and then removing the pssl
offset. The _from_numpy dtype rule still gives floating input its own
dtype and integer input float32, which the dtype tests pin.

The tests for the exact unflagged pixels and for computing a lazy input
once go with it (the second existed only because the merge used the
input a second time), as do the cosmicray_lacosmic Notes paragraph, the
_lacosmic_on_host docstring lines and the CHANGES sentence that
described the behaviour. Part of #1023.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RgzdGvWt8d9CYfT2ayUpxz
---
 CHANGES.rst                     |  3 --
 ccdproc/core.py                 | 39 +++-------------
 ccdproc/tests/test_cosmicray.py | 83 ---------------------------------
 3 files changed, 7 insertions(+), 118 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index 3fab4dcd..e1b3b190 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -137,9 +137,6 @@ Bug Fixes
   than in the dtype of the NumPy-only library that does the work:
   ``float32`` data no longer becomes ``float64`` in the first two, and
   ``float64`` data no longer becomes ``float32`` in ``cosmicray_lacosmic``.
-  ``cosmicray_lacosmic`` also returns the pixels it does not flag as cosmic
-  rays exactly as they were (multiplied by the gain when ``gain_apply`` is
-  set); previously they came back rounded to ``float32`` by astroscrappy.
   Integer data still gives a floating result: ``float32`` from
   ``cosmicray_lacosmic``, ``float64`` from ``wcs_project``, and from
   ``subtract_overscan`` the same dtype with or without a ``model``. [#1023]
diff --git a/ccdproc/core.py b/ccdproc/core.py
index d401bfcf..421feab0 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -2734,13 +2734,6 @@ def cosmicray_lacosmic(
     `HostCopyWarning`. The cleaned data and cosmic-ray mask are returned in
     the array namespace and on the device of the input.
 
-    astroscrappy works in float32 and rounds every pixel to float32, so
-    only the pixels it flags as cosmic rays are taken from its output; the
-    rest keep their input values exactly, apart from the multiplication by
-    the gain when ``gain_apply`` is set. Floating data keeps its dtype,
-    with float32 precision in the replaced pixels only. Integer data gives
-    a float32 result, the dtype astroscrappy provides, whatever the gain.
-
     Returns
     -------
     nccd : `~astropy.nddata.CCDData` or array
@@ -2985,16 +2978,12 @@ def _lacosmic_on_host(
     cleanarr, crmask : array
         The cleaned data, not yet gain-corrected, and the boolean
         cosmic-ray mask, in the namespace of ``data`` and on its device.
-        ``cleanarr`` has the dtype of ``data`` if that is floating, and is
-        float32 otherwise; only the pixels in ``crmask`` differ from
-        ``data`` converted to that dtype.
 
     Notes
     -----
     This is the part of `cosmicray_lacosmic` shared by its ``CCDData`` and
     bare-array input: the `HostCopyWarning`, the copy of every input array
-    to the host, the astroscrappy call, the merge of its cleaned pixels with
-    the input, and the copy of the results back.
+    to the host, the astroscrappy call, and the copy of the results back.
     Validation that should fail before any copy stays in the caller, and
     so does everything that depends on the astroscrappy version: the old
     interface's ``pssl`` goes in ``detect_kwargs``, and its ``inbkg`` and
@@ -3024,35 +3013,21 @@ def _lacosmic_on_host(
     # to add a zero offset just to copy the data: astroscrappy only reads
     # its input and works on a copy of its own.
     host_data = _to_numpy(data)
+    if data_offset != 0:
+        host_data = host_data + data_offset
     crmask, cleanarr = detect_cosmics(
-        host_data + data_offset if data_offset != 0 else host_data,
+        host_data,
         inmask=None if mask is None else _to_numpy(mask),
         gain=gain,
         **detect_kwargs,
         **background_kwargs,
     )
 
-    # astroscrappy rounds every pixel to float32, and scales by the gain and
-    # back, so take only the pixels it replaced from its output and keep the
-    # rest of the input exactly. The dtype is named rather than left to type
-    # promotion, which for integer input depends on the integer type (NumPy
-    # gives float32 for uint16 but float64 for int32): floating input keeps
-    # its own dtype, and integer input gets astroscrappy's float32. This is
-    # done on the host copy of the input, which already exists: using
-    # ``data`` itself would put its graph into the result a second time for
-    # a lazy namespace such as dask, so computing the result would compute
-    # everything upstream of this call again.
-    if np.isdtype(host_data.dtype, "real floating"):
-        dtype = host_data.dtype
-    else:
-        dtype = np.float32
-    cleanarr = cleanarr.astype(dtype, copy=False)
+    # Back to the caller's namespace and device before any arithmetic, so
+    # that everything below runs natively.
+    cleanarr = _from_numpy(cleanarr, like=data, xp=xp)
     if data_offset != 0:
         cleanarr = cleanarr - data_offset
-    cleanarr = np.where(crmask, cleanarr, host_data.astype(dtype, copy=False))
-
-    # Back to the caller's namespace and device, once, at the end.
-    cleanarr = _from_numpy(cleanarr, like=data, xp=xp)
     crmask = _from_numpy(crmask, like=data, xp=xp)
     return cleanarr, crmask
 
diff --git a/ccdproc/tests/test_cosmicray.py b/ccdproc/tests/test_cosmicray.py
index 38d5deca..f65483b1 100644
--- a/ccdproc/tests/test_cosmicray.py
+++ b/ccdproc/tests/test_cosmicray.py
@@ -987,89 +987,6 @@ def test_cosmicray_lacosmic_integer_input_gives_float32(dtype, gain):
     assert result.data.dtype == xp.float32
 
 
-@pytest.mark.parametrize("pssl", [0.0, 100.0])
-@pytest.mark.parametrize("dtype", ["float32", "float64"])
-def test_cosmicray_lacosmic_keeps_unflagged_pixels_exactly(dtype, pssl):
-    """
-    Only the pixels flagged as cosmic rays differ from the input; every
-    other pixel comes back bit for bit (#1023).
-
-    Notes
-    -----
-    astroscrappy multiplies the data by the gain, rounds every pixel to
-    float32, and divides by the gain again at the end, so taking its whole
-    output changed every pixel of float64 data, by up to 6e-5 at 1000 ADU,
-    and even pixels of float32 data when the gain is not a power of two, as
-    here. ``gain_apply=False`` keeps ccdproc's own gain multiplication out
-    of the comparison.
-
-    A ``pssl`` (a sky level of 100 ADU here) is added to the data before
-    astroscrappy runs and removed from its output afterwards. The unflagged
-    pixels must come straight from the input, not from the offset data with
-    the offset taken off again, which is not exact in floating point; the
-    ``pssl`` case pins that. ``pssl`` is deprecated, so passing it warns
-    even when it is zero, which is the case without an offset.
-    """
-    ccd_data = _lacosmic_frame(dtype)
-
-    with pytest.warns(AstropyDeprecationWarning):
-        cleaned, crmask = cosmicray_lacosmic(
-            ccd_data.data, sigclip=5.9, gain=1.7, gain_apply=False, pssl=pssl
-        )
-
-    np_cleaned = _to_numpy(cleaned)
-    np_data = _to_numpy(ccd_data.data)
-    np_crmask = _to_numpy(crmask)
-    # Every planted cosmic ray is found (with some of its neighbours at
-    # this gain), and the flagged pixels really were replaced.
-    assert np_crmask.sum() >= NCRAYS
-    assert not np.array_equal(np_cleaned[np_crmask], np_data[np_crmask])
-    np.testing.assert_array_equal(np_cleaned[~np_crmask], np_data[~np_crmask])
-
-
-@pytest.mark.backend_skip(
-    "numpy",
-    "jax",
-    "array-api-strict",
-    "cupy",
-    reason="only dask builds its arrays lazily",
-)
-def test_cosmicray_lacosmic_computes_lazy_input_once():
-    """
-    A lazy (dask) input is computed once, for the host copy astroscrappy
-    needs, and computing the result does not compute it again.
-
-    Notes
-    -----
-    Keeping the unflagged pixels exactly (#1023) needs the input as well as
-    astroscrappy's output. Taking them from the caller's array rather than
-    from the host copy put the whole graph upstream of the call into the
-    result a second time, so computing the cleaned image re-ran every
-    earlier step of the pipeline, such as reading the files or subtracting
-    the bias, which ``earlier_step`` stands in for here.
-    """
-    evaluated = []
-
-    def earlier_step(block):
-        evaluated.append(block.shape)
-        return block
-
-    ccd_data = _lacosmic_frame("float64")
-    # Several chunks, as a real pipeline would have. Giving the meta stops
-    # dask from calling earlier_step on small dummy arrays to work out what
-    # it returns, which would count as evaluations.
-    data = ccd_data.data.rechunk(50).map_blocks(
-        earlier_step, meta=np.empty((0, 0), dtype=ccd_data.data.dtype)
-    )
-
-    cleaned, crmask = cosmicray_lacosmic(data, sigclip=5.9, gain=2.0)
-    assert len(evaluated) == data.npartitions
-
-    _to_numpy(cleaned)
-    _to_numpy(crmask)
-    assert len(evaluated) == data.npartitions
-
-
 def test_cosmicray_median_mask_shape_mismatch():
     # CCDData and MaskedArray validate the mask shape themselves, so exercise
     # the shared helper directly.

From d4fd577fad54f1ac8b57c6f8a64691bf6cc4dc97 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 27 Sep 2026 19:51:16 -0500
Subject: [PATCH 07/13] Trim the dtype notes in subtract_overscan and
 wcs_project docstrings

The paragraph on how the fitted overscan gets its dtype is not something
a caller of subtract_overscan needs, and CHANGES already records the
change. wcs_project keeps only the sentence saying floating data keeps
its dtype; why integer data is not cast back is for the code, not the
docstring. Part of #1023.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RgzdGvWt8d9CYfT2ayUpxz
---
 ccdproc/core.py | 9 +--------
 1 file changed, 1 insertion(+), 8 deletions(-)

diff --git a/ccdproc/core.py b/ccdproc/core.py
index 421feab0..8ba3ea2b 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -1258,11 +1258,6 @@ def subtract_overscan(
     it is subtracted. The median and mean paths stay in the input's array
     namespace.
 
-    The result has the same dtype with or without a ``model``, even though
-    ``astropy.modeling`` fits in float64: the fitted overscan is converted
-    to the dtype of the mean or median of the overscan. Floating data
-    therefore keeps its dtype.
-
     Examples
     --------
     Creating a 100x100 array containing ones just for demonstration purposes::
@@ -1911,9 +1906,7 @@ def wcs_project(ccd, target_wcs, target_shape=None, order="bilinear"):
     copy is announced with a `HostCopyWarning`.
 
     Floating data keeps its dtype, even though ``reproject`` works in
-    float64. Integer data gives reproject's float64 result, which is not
-    cast back to an integer dtype because pixels that fall outside the
-    input image are NaN.
+    float64.
     """
     from astropy.nddata.ccddata import _generate_wcs_and_update_header
     from reproject import reproject_interp

From c62c932cf9bf1319c765480c38d850b77ab9e21f Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 27 Sep 2026 19:51:24 -0500
Subject: [PATCH 08/13] Shorten the #1023 CHANGES entry

The entry says what changed for floating data and nothing more. The
integer behaviour is left out because it has not changed.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RgzdGvWt8d9CYfT2ayUpxz
---
 CHANGES.rst | 9 ++-------
 1 file changed, 2 insertions(+), 7 deletions(-)

diff --git a/CHANGES.rst b/CHANGES.rst
index e1b3b190..e4ee661a 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -133,13 +133,8 @@ Bug Fixes
   input: the reprojected data and mask are converted back to the array
   namespace and device of the input image. [#930, #1011]
 - ``wcs_project``, ``subtract_overscan`` with a ``model`` and
-  ``cosmicray_lacosmic`` now return floating data in its own dtype rather
-  than in the dtype of the NumPy-only library that does the work:
-  ``float32`` data no longer becomes ``float64`` in the first two, and
-  ``float64`` data no longer becomes ``float32`` in ``cosmicray_lacosmic``.
-  Integer data still gives a floating result: ``float32`` from
-  ``cosmicray_lacosmic``, ``float64`` from ``wcs_project``, and from
-  ``subtract_overscan`` the same dtype with or without a ``model``. [#1023]
+  ``cosmicray_lacosmic`` now return floating data in its input dtype rather
+  than in the dtype of the NumPy-only library that does the work. [#1023]
 - ``flat_correct`` now checks the flat's mask with ``xp.any`` instead of the
   ``mask.any()`` method, which arrays of spec-only array-API namespaces do
   not have. [#1005]

From b52bb1725ce90ade7e9d78b993a2b547a5195f9a Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 27 Sep 2026 19:52:15 -0500
Subject: [PATCH 09/13] Let the ccd_data test helper make an image of a given
 dtype

ccd_data gains a dtype keyword that casts its random values to that
dtype of the array library under test. This replaces the two
near-identical helpers the #1023 tests added, _ccd_of_dtype in
test_hostcopy.py and _lacosmic_frame in test_cosmicray.py, which built
the same image and cast it. The lacosmic dtype tests add their cosmic
rays inline.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RgzdGvWt8d9CYfT2ayUpxz
---
 ccdproc/tests/pytest_fixtures.py | 11 ++++++++++-
 ccdproc/tests/test_cosmicray.py  | 18 ++++--------------
 ccdproc/tests/test_hostcopy.py   | 20 ++++++--------------
 3 files changed, 20 insertions(+), 29 deletions(-)

diff --git a/ccdproc/tests/pytest_fixtures.py b/ccdproc/tests/pytest_fixtures.py
index 7cd8c041..aa271536 100644
--- a/ccdproc/tests/pytest_fixtures.py
+++ b/ccdproc/tests/pytest_fixtures.py
@@ -35,6 +35,7 @@ def ccd_data(
     data_scale=DEFAULT_DATA_SCALE,
     data_mean=DEFAULT_DATA_MEAN,
     rng_seed=DEFAULT_SEED,
+    dtype=None,
 ):
     """
     Return a CCDData object with units of ADU.
@@ -51,6 +52,10 @@ def ccd_data(
 
     The mean can be changed with the marker @pytest.marker.scale(m) on the
     test function, where m is the desired mean.
+
+    If ``dtype``, the name of a dtype of the array library under test such
+    as ``"float32"`` or ``"uint16"``, is given, the random values are cast
+    to it.
     """
     # Need the import here to avoid circular import issues
     from ..conftest import testing_array_device as xp_device
@@ -65,8 +70,12 @@ def ccd_data(
 
     data = rng.normal(loc=mean, size=[size, size], scale=scale)
 
+    data = xp.asarray(data, device=xp_device)
+    if dtype is not None:
+        data = xp.astype(data, getattr(xp, dtype))
+
     fake_meta = {"my_key": 42, "your_key": "not 42"}
-    ccd = CCDData(xp.asarray(data, device=xp_device), unit=u.adu)
+    ccd = CCDData(data, unit=u.adu)
     ccd.header = fake_meta
     return ccd
 
diff --git a/ccdproc/tests/test_cosmicray.py b/ccdproc/tests/test_cosmicray.py
index f65483b1..5c80094c 100644
--- a/ccdproc/tests/test_cosmicray.py
+++ b/ccdproc/tests/test_cosmicray.py
@@ -929,18 +929,6 @@ def no_cosmics(data, **_kwargs):
     assert_allclose(_to_numpy(cleaned), frame, rtol=1e-6)
 
 
-def _lacosmic_frame(dtype):
-    """
-    A raw-looking frame around 1000 ADU in ``dtype``, with cosmic rays.
-
-    Returns the frame as a ``CCDData`` on the test backend.
-    """
-    ccd_data = ccd_data_func(data_scale=DATA_SCALE, data_mean=1000.0)
-    add_cosmicrays(ccd_data, DATA_SCALE, threshold=10, ncrays=NCRAYS)
-    ccd_data.data = xp.astype(ccd_data.data, getattr(xp, dtype))
-    return ccd_data
-
-
 @pytest.mark.parametrize("array_input", [True, False])
 @pytest.mark.parametrize("dtype", ["float32", "float64"])
 def test_cosmicray_lacosmic_keeps_floating_dtype(dtype, array_input):
@@ -953,7 +941,8 @@ def test_cosmicray_lacosmic_keeps_floating_dtype(dtype, array_input):
     astroscrappy always returns float32, and its result used to be returned
     as it was, so float64 data came back as float32.
     """
-    ccd_data = _lacosmic_frame(dtype)
+    ccd_data = ccd_data_func(data_scale=DATA_SCALE, data_mean=1000.0, dtype=dtype)
+    add_cosmicrays(ccd_data, DATA_SCALE, threshold=10, ncrays=NCRAYS)
     ccd = ccd_data.data if array_input else ccd_data
 
     result = cosmicray_lacosmic(ccd, sigclip=5.9, gain=2.0)
@@ -980,7 +969,8 @@ def test_cosmicray_lacosmic_integer_input_gives_float32(dtype, gain):
     result to float64 by accident, through a NumPy float64 scalar; the
     ``gain`` parametrization pins that it no longer does.
     """
-    ccd_data = _lacosmic_frame(dtype)
+    ccd_data = ccd_data_func(data_scale=DATA_SCALE, data_mean=1000.0, dtype=dtype)
+    add_cosmicrays(ccd_data, DATA_SCALE, threshold=10, ncrays=NCRAYS)
 
     result = cosmicray_lacosmic(ccd_data, sigclip=5.9, gain=gain)
 
diff --git a/ccdproc/tests/test_hostcopy.py b/ccdproc/tests/test_hostcopy.py
index b495ea85..1871d9c2 100644
--- a/ccdproc/tests/test_hostcopy.py
+++ b/ccdproc/tests/test_hostcopy.py
@@ -470,18 +470,6 @@ def test_scalar_quantity_gain_keeps_namespace_and_device():
     assert array_api_compat.device(cleaned.data) == array_api_compat.device(ccd.data)
 
 
-def _ccd_of_dtype(dtype):
-    """
-    A ``DATA_SIZE`` square image around 1000 ADU, in ``dtype``.
-
-    The level is that of a raw frame, so that the same values make sense
-    as floating and as unsigned integer data.
-    """
-    ccd = ccd_data_func(data_size=DATA_SIZE, data_mean=1000.0, data_scale=10.0)
-    ccd.data = xp.astype(ccd.data, getattr(xp, dtype))
-    return ccd
-
-
 _NO_INTEGER_OVERSCAN_ON_STRICT = pytest.mark.backend_xfail(
     "array-api-strict",
     reason="array-api-strict allows the mean only of floating data, so "
@@ -513,7 +501,9 @@ def test_subtract_overscan_model_keeps_dtype_of_plain_path(dtype):
     but float32 on JAX, so the two paths are compared with each other rather
     than with a fixed dtype.
     """
-    ccd = _ccd_of_dtype(dtype)
+    ccd = ccd_data_func(
+        data_size=DATA_SIZE, data_mean=1000.0, data_scale=10.0, dtype=dtype
+    )
     overscan = ccd[:, :5]
 
     plain = subtract_overscan(ccd, overscan=overscan, overscan_axis=1)
@@ -548,7 +538,9 @@ def test_wcs_project_keeps_floating_dtype(dtype, expected_dtype):
     image are NaN, which an integer dtype cannot hold. The shifted target
     WCS makes sure there are such pixels.
     """
-    ccd = _ccd_of_dtype(dtype)
+    ccd = ccd_data_func(
+        data_size=DATA_SIZE, data_mean=1000.0, data_scale=10.0, dtype=dtype
+    )
     ccd.wcs = wcs_for_testing(ccd.shape)
     target_wcs = wcs_for_testing(ccd.shape)
     target_wcs.wcs.crpix += [1, 1]

From ed919430a4e3eeb3e1b086a999847ed30184df0e Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Sun, 27 Sep 2026 19:53:26 -0500
Subject: [PATCH 10/13] Keep the uncertainty's dtype when cosmicray_lacosmic
 applies the gain

With gain_apply=True the uncertainty is multiplied by the gain, which
was made a float64 array, so a float32 uncertainty came back as float64
on every backend except JAX. The gain now has the uncertainty's dtype
(the cleaned data's dtype if the uncertainty is not floating), so the
uncertainty keeps its own dtype and a float64 uncertainty keeps a
float64 gain.

A new test pins that for float32, float64 and uint16 data with a float32
uncertainty. array-api-strict still refuses uint16 data, because the
integer data itself is multiplied by a floating scalar; that case is a
non-strict xfail pointing at #1013. Part of #1023.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RgzdGvWt8d9CYfT2ayUpxz
---
 ccdproc/core.py                 |  9 +++++++-
 ccdproc/tests/test_cosmicray.py | 41 +++++++++++++++++++++++++++++++++
 2 files changed, 49 insertions(+), 1 deletion(-)

diff --git a/ccdproc/core.py b/ccdproc/core.py
index 8ba3ea2b..d792382e 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -2897,8 +2897,15 @@ def cosmicray_lacosmic(
 
         if gain_apply:
             if nccd.uncertainty is not None:
+                # Give the gain the uncertainty's dtype so that multiplying
+                # by it does not change that dtype.
+                unc_dtype = nccd.uncertainty.array.dtype
+                if not xp.isdtype(unc_dtype, "real floating"):
+                    unc_dtype = cleanarr.dtype
                 gain_value = xp.asarray(
-                    float(gain.value), device=array_api_compat.device(_ccd.data)
+                    float(gain.value),
+                    dtype=unc_dtype,
+                    device=array_api_compat.device(_ccd.data),
                 )
                 gain_corrected = _ccd.multiply(
                     gain_value, xp=xp, handle_mask=xp.logical_or
diff --git a/ccdproc/tests/test_cosmicray.py b/ccdproc/tests/test_cosmicray.py
index 5c80094c..92087447 100644
--- a/ccdproc/tests/test_cosmicray.py
+++ b/ccdproc/tests/test_cosmicray.py
@@ -977,6 +977,47 @@ def test_cosmicray_lacosmic_integer_input_gives_float32(dtype, gain):
     assert result.data.dtype == xp.float32
 
 
+@pytest.mark.parametrize(
+    "dtype",
+    [
+        "float32",
+        "float64",
+        pytest.param(
+            "uint16",
+            marks=pytest.mark.backend_xfail(
+                "array-api-strict",
+                reason="array-api-strict does not promote integer data with a "
+                "floating scalar, so the gain cannot be applied to it "
+                "(https://github.com/astropy/ccdproc/issues/1013)",
+            ),
+        ),
+    ],
+)
+def test_cosmicray_lacosmic_gain_apply_keeps_uncertainty_dtype(dtype):
+    """
+    With ``gain_apply=True`` the uncertainty keeps its own dtype, whatever
+    the dtype of the data (#1023).
+
+    Notes
+    -----
+    The uncertainty is multiplied by the gain, which used to be made a
+    float64 array, so a float32 uncertainty came back as float64 on every
+    backend except JAX. The gain now has the uncertainty's dtype.
+    """
+    ccd_data = ccd_data_func(
+        data_size=40, data_scale=DATA_SCALE, data_mean=1000.0, dtype=dtype
+    )
+    ccd_data.uncertainty = StdDevUncertainty(
+        xp.full(ccd_data.shape, DATA_SCALE, dtype=xp.float32, device=xp_device)
+    )
+
+    result = cosmicray_lacosmic(
+        ccd_data, gain=2.0 * u.electron / u.adu, gain_apply=True
+    )
+
+    assert result.uncertainty.array.dtype == xp.float32
+
+
 def test_cosmicray_median_mask_shape_mismatch():
     # CCDData and MaskedArray validate the mask shape themselves, so exercise
     # the shared helper directly.

From 185b745da0fdbf5393e48c387e187dafee4cebe8 Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Mon, 28 Sep 2026 09:42:42 -0500
Subject: [PATCH 11/13] Make the lacosmic gain without a branch for integer
 uncertainties

The gain now has xp.result_type(uncertainty dtype, float32): a float32
or float64 uncertainty keeps its dtype, and an integer one still gets a
floating gain rather than one rounded to an integer. This replaces the
fallback branch for non-floating uncertainties, which no test reached
and which codecov flagged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RgzdGvWt8d9CYfT2ayUpxz
---
 ccdproc/core.py | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/ccdproc/core.py b/ccdproc/core.py
index d792382e..767f0384 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -2899,12 +2899,9 @@ def cosmicray_lacosmic(
             if nccd.uncertainty is not None:
                 # Give the gain the uncertainty's dtype so that multiplying
                 # by it does not change that dtype.
-                unc_dtype = nccd.uncertainty.array.dtype
-                if not xp.isdtype(unc_dtype, "real floating"):
-                    unc_dtype = cleanarr.dtype
                 gain_value = xp.asarray(
                     float(gain.value),
-                    dtype=unc_dtype,
+                    dtype=xp.result_type(nccd.uncertainty.array.dtype, xp.float32),
                     device=array_api_compat.device(_ccd.data),
                 )
                 gain_corrected = _ccd.multiply(

From c252b05d80859731b2ef2c33beee492c51c0a2da Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Mon, 28 Sep 2026 10:02:54 -0500
Subject: [PATCH 12/13] Point the strict integer xfails at the new issue #1034

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RgzdGvWt8d9CYfT2ayUpxz
---
 ccdproc/tests/test_cosmicray.py | 2 +-
 ccdproc/tests/test_hostcopy.py  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/ccdproc/tests/test_cosmicray.py b/ccdproc/tests/test_cosmicray.py
index 92087447..d54b64c9 100644
--- a/ccdproc/tests/test_cosmicray.py
+++ b/ccdproc/tests/test_cosmicray.py
@@ -988,7 +988,7 @@ def test_cosmicray_lacosmic_integer_input_gives_float32(dtype, gain):
                 "array-api-strict",
                 reason="array-api-strict does not promote integer data with a "
                 "floating scalar, so the gain cannot be applied to it "
-                "(https://github.com/astropy/ccdproc/issues/1013)",
+                "(https://github.com/astropy/ccdproc/issues/1034)",
             ),
         ),
     ],
diff --git a/ccdproc/tests/test_hostcopy.py b/ccdproc/tests/test_hostcopy.py
index 1871d9c2..f8c244a9 100644
--- a/ccdproc/tests/test_hostcopy.py
+++ b/ccdproc/tests/test_hostcopy.py
@@ -474,7 +474,7 @@ def test_scalar_quantity_gain_keeps_namespace_and_device():
     "array-api-strict",
     reason="array-api-strict allows the mean only of floating data, so "
     "subtract_overscan cannot reduce an integer overscan there yet "
-    "(https://github.com/astropy/ccdproc/issues/971)",
+    "(https://github.com/astropy/ccdproc/issues/1034)",
 )
 
 

From 442314642f58ee1f311ff7099c947339c47ecc0c Mon Sep 17 00:00:00 2001
From: Matt Craig <mattwcraig@gmail.com>
Date: Mon, 28 Sep 2026 10:05:47 -0500
Subject: [PATCH 13/13] Pass a float32 invar to astroscrappy in
 cosmicray_lacosmic

astroscrappy copies inbkg into a float32 array but uses invar as it is,
and its cleaning routines accept only float32, so a float64 invar, such
as a variance image built with NumPy defaults, made every cleantype
raise "Item size of buffer (8 bytes) does not match size of
'float32_t'". The host copy of an array invar is now cast to float32
before it is handed to astroscrappy. A float64 inbkg already worked and
is left alone.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01RgzdGvWt8d9CYfT2ayUpxz
---
 CHANGES.rst                     |  3 +++
 ccdproc/core.py                 |  6 ++++++
 ccdproc/tests/test_cosmicray.py | 22 ++++++++++++++++++++++
 3 files changed, 31 insertions(+)

diff --git a/CHANGES.rst b/CHANGES.rst
index e4ee661a..31349f6d 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -135,6 +135,9 @@ Bug Fixes
 - ``wcs_project``, ``subtract_overscan`` with a ``model`` and
   ``cosmicray_lacosmic`` now return floating data in its input dtype rather
   than in the dtype of the NumPy-only library that does the work. [#1023]
+- ``cosmicray_lacosmic`` now accepts a ``float64`` ``invar``, which it
+  passes to astroscrappy as ``float32`` instead of letting astroscrappy
+  raise ``ValueError``. [#1032]
 - ``flat_correct`` now checks the flat's mask with ``xp.any`` instead of the
   ``mask.any()`` method, which arrays of spec-only array-API namespaces do
   not have. [#1005]
diff --git a/ccdproc/core.py b/ccdproc/core.py
index 767f0384..16ff9f61 100644
--- a/ccdproc/core.py
+++ b/ccdproc/core.py
@@ -3003,6 +3003,12 @@ def _lacosmic_on_host(
         for name, value in (("inbkg", inbkg), ("invar", invar))
         if value is not None
     }
+    # astroscrappy copies inbkg to float32 but not invar, and its cleaning
+    # routines only accept float32, so a float64 invar would make it raise.
+    if _is_array(invar):
+        background_kwargs["invar"] = background_kwargs["invar"].astype(
+            np.float32, copy=False
+        )
 
     # pssl is added on the host: an integer array plus a float pssl is
     # refused by some namespaces (array-api-strict), and adding it on the
diff --git a/ccdproc/tests/test_cosmicray.py b/ccdproc/tests/test_cosmicray.py
index d54b64c9..fea39a49 100644
--- a/ccdproc/tests/test_cosmicray.py
+++ b/ccdproc/tests/test_cosmicray.py
@@ -404,6 +404,28 @@ def test_cosmicray_lacosmic_invar_inbkg(new_args):
         cosmicray_lacosmic(ccd_data, sigclip=5.9, **new_args)
 
 
+def test_cosmicray_lacosmic_float64_invar():
+    """
+    A float64 ``invar`` finds the cosmic rays instead of making astroscrappy
+    raise.
+
+    Notes
+    -----
+    A variance image built with NumPy defaults is float64. astroscrappy
+    works in float32 and does not convert ``invar`` itself, so it raised a
+    ``ValueError`` about the item size of the buffer.
+    """
+    ccd_data = ccd_data_func(data_scale=DATA_SCALE)
+    add_cosmicrays(ccd_data, DATA_SCALE, threshold=10, ncrays=NCRAYS)
+    # The data are near zero, so their variance is the read noise squared.
+    readnoise = 6.5
+    invar = xp.full(ccd_data.shape, readnoise**2, dtype=xp.float64, device=xp_device)
+
+    result = cosmicray_lacosmic(ccd_data, sigclip=5.9, readnoise=readnoise, invar=invar)
+
+    assert count_true(result.mask) == NCRAYS
+
+
 def test_cosmicray_median_check_data():
     with pytest.raises(TypeError):
         ndata, crarr = cosmicray_median(10, thresh=5, mbox=11, error_image=DATA_SCALE)