Lab 1.0: Flats, Darks, and Biases#

Due: Wednesday, April 16 by 5:00 pm.

This lab walks through calibration of raw images from the telescope. By the end of this lab, you will be able to

  • perform image calibration using darks, flats, and biases

  • analyze the performance of the camera and telescope using these files

You should then stack these calibrated science images and produced final coadded images for both M35 and M3 in B and V filters on both telescopes (eight images in total).

Each team will submit one PDF of their notebook together. Please make sure all the code is visible.

# If using Colab, mount your Google Drive to access data in the shared folder
from google.colab import drive
drive.mount('/content/drive')
# If using Colab, need to install the reproject package
%pip install reproject
# Import packages!

# numpy is a core package for numerical computing in Python, mostly it does fast array operations
import numpy as np

# matplotlib is a common plotting library (plotly is another good one)
import matplotlib.pyplot as plt
%matplotlib inline

# astropy is useful for reading FITS files, doing coordinate transformations, converting units, and much more
import astropy
import astropy.stats
from astropy.io import fits
from astropy.wcs import WCS # World Coordinate System

# reproject allows us to reproject images onto the same coordinates before stacking
from reproject import reproject_interp

# glob and os are useful for navigating your file system
import glob
import os

# miscellaneous
from tqdm.notebook import tqdm
import warnings
from collections import defaultdict

Calibrating Astronomical Images#

Before we jump into the code, here is a quick primer on calibrations for any charge coupled device (or CCD) camera and why they are necessary. The text has been lifted heavily from Howell (2006, chapter 4).

Bias#

Every camera offsets the zero level so that the A/D converter never outputs a negative number. This offset is different for every camera and every camera setting (for those cameras with changeable settings). This “bias” value is not necessarily constant for every pixel and must be measured so that it can be subtracted from the raw science frame.

The bias image has an exposure time of zero seconds. The shutter remains closed and the CCD is simply read out. The purpose of a bias or zero frame is to allow the user to determine the underlying noise level within each data frame. The bias value in a CCD image is usually a low spatial frequency variation throughout the array, caused by the CCD on-chip amplifiers. This variation should remain constant with time. The rms value of the bias level is the CCD read noise. A bias frame contains both the DC offset level (overscan) and the variations on that level. The nature of the bias variations for a given CCD are usually column-wise variations, but may also have small row-wise components as well. Thus, a 2-D, pixel-by-pixel subtraction is often required. A single bias frame will not sample these variations well in a statistical fashion, so an average bias image of 10 or more single bias frames is recommended.

Dark#

Photons are emitted from every body with a non-zero temperature (on the Kelvin scale). Therefore everything emits light, even your camera. These thermal photons can cause a signal in a sensitive CCD camera and therefore the rate at which this occurs needs to be accounted for in a fully calibrated image.

CCD dark frames are images taken with the shutter closed but for some time period, usually equal to that of your target frames. That is, if one is planning to dark correct a 45 second exposure, a 45 second dark frame would typically be obtained. Longer dark frames can often be avoided using the assumption that the dark current increases linearly with time and a simple scaling can be applied. However, this is not always true. Dark frames are a method by which the thermal noise (dark current) in a CCD can be measured. They also can give you information about bad or “hot” pixels that exist as well as provide an estimate of the rate of cosmic ray strikes at your observing site. Observatory class CCD cameras are usually cooled with liquid nitrogen to temperatures at which the dark current is essentially zero. Many of these systems therefore do not require the use of dark exposure CCD frames in the calibration process. Thermoelectrically cooled systems are typically not cooled to low enough temperatures such that one may ignore the dark current, but they are getting better. In addition, these less expensive models often have poor temperature stability allowing the dark current to wander a bit with time. Multiple darks averaged together are the best way to produce the final dark calibration frame. Note that the bias is also present in dark frames. To get an accurate measure of the dark current, the bias must be subtracted.

Flat#

Not every pixel in a CCD camera responds to the sky in the exact same way. There are many reasons for this. For one, not every pixel is created equally, so there may be some intrinsic sensitivity difference from pixel to pixel. However, the CCD may be obscured slightly by the optics of the telescope, or by dust on the camera or filters. Therefore this pixel-to-pixel variation must also be accounted for.

Flat field exposures are used to correct for pixel-to-pixel variations in the CCD response as well as any nonuniform illumination of the detector itself. Flat fields expose the CCD to light from either a dome screen, the twilight sky, the nighttime sky, or a projector lamp in an attempt to provide a high signal-to-noise ratio (SNR) uniformly illuminated calibration image. For narrow-band imaging, flats are very helpful in removing fringing, which may occur in object frames. Flat field calibration frames are needed for each color, wavelength region, or different instrumental setup used in which object frames are to be taken. A good flat should remain constant to about 1%, with 2% or larger changes being indicators of a possible problem. As with the other calibration frames, at least 5 or more flat fields should be taken and averaged to produce the final flat used for image calibration.

# Specify data directory, and which object we want from which date
data_dir = '/content/drive/MyDrive/' # add a shortcut to PHYSICS100_S2025 in your Google Drive for easy access
object_name = # TODO
date_string = # TODO
telescope_name = # TODO

# Gather all the filenames for the object and dates
filenames = {k: glob.glob(os.path.join(data_dir, 'calibration', telescope_name, date_string, f'{k}*.fit')) for k in ['flat', 'dark', 'bias']} # calibration exposures
filenames['light'] = glob.glob(os.path.join(data_dir, object_name, telescope_name, date_string, '*.fit')) # raw science exposures
# Sort the data and headers into dictionaries (defaultdicts are like dictionaries, but they create new keys automatically)
data = defaultdict(lambda: defaultdict(lambda: defaultdict(list)))
headers = defaultdict(lambda: defaultdict(lambda: defaultdict(list)))

# For each type of exposure, load all the files and append them to the appropriate lists
for k in filenames.keys():
    for fn in filenames[k]:
        with fits.open(fn) as hdul:
            filter = hdul[0].header['FILTER'] if 'FILTER' in hdul[0].header else None # get filter from header
            exp_time = hdul[0].header['EXPTIME'] # get exposure time from header
            data[k][filter][exp_time].append(hdul[0].data) # append data to the appropriate list
            headers[k][filter][exp_time].append(hdul[0].header) # append header to the appropriate list

Question: Why do we need all these exposure times for darks? Why do biases not have a filter or exposure time?

Answer:

# Compute average bias
bias_frames = np.array(data['bias'][None][0.0]) # Bias has no filter and 0.0s exposure time
avg_bias = # TODO: Compute the average bias frame (hint: use np.mean with the axis argument)

Question: What is the shape of the bias_frames array? What does each dimension represent? (Hint: use .shape)

Answer:

# Look at the average bias
plt.imshow(avg_bias, cmap='gray', vmin=, vmax=) # TODO: Set vmin and vmax so that you can see the structure in the image
plt.show()

Question: What is the typical value of the bias? Why does it look noisy? Is there any structure to the noise?

Answer:

# Compute average dark for each exposure time (note: if some exposure times were missing, you would have to scale an existing dark appropriately)
avg_darks = {}
for time, darks in data['dark'][None].items(): # Dark has no filter
    avg_darks[time] = # TODO: compute the average dark and subtract the average bias
# Make a histogram of the dark current at different exposure times
hist_range = # TODO: Set the range of the histogram to ignore the outliers
plt.hist(avg_darks[].flatten(), range=hist_range, bins=50, histtype='step', label='') # TODO: plot a histogram for each exposure time
plt.legend()
plt.show()

Question: What is the typical value of the dark current and do you see any difference between exposure times? What does this tell you? What do you think the outliers are?

Answer:

# Compute mean flat for each filter
mean_flats = {}

for filter in data['flat'].keys():
    all_normed_flats = []
    for time, flats in data['flat'][filter].items():
        corrected_flats = # TODO: subtract average bias and corresponding dark from each flat
        normed_flats = # TODO: normalize flats so that they have a sigma-clipped mean value of 1.0 (hint: use astropy.stats.sigma_clipped_stats)
        all_normed_flats.append(normed_flats)
    all_normed_flats = np.concatenate(all_normed_flats, axis=0)
    mean_flats[filter] = # TODO: compute the sigma-clipped mean flat for each filter

Question: Why do we use a sigma-clipped mean for the flats?

Answer:

fig, axes = plt.subplots(ncols=3, figsize=(15,5))
# TODO: Plot the mean flat for each filter, like we did for the bias
plt.show()

Question: What do you see in the flats? Explain what causes this and why it has that shape.

Answer:

# Calibrate all science images
calibrated_data = defaultdict(lambda: defaultdict())

for filter, dict in data['light'].items():
    for time, lights in dict.items():
        dark = # TODO: get the right dark for this exposure time
        flat = # TODO: get the right flat for this filter
        calibrated_lights = # TODO: correct for bias, dark, and flat (hint: two are subtracted, one is divided)
        calibrated_data[filter][time] = calibrated_lights

Question: In the cell above, which are subtracted, which are divided, and why?

Stack your images#

We’ll leave this section up to you. You should be able to use the code from Lab 0.2. You should plot and save your coadded images for M35 and M3 in B and V filters.

# TODO: stack your images
# TODO: visually inspect your coadded images
# TODO: save your coadded images

Acknowledgements#

TODO: Please write a short team contribution statement together explaining what each member did for the lab!