Image pre-reduction

Image pre-reduction#

This page will take us through the full pre-reduction of the science images.

You could copy or write the code shown below in a script, or alternatively directly download this page as a jupyter notebook file.

Before running the code, make sure you’ve created and adjusted the module dataredconfig.py, as explained here.

import dataredconfig

from pathlib import Path
import numpy as np
import astropy
from astropy import units as u
%matplotlib ipympl
import matplotlib
from matplotlib import pyplot as plt
import ccdproc
import photutils.background

print("numpy", np.__version__, ", astropy", astropy.__version__, ", matplotlib", matplotlib.__version__, ", ccdproc", ccdproc.__version__, ", photutils", photutils.__version__)
# Developed with (Feb 2025): numpy 2.2.2 , astropy 7.0.1 , matplotlib 3.10.0 , ccdproc 2.4.3 , photutils 2.1.0

# We'll ignore some astropy warnings that get raised as our FITS headers (from NINA) are not 100% standards compliant.
import warnings
from astropy.utils.exceptions import AstropyWarning
warnings.simplefilter('ignore', category=AstropyWarning)

Creating a master-bias#

bias_files = ccdproc.ImageFileCollection(dataredconfig.data_dir / "BIAS", keywords=dataredconfig.ifc_header_keywords)
bias_files.summary

# We check the level and noise in these files:
for ccd, filename in bias_files.ccds(ccd_kwargs={'unit': 'adu'}, return_fname=True):
    print(f"{filename}: standard deviation {ccd.data.std():.2f} ADU, median level {np.median(ccd.data):.2f} ADU")

Question

We’ll now build our master-bias by averaging over the bias frames. Assuming that the bias frames contain pure uncorrelated noise (in addition to a prefectly constant bias level), what pixel-to-pixel standard deviation would you expect to observe in the master-bias, and why? Compare this with the actual measurement.

combiner = ccdproc.Combiner(bias_files.ccds(ccd_kwargs={"unit":"adu"}))
masterbias = combiner.average_combine()

masterbias.meta['combined'] = True
masterbias.data = masterbias.data.astype('float32') # Converts to float32 to save space
masterbias.write(dataredconfig.work_dir / 'masterbias.fits', overwrite=True)

print(f"masterbias: standard deviation {masterbias.data.std():.2f} ADU, median level {np.median(masterbias.data):.2f} ADU")

Let’s now inspect this master bias visually, to see if it has any structure.

With the code below, you could zoom around in the masterbias. Does it look homogeneous? We also create a smoothed version of the masterbias, to highlight a bit some (potential) larger scale structures. We will of course not use this smoothed version in the pre-reduction, but it might help to visualize some of the structure.

Note

When displaying any image, always make sure to adjust the so called greyscale (or “brightness” scale, sometimes also called “z-scale” in astronomical software) to properly see whatever structure you want to inspect. In the present case, add the keyword arguments vmin=?, vmax=? (with properly chosen values instead of ?, based for example on the median level computed above) to the call of plt.imshow(). These are “cuts” for a linear grey level mapping: pixels with value vmin (or below) will be shown as black, and those with value vmax (or above) will be shown as white, and everything in between will be mapped to some grey level. If you don’t do this, these “cuts” will be automatically selected by some algorithm (often as the minimum and maximum values of the image), and might not be adequate to visualize the details that you want to inspect.

# Median-smoothing the masterbias (only for display purposes):
smooth_masterbias = photutils.background.Background2D(masterbias.data, (100, 100), filter_size=(3, 3), bkg_estimator=photutils.background.MedianBackground()).background

plt.figure(figsize=(8, 4))
plt.title("Masterbias")
cbar = plt.imshow(masterbias, origin='lower', cmap='Greys_r', interpolation='nearest')
plt.colorbar(cbar, label="Pixel value [ADU]")
plt.xlabel("x [pixel]")
plt.ylabel("y [pixel]")
plt.tight_layout()
plt.show()

Creating a master-dark#

dark_files = ccdproc.ImageFileCollection(dataredconfig.data_dir / "DARK", keywords=dataredconfig.ifc_header_keywords)

# See below if such a filter is needed:
#dark_files = dark_files.filter(exptime=60.0)

dark_files.summary

# For this simple example we want all darks to have the same exptime.
exptimes = list(set(dark_files.summary["exptime"]))
assert(len(exptimes) == 1) # If this fails, apply some selection to your files.
dark_exptime = exptimes[0]
print("Exposure times of darks: ", dark_exptime)

# As done for the bias frames, let's check the level and noise in these darks:

for ccd, filename in dark_files.ccds(ccd_kwargs={'unit': 'adu'}, return_fname=True):
    print(f"{filename}: standard deviation {ccd.data.std():.2f} ADU, median level {np.median(ccd.data):.2f} ADU")

Question

Why is the standard deviation on these dark frames much larger than for bias frames, while the level is almost the same?

# Create the masterdark

dest_dir = dataredconfig.work_dir / "DARK_BIASSUB"
dest_dir.mkdir(exist_ok=True)

# First loop over the dark frames to subtract the bias: 
for ccd, filename in dark_files.ccds(ccd_kwargs={'unit': 'adu'}, return_fname=True):
    print(f"Processing {filename}...")
    ccd = ccdproc.subtract_bias(ccd, masterbias)
    # Write to disk:
    ccd.data = ccd.data.astype('float32')
    ccd.write(dest_dir / filename, overwrite=True)

# Note: depending on the available memory,
# it might be possible to do this without writing intermediary files,
# following something like:
#dark_ccds = [ccdproc.subtract_bias(ccd, masterbias) for ccd in dark_files.ccds(ccd_kwargs={'unit': 'adu'})]

# Then we combine the files:
files_to_combine = ccdproc.ImageFileCollection(dest_dir).files_filtered(include_path=True)
masterdark = ccdproc.combine(files_to_combine,
                                 method='average',
                                 sigma_clip=True, sigma_clip_low_thresh=5, sigma_clip_high_thresh=5,
                                 sigma_clip_func=np.ma.median, signma_clip_dev_func=astropy.stats.mad_std,
                                 mem_limit=1e9
                                )

masterdark.meta['combined'] = True
masterdark.data = masterdark.data.astype('float32') # Converts to float32 to save space
masterdark.write(dataredconfig.work_dir / 'masterdark.fits', overwrite=True)

Same as for the bias: let’s check this for some structure. The median-filtering helps to hide the hot pixels and highlight the larger spatial scales instead. Display the plain masterdark to see the hot pixels. Of course we’ll later use the plain masterdark without any filtering!

smooth_masterdark = photutils.background.Background2D(masterdark.data, (200, 200), filter_size=(1, 1), bkg_estimator=photutils.background.MedianBackground()).background

plt.figure(figsize=(8, 4))
cbar = plt.imshow(smooth_masterdark, origin='lower', cmap='Greys_r', interpolation='nearest')
plt.colorbar(cbar, label="Pixel value [ADU]")
plt.tight_layout()
plt.show()

Note

As you saw in the code, the masterbias has already been subtracted from this masterdark. So in this figure you should see the effect of the dark current alone.

Creating the master-flats#

flat_files = ccdproc.ImageFileCollection(dataredconfig.data_dir / "FLAT", keywords=dataredconfig.ifc_header_keywords)
flat_files.summary

# Scaling function for the flats: they get "normalized" (i.e., divided) by their median:
def inv_median(a):
    return 1 / np.median(a)

# We'll simply loop over these filters:
flat_filters_to_run_on = ["g", "r", "i"]

for selected_filter in flat_filters_to_run_on:

    selected_flat_files = flat_files.filter(filter=selected_filter)
    
    dest_dir = dataredconfig.work_dir / f"FLAT_{selected_filter}_BIASDARKSUB"
    dest_dir.mkdir(exist_ok=True)

    # First subtract bias and dark from every flat:
    for ccd, filename in selected_flat_files.ccds(ccd_kwargs={'unit': 'adu'}, return_fname=True):
        print(f"Processing {filename}...")
        ccd = ccdproc.subtract_bias(ccd, masterbias)
        ccd = ccdproc.subtract_dark(ccd, masterdark, exposure_time='exptime', exposure_unit=u.second, scale=True)
    
        # Write to disk:
        ccd.data = ccd.data.astype('float32') # Converts to float32 to save space
        ccd.write(dest_dir / filename, overwrite=True)


    # And now combine to a masterflat:
    files_to_combine = ccdproc.ImageFileCollection(dest_dir).files_filtered(include_path=True)
    masterflat = ccdproc.combine(files_to_combine,
                                 method='average', scale=inv_median,
                                 sigma_clip=True, sigma_clip_low_thresh=5, sigma_clip_high_thresh=5,
                                 sigma_clip_func=np.ma.median, signma_clip_dev_func=astropy.stats.mad_std,
                                 mem_limit=1e9
                                )

    masterflat.meta['combined'] = True
    masterflat.data = masterflat.data.astype('float32') # Converts to float32 to save space
    masterflat.write(dataredconfig.work_dir / f"masterflat_{selected_filter}.fits", overwrite=True)

# Checking the masterflats

masterflat = ccdproc.CCDData.read(dataredconfig.work_dir / "masterflat_r.fits")
plt.figure(figsize=(8, 4))
cbar = plt.imshow(masterflat.data, vmin=0.9, vmax=1.05, origin='lower', cmap='Greys_r', interpolation='nearest')
plt.colorbar(cbar)
plt.tight_layout()
plt.show()

Processing the science images#

science_files = ccdproc.ImageFileCollection(dataredconfig.data_dir / "LIGHT", keywords=dataredconfig.ifc_header_keywords)
science_files.summary

# Defining the directory where pre-reduced images should be written:
dest_dir = dataredconfig.work_dir / "LIGHT_PRERED"
dest_dir.mkdir(exist_ok=True)

# Reading in from disk the master calibration frames:
masterbias = ccdproc.CCDData.read(dataredconfig.work_dir / "masterbias.fits")
masterdark = ccdproc.CCDData.read(dataredconfig.work_dir / "masterdark.fits")
# For the flats, we build a dictionnary that holds the masterflat for each filter name:
available_flat_filters = ["g", "r", "i"]
masterflats = {filtername: ccdproc.CCDData.read(dataredconfig.work_dir / f"masterflat_{filtername}.fits") for filtername in available_flat_filters}

n_files = len(science_files.summary)
# We loop over all files, selecting the right masterflat on the fly
for (i, (ccd, filename)) in enumerate(science_files.ccds(ccd_kwargs={'unit': 'adu'}, return_fname=True)):
    print(f"=== Processing image {i+1}/{n_files}: {filename} ===")
    ccd = ccdproc.subtract_bias(ccd, masterbias)
    ccd = ccdproc.subtract_dark(ccd, masterdark, exposure_time='exptime', exposure_unit=u.second, scale=True)
    
    masterflat_to_use = masterflats[ccd.header['filter']]
    ccd = ccdproc.flat_correct(ccd, masterflat_to_use)
    
    # Write to disk:
    ccd.data = ccd.data.astype('float32') # Converts to float32 to save space
    ccd.write(dest_dir / filename, overwrite=True)