6.2. Creating an image mask#

Calibration cannot compensate for every defect in a CCD. Some examples (a non-exhaustive list):

  • Some hot pixels are not actually linear with exposure time.

  • Some pixels in the CCD may respond less to light than others in a way that flat frames cannot compensate for.

  • There may be defects in all or part of a row or column of the chip.

  • Cosmic rays strike the CCD during every exposure. While those are eliminated in the combined calibrated frames with the proper choice of combination parameters, they are not removed from science images.

The first three points are discussed in this notebook. Removal of cosmic rays from science images is discussed in the cosmic ray notebook.

from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt

from astropy import units as u
from astropy.nddata import CCDData

import ccdproc as ccdp

from convenience_functions import show_image, image_snippet

# Use custom style for larger fonts and figures
plt.style.use('guide.mplstyle')

6.2.1. Detecting bad pixels with ccdmask#

The ccdproc function ccdmask uses a method that is based on the IRAF task ccdmask. The method works best when the input image used to detect flaws in the CCD is the ratio of two flat frames with different counts. That may or may not be available depending on what images are collected.

In the example below, which uses images from Example 2 in the reduction notebooks, the two extreme exposure times are 1 sec and 1.2 sec, but the average counts in the images differ by 10,000. These were twilight flats taken just after sunset.

Even with dome flats where the illumination is supposed to be constant, the counts may actually vary. If they do not, use a single flat for identifying bad pixels instead of a ratio.

We begin by creating an image collection and then the information for all of the calibrated, uncombined, flat images.

ex2_path = Path('example2-reduced')

ifc = ccdp.ImageFileCollection(ex2_path)

for long_values in ['history', 'comment']:
    try:
        ifc.summary.remove_column(long_values)
    except KeyError:
        # These two columns were not present, so removing them failed.
        # Just keep going.
        pass

flats = (ifc.summary['imagetyp'] == 'FLAT') & (ifc.summary['combined'] != True)
ifc.summary[flats]
Table masked=True length=0
filesimplebitpixnaxisnaxis1naxis2date-obsexptimeexposureset-tempccd-tempxpixszypixszxbinningybinningxorgsubfyorgsubfreadoutmimagetypfocallenaptdiaaptareaswcreateswserialsitelatsitelongtelescopinstrumenotesflipstatcstretchcblackcwhitepedestalswownerpurgedlatitudelongitudaltitudelstjd-obsmjd-obsbiassectrimsecbunittrim_imagetrimimbscalebzeroextendcombined
str30boolint64int64int64int64str19float64float64float64float64float64float64int64int64int64int64str10str4float64float64float64str34str32str13str14str1str9str1str1str6int64int64int64str10boolstr13str14float64str13float64float64str11str11str3str6str13objectobjectobjectobject

The best we can do here is the ratio of the first and last of the flat images listed above.

first = ifc.summary['file'][flats][0]
last = ifc.summary['file'][flats][-1]

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[5], line 1
----> 1 first = ifc.summary['file'][flats][0]
      2 last = ifc.summary['file'][flats][-1]

File ~/miniconda3/envs/mamba/lib/python3.11/site-packages/astropy/table/_column_mixins.pyx:92, in astropy.table._column_mixins._MaskedColumnGetitemShim.__getitem__()

File ~/miniconda3/envs/mamba/lib/python3.11/site-packages/astropy/table/_column_mixins.pyx:62, in astropy.table._column_mixins.base_getitem()

File ~/miniconda3/envs/mamba/lib/python3.11/site-packages/astropy/table/_column_mixins.pyx:86, in astropy.table._column_mixins.masked_column_getitem()

File ~/miniconda3/envs/mamba/lib/python3.11/site-packages/numpy/ma/core.py:3305, in MaskedArray.__getitem__(self, indx)
   3295 """
   3296 x.__getitem__(y) <==> x[y]
   3297 
   3298 Return the item described by i, as a masked array.
   3299 
   3300 """
   3301 # We could directly use ndarray.__getitem__ on self.
   3302 # But then we would have to modify __array_finalize__ to prevent the
   3303 # mask of being reshaped if it hasn't been set up properly yet
   3304 # So it's easier to stick to the current version
-> 3305 dout = self.data[indx]
   3306 _mask = self._mask
   3308 def _is_scalar(m):

File ~/miniconda3/envs/mamba/lib/python3.11/site-packages/numpy/ma/core.py:3305, in MaskedArray.__getitem__(self, indx)
   3295 """
   3296 x.__getitem__(y) <==> x[y]
   3297 
   3298 Return the item described by i, as a masked array.
   3299 
   3300 """
   3301 # We could directly use ndarray.__getitem__ on self.
   3302 # But then we would have to modify __array_finalize__ to prevent the
   3303 # mask of being reshaped if it hasn't been set up properly yet
   3304 # So it's easier to stick to the current version
-> 3305 dout = self.data[indx]
   3306 _mask = self._mask
   3308 def _is_scalar(m):

IndexError: index 0 is out of bounds for axis 0 with size 0

ccd1 = CCDData.read(ex2_path / first)
ccd2 = CCDData.read(ex2_path / last)

ratio = ccd2.divide(ccd1)

The ratio is roughly 0.85:

ratio.data.mean()

Running ccdmask takes a little time but only needs to be done once, not once for each image.

%%time
maskr = ccdp.ccdmask(ratio)

The result of ccdmask is one where there is a defect and zero where the chip is good, which matches the format of the mask NumPy expects.

The input image and derived mask are shown below:

fig, axes = plt.subplots(1, 2, figsize=(20, 10))

show_image(ratio, cmap='gray', fig=fig, ax=axes[0], show_colorbar=False)
axes[0].set_title('Ratio of two flats')

show_image(maskr, cmap='gray', fig=fig, ax=axes[1], show_colorbar=False, percl=99.95)
axes[1].set_title('Derived mask')

Two comments are in order:

  • The “starfish” pattern in the first image is an artifact of the camera shutter. Ideally, a longer exposure time would be used for the flats to avoid this.

  • It appear at first glance that there were no pixels masked. The problem is that the masked regions are very small and, at the scale shown, happen to not be visible.

Two defects in this CCD are shown below. The first is a small patch of pixels that are vastly less sensitive than the rest. The second is a column on the left edge of the CCD. It turns out this column is not actually exposed to light. ccdmask correctly identifies both patches as bad.

fig, axes = plt.subplots(2, 2, figsize=(10, 10))

width = 100
center = (3823, 2446)
plot_row = 0

image_snippet(ccd1, center, width=width, fig=fig, axis=axes[plot_row, 0])
axes[plot_row, 0].set_title('Flat, camera defect')

image_snippet(maskr, center, width=width, fig=fig, axis=axes[plot_row, 1], is_mask=True)
axes[plot_row, 1].set_title('Mask, same center')

center = (0, 2048)
plot_row = 1

image_snippet(ccd1, center, width=width, fig=fig, axis=axes[plot_row, 0], percu=99.9, percl=70)
axes[plot_row, 0].set_title('Flat, bad column')

image_snippet(maskr, center, width=width, fig=fig, axis=axes[plot_row, 1], is_mask=True)
axes[plot_row, 1].set_title('Mask, same center')

6.2.1.1. Saving the mask#

The mask can be saved in a FITS file as an image. We will see in the summary notebook on masking how to combine the mask generated here with a mask generated from the dark current and with a cosmic ray mask for each science image.

mask_as_ccd = CCDData(data=maskr.astype('uint8'), unit=u.dimensionless_unscaled)
mask_as_ccd.header['imagetyp'] = 'flat mask'
mask_as_ccd.write(ex2_path / 'mask_from_ccdmask.fits')

6.2.2. Making the mask with a single flat#

The flats we used in Example 1, taken with the Large Format Camera at Palomar, are dome flats taken with nearly constant illumination. In that case the best we can do is run ccdmask on a single flat image. As we will see, this still allows the identification of several clearly bad areas of the chip.

First, a look at the calibratted, but not combined, flat images.

ex1_path = Path('example1-reduced')

ifc1 = ccdp.ImageFileCollection(ex1_path)

flats = (ifc1.summary['imagetyp'] == 'FLATFIELD') & (ifc1.summary['combined'] != True)
ifc1.summary[flats]

We can double check that a ratio of flats will not be useful by calculating the mean counts in each flat image:

ccs = []

for c in ifc1.ccds(imagetyp='flatfield', filter="g'"):
    if 'combined' in c.header:
        continue
    print(c.data.mean())
    ccs.append(c)

The variation in counts is so small that the ratio of two flats will not be useful.

Instead, we run ccdmask on the first flat. There is nothing special about that one. The kind of defects that ccdmask tries to identify are in the CCD sensor itself and should be the same for all filters.

%%time
ccs1_mask = ccdp.ccdmask(ccs[0])

Displaying the flat we used and the mask side by side demonstrates that the defects which are clear in the flat are picked up in the mask.

fig, axes = plt.subplots(1, 2, figsize=(15, 10))

show_image(ccs[0], cmap='gray', fig=fig, ax=axes[0])
axes[0].set_title('Single calibrated flat')

show_image(ccs1_mask, cmap='gray', fig=fig, ax=axes[1], is_mask=False)
axes[1].set_title('Derived mask');

A couple of cutouts are shown below illustrating some of the individual defects identified.

fig, axes = plt.subplots(2, 2, figsize=(10, 10))

width = 300
center = (512, 3976)
plot_row = 0

image_snippet(ccs[0], center, width=width, fig=fig, axis=axes[plot_row, 0])
axes[plot_row, 0].set_title('Flat, partial bad column')

image_snippet(ccs1_mask, center, width=width, fig=fig, axis=axes[plot_row, 1], is_mask=True)
axes[plot_row, 1].set_title('Mask, same center')

center = (420, 3250)
width = 100
plot_row = 1

image_snippet(ccs[0], center, width=width, fig=fig, axis=axes[plot_row, 0])
axes[plot_row, 0].set_title('Flat, bad patch')

image_snippet(ccs1_mask, center, width=width, fig=fig, axis=axes[plot_row, 1], is_mask=True)
axes[plot_row, 1].set_title('Mask, same center')