Merge pull request #716 from mwcraig/fix-lacosmic-gain

Update interface to cosmic ray functions

CHANGES.rst

@@ -13,6 +13,8 @@ Other Changes and Additions
 Bug Fixes
 ^^^^^^^^^
 
+- Update units if gain is applied in ``cosmicray_lacosmic``. [#716, #705]
+
 2.0.1 (2019-09-05)
 ------------------
 

ccdproc/core.py

@@ -5,6 +5,7 @@
 import math
 import numbers
 import logging
+import warnings
 
 import numpy as np
 from scipy import ndimage
@@ -1306,19 +1307,27 @@ def cosmicray_lacosmic(ccd, sigclip=4.5, sigfrac=0.3,
                        satlevel=65535.0, pssl=0.0, niter=4,
                        sepmed=True, cleantype='meanmask', fsmode='median',
                        psfmodel='gauss', psffwhm=2.5, psfsize=7,
-                       psfk=None, psfbeta=4.765, verbose=False):
+                       psfk=None, psfbeta=4.765, verbose=False,
+                       gain_apply=True):
     r"""
     Identify cosmic rays through the L.A. Cosmic technique. The L.A. Cosmic
     technique identifies cosmic rays by identifying pixels based on a variation
     of the Laplacian edge detection. The algorithm is an implementation of the
     code describe in van Dokkum (2001) [1]_ as implemented by McCully (2014)
     [2]_. If you use this algorithm, please cite these two works.
 
+
+
     Parameters
     ----------
     ccd : `~astropy.nddata.CCDData` or `numpy.ndarray`
         Data to have cosmic ray cleaned.
 
+    gain_apply : bool, optional
+        If ``True``, **return gain-corrected data**, with correct units,
+        otherwise do not gain-correct the data. Default is ``True`` to
+        preserve backwards compatibility.
+
     sigclip : float, optional
         Laplacian-to-noise limit for cosmic ray detection. Lower values will
         flag more pixels as cosmic rays. Default: 4.5.
@@ -1338,7 +1347,7 @@ def cosmicray_lacosmic(ccd, sigclip=4.5, sigfrac=0.3,
         electrons for cosmic ray detection, so we need to know the sky level
         that has been subtracted so we can add it back in. Default: 0.0.
 
-    gain : float, optional
+    gain : float or `~astropy.units.Quantity`, optional
         Gain of the image (electrons / ADU). We always need to work in
         electrons for cosmic ray detection. Default: 1.0
 
@@ -1422,7 +1431,9 @@ def cosmicray_lacosmic(ccd, sigclip=4.5, sigfrac=0.3,
     nccd : `~astropy.nddata.CCDData` or `numpy.ndarray`
         An object of the same type as ccd is returned. If it is a
         `~astropy.nddata.CCDData`, the mask attribute will also be updated with
-        areas identified with cosmic rays masked.
+        areas identified with cosmic rays masked. **By default, the image is
+        multiplied by the gain.** You can control this behavior with the
+        ``gain_apply`` argument.
 
     crmask : `numpy.ndarray`
         If an `numpy.ndarray` is provided as ccd, a boolean ndarray with the
@@ -1460,32 +1471,73 @@ def cosmicray_lacosmic(ccd, sigclip=4.5, sigfrac=0.3,
        updated with the detected cosmic rays.
     """
     from astroscrappy import detect_cosmics
+
+    # If we didn't get a quantity, put them in, with unit specified by the
+    # documentation above.
+    if not isinstance(gain, u.Quantity):
+        # Gain will change the value, so use the proper units
+        gain = gain * u.electron / u.adu
+
+    # Set the units of readnoise to electrons, as specified in the
+    # documentation, if no unit is present.
+    if not isinstance(readnoise, u.Quantity):
+        readnoise = readnoise * u.electron
+
     if isinstance(ccd, np.ndarray):
         data = ccd
 
         crmask, cleanarr = detect_cosmics(
             data, inmask=None, sigclip=sigclip,
-            sigfrac=sigfrac, objlim=objlim, gain=gain,
-            readnoise=readnoise, satlevel=satlevel, pssl=pssl,
+            sigfrac=sigfrac, objlim=objlim, gain=gain.value,
+            readnoise=readnoise.value, satlevel=satlevel, pssl=pssl,
             niter=niter, sepmed=sepmed, cleantype=cleantype,
             fsmode=fsmode, psfmodel=psfmodel, psffwhm=psffwhm,
             psfsize=psfsize, psfk=psfk, psfbeta=psfbeta,
             verbose=verbose)
 
+        if not gain_apply and gain != 1.0:
+            cleanarr = cleanarr / gain
         return cleanarr, crmask
 
     elif isinstance(ccd, CCDData):
+        # Start with a check for a special case: ccd is in electron, and
+        # gain and readnoise have no units. In that case we issue a warning
+        # instead of raising an error to avoid crashing user's pipelines.
+        if ccd.unit.is_equivalent(u.electron) and gain.value != 1.0:
+            warnings.warn("Image unit is electron but gain value "
+                          "is not 1.0. Data maybe end up being gain "
+                          "corrected twice.")
+
+        else:
+            if ((readnoise.unit == u.electron)
+                and (ccd.unit == u.electron)
+                and (gain.value == 1.0)):
+                gain = gain.value * u.one
+            # Check unit consistency before taking the time to check for
+            # cosmic rays.
+            if not (gain * ccd).unit.is_equivalent(readnoise.unit):
+                raise ValueError('Inconsistent units for gain ({}) '.format(gain.unit) +
+                                 ' ccd ({}) and readnoise ({}).'.format(ccd.unit,
+                                                                        readnoise.unit))
 
         crmask, cleanarr = detect_cosmics(
             ccd.data, inmask=ccd.mask,
-            sigclip=sigclip, sigfrac=sigfrac, objlim=objlim, gain=gain,
-            readnoise=readnoise, satlevel=satlevel, pssl=pssl,
+            sigclip=sigclip, sigfrac=sigfrac, objlim=objlim, gain=gain.value,
+            readnoise=readnoise.value, satlevel=satlevel, pssl=pssl,
             niter=niter, sepmed=sepmed, cleantype=cleantype,
             fsmode=fsmode, psfmodel=psfmodel, psffwhm=psffwhm,
             psfsize=psfsize, psfk=psfk, psfbeta=psfbeta, verbose=verbose)
 
         # create the new ccd data object
         nccd = ccd.copy()
+
+        # Remove the gain scaling if it wasn't desired
+        if not gain_apply and gain != 1.0:
+            cleanarr = cleanarr / gain.value
+
+        # Fix the units if the gain is being applied.
+        nccd.unit = ccd.unit * gain.unit
+
         nccd.data = cleanarr
         if nccd.mask is None:
             nccd.mask = crmask

ccdproc/tests/test_cosmicray.py

@@ -6,7 +6,7 @@
 import pytest
 from astropy.utils import NumpyRNGContext
 from astropy.nddata import StdDevUncertainty
-
+from astropy import units as u
 
 from ..core import (cosmicray_lacosmic, cosmicray_median,
                     background_deviation_box, background_deviation_filter)
@@ -64,6 +64,122 @@ def test_cosmicray_lacosmic_check_data():
         cosmicray_lacosmic(10, noise)
 
 
+@pytest.mark.parametrize('array_input', [True, False])
+@pytest.mark.parametrize('gain_correct_data', [True, False])
+def test_cosmicray_gain_correct(array_input, gain_correct_data):
+    # Add regression check for #705 and for the new gain_correct
+    # argument.
+    # The issue is that cosmicray_lacosmic gain-corrects the
+    # data and returns that gain corrected data. That is not the
+    # intent...
+    ccd_data = ccd_data_func(data_scale=DATA_SCALE)
+    threshold = 5
+    add_cosmicrays(ccd_data, DATA_SCALE, threshold, ncrays=NCRAYS)
+    noise = DATA_SCALE * np.ones_like(ccd_data.data)
+    ccd_data.uncertainty = noise
+    # No units here on purpose.
+    gain = 2.0
+    # Don't really need to set this (6.5 is the default value) but want to
+    # make lack of units explicit.
+    readnoise = 6.5
+    if array_input:
+        new_data, cr_mask = cosmicray_lacosmic(ccd_data.data,
+                                               gain=gain,
+                                               gain_apply=gain_correct_data)
+    else:
+        new_ccd = cosmicray_lacosmic(ccd_data,
+                                     gain=gain,
+                                     gain_apply=gain_correct_data)
+        new_data = new_ccd.data
+        cr_mask = new_ccd.mask
+    # Fill masked locations with 0 since there is no simple relationship
+    # between the original value and the corrected value.
+    orig_data = np.ma.array(ccd_data.data, mask=cr_mask).filled(0)
+    new_data = np.ma.array(new_data.data, mask=cr_mask).filled(0)
+    if gain_correct_data:
+        gain_for_test = gain
+    else:
+        gain_for_test = 1.0
+
+    np.testing.assert_allclose(gain_for_test * orig_data, new_data)
+
+
+def test_cosmicray_lacosmic_accepts_quantity_gain():
+    ccd_data = ccd_data_func(data_scale=DATA_SCALE)
+    threshold = 5
+    add_cosmicrays(ccd_data, DATA_SCALE, threshold, ncrays=NCRAYS)
+    noise = DATA_SCALE * np.ones_like(ccd_data.data)
+    ccd_data.uncertainty = noise
+    # The units below are the point of the test
+    gain = 2.0 * u.electron / u.adu
+
+    # Since gain and ccd_data have units, the readnoise should too.
+    readnoise = 6.5 * u.electron
+    new_ccd = cosmicray_lacosmic(ccd_data,
+                                 gain=gain,
+                                 gain_apply=True)
+
+
+def test_cosmicray_lacosmic_accepts_quantity_readnoise():
+    ccd_data = ccd_data_func(data_scale=DATA_SCALE)
+    threshold = 5
+    add_cosmicrays(ccd_data, DATA_SCALE, threshold, ncrays=NCRAYS)
+    noise = DATA_SCALE * np.ones_like(ccd_data.data)
+    ccd_data.uncertainty = noise
+    gain = 2.0 * u.electron / u.adu
+    # The units below are the point of this test
+    readnoise = 6.5 * u.electron
+    new_ccd = cosmicray_lacosmic(ccd_data,
+                                 gain=gain,
+                                 gain_apply=True,
+                                 readnoise=readnoise)
+
+
+def test_cosmicray_lacosmic_detects_inconsistent_units():
+    # This is intended to detect cases like a ccd with units
+    # of adu, a readnoise in electrons and a gain in adu / electron.
+    # That is not internally inconsistent.
+    ccd_data = ccd_data_func(data_scale=DATA_SCALE)
+    ccd_data.unit = 'adu'
+    threshold = 5
+    add_cosmicrays(ccd_data, DATA_SCALE, threshold, ncrays=NCRAYS)
+    noise = DATA_SCALE * np.ones_like(ccd_data.data)
+    ccd_data.uncertainty = noise
+    readnoise = 6.5 * u.electron
+
+    # The units below are deliberately incorrect.
+    gain = 2.0 * u.adu / u.electron
+    with pytest.raises(ValueError) as e:
+        cosmicray_lacosmic(ccd_data,
+                           gain=gain,
+                           gain_apply=True,
+                           readnoise=readnoise)
+    assert 'Inconsistent units' in str(e.value)
+
+
+def test_cosmicray_lacosmic_warns_on_ccd_in_electrons(recwarn):
+    # Check that an input ccd in electrons raises a warning.
+    ccd_data = ccd_data_func(data_scale=DATA_SCALE)
+    # The unit below is important for the test; this unit on
+    # input is supposed to raise an error.
+    ccd_data.unit = u.electron
+    threshold = 5
+    add_cosmicrays(ccd_data, DATA_SCALE, threshold, ncrays=NCRAYS)
+    noise = DATA_SCALE * np.ones_like(ccd_data.data)
+    ccd_data.uncertainty = noise
+    # No units here on purpose.
+    gain = 2.0
+    # Don't really need to set this (6.5 is the default value) but want to
+    # make lack of units explicit.
+    readnoise = 6.5
+    new_ccd = cosmicray_lacosmic(ccd_data,
+                                 gain=gain,
+                                 gain_apply=True,
+                                 readnoise=readnoise)
+    assert len(recwarn) == 1
+    assert "Image unit is electron" in str(recwarn.pop())
+
+
 def test_cosmicray_median_check_data():
     with pytest.raises(TypeError):
         ndata, crarr = cosmicray_median(10, thresh=5, mbox=11,

docs/reduction_toolbox.rst

@@ -99,6 +99,21 @@ Use this technique with `~ccdproc.cosmicray_lacosmic`:
 
     >>> cr_cleaned = ccdproc.cosmicray_lacosmic(gain_corrected, sigclip=5)
 
+.. note::
+
+    By default, `~ccdproc.cosmicray_lacosmic` multiplies the image by
+    the gain; prior to version 2.1 it did so without changing the units of
+    the image which could result in incorrect results.
+
+    There are two ways to correctly invoke `~ccdproc.cosmicray_lacosmic`:
+
+    + Supply a gain-corrected image, in units of ``electron``, and set ``gain=1.0``
+      (the default value) in `~ccdproc.cosmicray_lacosmic`.
+    + Supply an image in ``adu`` and set the ``gain`` argument of
+      `~ccdproc.cosmicray_lacosmic` to the appropriate value for your
+      instrument. Ideally, pass in a ``gain`` with units, but if units are
+      omitted the will be assumed to be ``electron/adu``.
+
 median
 ++++++