Cosmic ray removal
Contents
6.3. Cosmic ray removal¶
Almost all images from a CCD will include some number of cosmic rays, charged particles which bombard the Earth’s upper atmosphere. Some of those will make it through the atmosphere and into your detector (the rate of cosmic rays will be much higher for cameras in space). Although the number of cosmic rays is roughly proportional to exposure time, there will be cosmic rays even in bias frames in which the chip is immediately read out.
This notebook explains how to remove cosmic rays from calibration images and science images.
6.3.1. Removal from calibration images¶
The most convenient way to remove cosmic rays from calibration images (bias, dark, and flat images) is to combine them properly. Cosmic rays are, by their nature, random events that will affect different parts of each of the calibration images. A pixel affected by a cosmic ray in one of the dark images, for example, will almost certainly not be affected by a cosmic ray in any of the other dark images.
Combining those images by averaging (to reduce noise as much as possible) and sigma clipping (to exclude extreme pixels in individual images like the one with a cosmic ray) will eliminate the cosmic ray from the combined dark image. An alternative would be to combine the images using a median. A detailed description of each option is discussed in the section on image combination.
The method described below for removing cosmic rays from science images will not work well for removing them from calibration images and is unnecessary because they can be removed by properly combining the images.
6.3.2. Removal from science images¶
One good technique for removing cosmic rays from an image is the LACosmic method originally developed and implemented for IRAF by Pieter G. van Dokkum. The original paper describing the method, which uses the sharp edges of cosmic rays to distinguish them from other sources in the image, is here.
The specific implementation of LACosmic used here is the astropy affiliated
package Astro-SCRAPPY. If you use this code to remove cosmic
rays you should cite both the original paper and
Astro-SCRAPPY (citation details are on its web site). The
code below never directly imports Astro-SCRAPPY because
ccdproc provides a wrapper for it, so we called attention to it here.
6.3.2.1. Very important notes about using LACosmic¶
There are a few things to be aware of before using the LACosmic technique. These are drawn from the advice van Dokkum provides under Notes for Users and the original paper.
The images must be bias and dark subtracted.
The images should be flat fielded, though the technique can be applied without flat fielding.
The images should not have the sky subtracted before detecting the cosmic rays.
The noise level in the image needs to be accurately measured.
The image and the noise have to be in the same units, typically electrons.
If there are pixels that are known to be bad (e.g., hot pixels, pixels identified by
ccdmask) they should be masked out before detecting cosmic rays.
from pathlib import Path
import numpy as np
from matplotlib import pyplot as plt
from astropy.nddata import CCDData
from astropy.nddata import block_replicate
from astropy import units as u
import ccdproc as ccdp
from photutils import detect_sources
from convenience_functions import show_image, display_cosmic_rays
# Use custom style for larger fonts and figures
plt.style.use('guide.mplstyle')
6.3.2.2. Input image¶
The image we will use in this notebook is one of the reduced images from Example 2 in the previous notebooks. It is an image of the field of the exoplanet KELT-16 b, taken with a thermoelectrically-cooled CCD with read noise of 10\(e^-\) and gain of \(1.5~e^-\)/ADU.
ex2_path = Path('example2-reduced')
ccd = CCDData.read(ex2_path / 'kelt-16-b-S001-R001-C084-r.fit')
WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Equatorial coordinate system
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
show_image(ccd, cmap='gray')
WARNING: The following attributes were set on the data object, but will be ignored by the function: uncertainty, mask, meta, unit, wcs [astropy.nddata.decorators]
The unit of this image is ADU, so we need to multiply by the gain to convert to electrons.
ccd = ccdp.gain_correct(ccd, 1.5 * u.electron / u.adu)
6.3.2.3. Reading in and applying masks¶
Two masks were calculated earlier. Hot pixels, whose dark current cannot be
corrected, were identified by comparing dark frames of different exposure time. The
function ccdmask was used to identify other bad pixels; one example was a
column of pixels on the left side of the image.
We read in both of the masks, which are the same for all images, and combine using logical “OR.”
dark_mask = CCDData.read(ex2_path / 'mask_from_dark_current.fits', unit=u.dimensionless_unscaled)
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Input In [6], in <cell line: 1>()
----> 1 dark_mask = CCDData.read(ex2_path / 'mask_from_dark_current.fits', unit=u.dimensionless_unscaled)
File /usr/share/miniconda/envs/test/lib/python3.8/site-packages/astropy/nddata/mixins/ndio.py:59, in NDDataRead.__call__(self, *args, **kwargs)
58 def __call__(self, *args, **kwargs):
---> 59 return self.registry.read(self._cls, *args, **kwargs)
File /usr/share/miniconda/envs/test/lib/python3.8/site-packages/astropy/io/registry/core.py:197, in UnifiedInputRegistry.read(self, cls, format, cache, *args, **kwargs)
195 try:
196 ctx = get_readable_fileobj(args[0], encoding='binary', cache=cache)
--> 197 fileobj = ctx.__enter__()
198 except OSError:
199 raise
File /usr/share/miniconda/envs/test/lib/python3.8/contextlib.py:113, in _GeneratorContextManager.__enter__(self)
111 del self.args, self.kwds, self.func
112 try:
--> 113 return next(self.gen)
114 except StopIteration:
115 raise RuntimeError("generator didn't yield") from None
File /usr/share/miniconda/envs/test/lib/python3.8/site-packages/astropy/utils/data.py:271, in get_readable_fileobj(name_or_obj, encoding, cache, show_progress, remote_timeout, sources, http_headers)
266 if is_url:
267 name_or_obj = download_file(
268 name_or_obj, cache=cache, show_progress=show_progress,
269 timeout=remote_timeout, sources=sources,
270 http_headers=http_headers)
--> 271 fileobj = io.FileIO(name_or_obj, 'r')
272 if is_url and not cache:
273 delete_fds.append(fileobj)
FileNotFoundError: [Errno 2] No such file or directory: 'example2-reduced/mask_from_dark_current.fits'
ccdmask_mask = CCDData.read(ex2_path / 'mask_from_ccdmask.fits', unit=u.dimensionless_unscaled)
combined_mask = dark_mask.data | ccdmask_mask.data
show_image(combined_mask, cmap='gray')
Excluding these pixels from cosmic ray detection ensures that only cosmic rays are identified.
The mask is now applied to the image of KELT-16 b.
ccd.mask = combined_mask
6.3.2.4. Running LACosmic¶
The actual invocation of LACosmic is fairly convenient. The key parameters
are readnoise, the read noise, and sigclip, which determines how far above
the background a pixel needs to be to consider it a cosmic ray. There is no
hard-and-fast rule for selecting the proper value of sigclip. In the original
paper a value of 5 is recommended, but for this image it finds several thousand
pixels contaminated by cosmic rays. That is not plausible for an image taken
with a camera 1,000 feet above sea level.
Higher values of sigclip reduce the number of cosmic rays found. The value
used below, 7, seemed to work well for this image, finding a total of roughly 70
pixels that are cosmic rays, and a couple dozen candidate cosmic rays that
extend across multiple pixels.
The function cosmicray_lacosmic from ccdproc returns a new image in which
the mask is True for pixels in which a cosmic ray was detected and False
otherwise. The data in the new image has values in the pixels in which cosmic
rays were identified replaced by interpolating the neighboring pixels.
We will take a look at the cosmic rays identified by LACosmic in a moment.
Expect the code below to take at least a few tens of seconds to run.
%%time
new_ccd = ccdp.cosmicray_lacosmic(ccd, readnoise=10, sigclip=7, verbose=True)
The mask of new_ccd includes both cosmic rays identified by
cosmicray_lacosmic and the mask we applied to ccd above. To get only the
cosmic rays, we set the mask to False for all of those pixels that were masked
before LACosmic ran.
cr_mask = new_ccd.mask
cr_mask[ccd.mask] = False
The sum of the mask indicates how many pixels have been identified as cosmic rays.
new_ccd.mask.sum()
6.3.2.5. Examining the cosmic rays identified by LACosmic¶
There are 70 pixels that have been flagged as cosmic rays. Looking through each of them individually would be tedious, at best. It would also presumably be difficult to decide visually if a single pixel tagged as a cosmic ray was actually a cosmic ray, but it would be helpful to look at the larger cosmic rays (i.e., those that span multiple pixels).
To identify those larger cosmic rays we will use the function detect_sources
from the package photutils, which identifies contiguous
pixels in an image via image segmentation. Though
detect_sources is intended for detecting extended or stellar sources in an image it happens to work very well for identifying extended cosmic rays in the mask generated by cosmicray_lacosmic.
The threshold below should be something less than 1 to ensure that only the masked pixels (i.e., those whose values are 1) are included as sources. The number of pixels is the number which must be adjacent (either by edge or by corner) to be considered a source.
threshold = 0.5
n_pixels = 3
crs = detect_sources(new_ccd.mask, threshold, n_pixels)
We first check to see how many of the pixels identified by LACosmic are part of these extended cosmic rays.
crs.areas
It looks like about 50% of the pixels marked as cosmic rays are extended over multiple adjacent pixels.
In this particular science image there are three things being identified as cosmic rays:
Actual cosmic rays.
Single hot pixels (i.e., pixels with unusually high dark current).
These conclusions are not at all clear from what we have discussed in this notebook so far. They are based on a detailed examination of the images after looking at thumbnail views of the cosmic ray mask, the science image in which the cosmic rays were detected and the combined dark frame that was used to calibrate this science image.
That thumbnail view turned out to be a useful enough comparison that a function
called display_cosmic_rays is provided in convenience_functions.py that will
display the cosmic ray mask and as many additional comparison images as you
would like.
Since one of the images that turns out to be useful to look at in this example is the combined dark used in calibrating the science image, we load it.
ccd_dark = CCDData.read(ex2_path / 'combined_dark_90.000.fit')
One note about the argument only_display_rays below, which restricts the
cosmic rays that are displayed to that list. The list was chosen by first
viewing all of the cosmic rays and choosing a representative sample for
inclusion in this discussion. Display them all by setting
only_display_rays=None in the call to display_cosmic_rays.
images_to_display = [new_ccd.mask, ccd, ccd_dark]
image_titles = ['Mask', 'Science image', 'Combined dark']
display_cosmic_rays(crs, images_to_display, titles=image_titles,
only_display_rays=[0, 1, 14, 18]
)
6.3.2.6. Discussion of sample cosmic rays¶
The first three examples above, labeled “Cosmic ray 0,” “Cosmic ray 1,” and “Cosmic ray 14,” are clear-cut; each is actually a cosmic ray.
The fourth, “Cosmic ray 18,” is not a cosmic ray, though it does correspond to a defect in the CCD. It is caused by a single hot pixel (dark current around 2\(e^-\)/sec) that has a high value in the combined dark frame. When that combined dark is subtracted from the science image it causes a large negative value in the value of the science image which ends up being identified as a cosmic ray.
6.3.2.7. Saving the mask with the image¶
To save the full mask, including cosmic rays, hot pixels, and pixels identified
by ccdmask, set the mask of ccd to the mask of new_ccd. In some use cases
you might prefer to save new_ccd itself. The difference between the two is
that the pixel values in which there are cosmic rays have been replaced in
new_ccd by values representative of the surrounding pixels.
ccd.mask = new_ccd.mask
# This saves both the image and the mask
ccd.write('example-with-cosmic-rays.fits')
6.3.3. What happens if you do not mask?¶
In the example above we masked the pixels known to be bad in all of the images before detecting cosmic rays. It is possible to do cosmic ray detection without prior masking. The downside of not masking is that many features that are not cosmic rays are identified as cosmic rays, and some real cosmic rays are not detected.
We begin with a fresh copy of the input image and run cosmic ray detection with the same parameters as above.
ccd = CCDData.read(ex2_path / 'kelt-16-b-S001-R001-C084-r.fit')
%%time
new_ccd_no_premask = ccdp.cosmicray_lacosmic(ccd, readnoise=10, sigclip=7, verbose=True)
print("Pre-masking detects {} cosmic ray pixels.\nNo pre-masking detects {} cosmic ray pixels.".format(new_ccd.mask.sum(), new_ccd_no_premask.mask.sum()))
Many of the additional pixels detected as cosmic rays are in the leftmost column of the image. The column is actually bad (it is covered on the CCD and receives no light).
crs_no_premask = detect_sources(new_ccd_no_premask.mask, threshold, n_pixels)
crs_no_premask.areas
6.3.3.1. Some of these are not cosmic rays¶
The display below shows a few examples of areas identified by
cosmicray_lacosmic that are not actually cosmic rays. These false positives
are harmless because they reflect real problems with the detector and should be
masked out anyway.
More problematic are the cosmic rays that are undetected if the image is not masked first. As an example, the cosmic ray labeled “Cosmic ray 0” in the example above in which a mask was applied before detecting cosmic rays is not detected at all when no masking is done.
images_to_display = [new_ccd_no_premask.mask, ccd, ccd_dark]
image_titles = ['Mask', 'Science image', 'Combined dark']
display_cosmic_rays(crs_no_premask, images_to_display, titles=image_titles,
only_display_rays=[0, 1, 2]
)
6.3.3.2. Discussion¶
The first example above is the bad column on the left side of the image.
The one label “Cosmic ray 1” is caused by a single hot pixel with a large count in the dark frames, which leads to a large negative value in the calibrated science image. LACosmic tags that pixel and many around it as cosmic rays. While the individual hot pixel should be masked, the ones around it do not need to be masked.
The final example, “Cosmic ray 2,” looks a lot like a star. It is, in fact, the after-image of one of the bright stars in the field of KELT-16 b. Bright stars in the field of view can deposit enough charge that it does not dissipate between images. Typically the effect can be avoided altogether by “pre-flashing” the CCD.
Since you cannot pre-flash after the images have been taken, we are stuck doing the next best thing: masking out this part of the images. The masking should be done during the stage at which hot pixels are being identified because all of these pixels will be identified as hot.