Calibration frames#
Prerequisites
Basic statistics
Knowledge about CCD, bias, dark, and flat frames.
We will learn
What do bias, dark, and flat frames look like
How to deal with outliers; sigma-clipping.
How to combine individual calibration frames to make “master” calibration frames.
Further familiarize with operations on numpy arrays and data I/O.
Reference
Here is a great tutorial on CCD data reduction. Please check out!!
Howell, S., Handbook of CCD Astronomy, Second Ed, Cambridge University Press 2006
%load_ext autoreload
%autoreload 2
import os, sys
import numpy as np
import matplotlib.pyplot as plt
# We can beautify our plots by changing the matpltlib setting a little
plt.rcParams['font.size'] = 18
plt.rcParams['figure.figsize'] = (8, 6)
plt.rcParams['figure.dpi'] = 80
plt.rcParams['axes.linewidth'] = 2
Show code cell source
# Don't worry about this cell, it's just to make the notebook look a little nicer.
try:
import google.colab
IN_COLAB = True
except ImportError:
IN_COLAB = False
# Get the directory right
if IN_COLAB:
# Download utils.py
!wget -q -O /content/utils.py https://raw.githubusercontent.com/AstroJacobLi/ObsAstGreene/refs/heads/main/book/docs/utils.py
# Function to check and install missing packages
from google.colab import drive
drive.mount('/content/drive/')
os.chdir('/content/drive/Shareddrives/AST207/data')
else:
# if not os.path.exists("../../../_static/ObsAstroData/"):
# os.makedirs("../../../_static/ObsAstroData/")
sys.path.append('../../docs/')
os.chdir('../../_static/ObsAstroData/')
1. Bias frames#
As discussed in class, a small voltage is applied to every pixel to ensure the pixel value is always positive. This voltage introduces a “bias” to the pixel values. However, the bias can vary slightly between pixels and may change over time. To account for this, we take zero-time exposures (where dark current and real photons do not contribute to the pixel values) to characterize the bias.
Below we load bias frames taken using the camera of IFU-M on the 6.5 m Magellan telescope in Chile.
1.1 Read a bias frame#
from astropy.io import fits
FITS (Flexible Image Transport System) format is widely used in astronomy for images and tables. The following code will read the file b0854c1.fits, which is a list containing only one HDU (Header Data Unit).
fits.open('./calibration_frames/bias/b0860c1.fits')
[<astropy.io.fits.hdu.image.PrimaryHDU object at 0x10e7606d0>]
fits.open('./calibration_frames/bias/b0860c1.fits').info()
Filename: ./calibration_frames/bias/b0860c1.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 143 (2176, 1156) int16 (rescales to uint16)
Now we know that this file has one HDU, which is an image containing 2176 * 1156 pixels. The data type in int16. We can check the header of the HDU.
hdu = fits.open('./calibration_frames/bias/b0860c1.fits')[0]
hdu.header
SIMPLE = T
BITPIX = 16
NAXIS = 2
NAXIS1 = 2176
NAXIS2 = 1156
BSCALE = 1.00000000 / real=bzero+bscale*value
BZERO = 32768.00000000
BUNIT = 'DU/PIXEL'
ORIGIN = 'LCO/OCIW'
OBSERVER= 'Bailey' / observer name
TELESCOP= 'Clay' / telescope
SITENAME= 'LCO'
SITEALT = 2405 / meters
SITELAT = -29.01423
SITELONG= -70.69242
TIMEZONE= 4
DATE-OBS= '2023-03-21T10:50:38' / (start)
UT-DATE = '2023-03-21' / UT date (start)
UT-TIME = '10:50:38' / UT time (start) 39038 seconds
UT-END = '10:50:39' / UT time (end) 39039 seconds
LC-TIME = '07:50:38' / local time (start) 28238 seconds
NIGHT = '20Mar2023' / local night
MJD = 60024.45183 / modified Julian day
INSTRUME= 'IFUM' / instrument name
SCALE = 0.250 / arcsec/pixel
EGAIN = 0.680 / electrons/DU
ENOISE = 2.700 / electrons/read
NOPAMPS = 4 / # of op-amps
OPAMP = 1 / this op-amp
CHOFFX = 256.000 / x-offset [arcsec]
CHOFFY = 257.000 / y-offset [arcsec]
DISPAXIS= 1 / disperser axis
RA = ' 17:57:11.9' / right ascension
RA-D = 269.2997917 / degrees
DEC = '-28:59:59.0' / declination
DEC-D = -28.9997222 / degrees
EQUINOX = 2000.00000 / of above coordinates
EPOCH = 2023.21768 / epoch (start)
AIRMASS = 1.000 / airmass (start)
ST = 64960.8 / sidereal time: 18:02:41 (start)
ROTANGLE= 273.69 / rotator offset angle
TEL-AZIM= 90.0 / telescope azimuth
TEL-ELEV= 90.0 / telescope elevation
T-DOME = 13.9 / dome temperature [C]
T-CELL = 12.4 / primary cell temperature [C]
T-TRUSS = 8.7 / telescope truss temperature [C]
WX-TEMP = 13.6 / weather station temperature [C]
FILENAME= 'b0860c1'
OBJECT = '1x2 bias' / object name
COMMENT = / comment
EXPTYPE = 'Bias' / exposure type
EXPTIME = 0.000 / exposure time
NLOOPS = 250 / # of loops per sequence
LOOP = 7 / # within this sequence
BINNING = '1x2' / binning
SPEED = 'Slow' / readout speed
NOVERSCN= 128 / overscan pixels
NBIASLNS= 128 / bias lines
BIASSEC = '[2049:2176,1029:1156]' / NOAO: bias section
DATASEC = '[1:2048,1:1028]' / NOAO: data section
TRIMSEC = '[1:2048,1:1028]' / NOAO: trim section
SUBRASTR= 'none' / full frame
TEMPCCD = -109.9 / CCD temperature [C]
SHOE = 'B'
IFU = 'LSB'
IFUENC = -1801 / IFU encoder
OCCULTOR= 'Science Gap'
OCCENC = 98000 / occultor encoder
SLIDE = 'LoRes'
SLIDEENC= 4061.0 / slide encoder
SLIDESTP= 285126 / slide step
LO-ELEV = '9.1301' / step= -107184
HI-AZIM = '-0.55001?' / step= -14143
HI-ELEV = '0.40001?' / step= 10286
FOCUS = 196.1 / step= 1961
FILTER = 'BK7'
SLITNAME= '80' / 607 526
BENEAR = 0 / lamp level (0..10)
LIHE = 0 / lamp level (0..10)
THXE = 0 / lamp level (0..10)
LAMPS = 'Off' / lamp combo
LED-UV = 0 / LED level (0..4096)
LED-BL = 0 / LED level (0..4096)
LED-VI = 0 / LED level (0..4096)
LED-NR = 0 / LED level (0..4096)
LED-FR = 0 / LED level (0..4096)
LED-IR = 0 / LED level (0..4096)
LEDS = 'Off' / LED combo
CONFIGFL= 'unknown' / configuration file
TEMP01 = 14.500 / IFU_Entrance
TEMP02 = 13.750 / IFU_Top
TEMP03 = 13.812 / Fiber_Exit
TEMP04 = 25.100 / IFU_Motor
TEMP05 = 24.600 / IFU_Drive
TEMP06 = 19.688 / IFU_Hoffman
TEMP07 = 20.562 / IFU_Shoebox
TEMP08 = 13.375 / Cradle_R
TEMP09 = 13.188 / Cradle_B
TEMP10 = 13.562 / Echelle_R
TEMP11 = 13.438 / Echelle_B
TEMP12 = 13.688 / Prism_R
TEMP13 = 13.375 / Prism_B
TEMP14 = 13.250 / LoResGrating_R
TEMP15 = 13.250 / LoResGrating_B
ADC-STGE= 'in' / ADC stage
ADC-CTRL= 'auto' / ADC control
MIRROR = 'in' / Mirror stage
DCU-FFS = 'retracted' / DCU FF-Screen
DCU-LMPS= 'Off' / DCU lamps
SOFTWARE= 'Version beta (0071) (Oct 5 2022, 17:02:27)'
FITSVERS= '0.070' / FITS header version
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
COMMENT
Let’s read the data array.
img = hdu.data
img
array([[817, 817, 817, ..., 816, 813, 819],
[816, 818, 820, ..., 817, 818, 828],
[809, 827, 817, ..., 812, 814, 816],
...,
[813, 810, 815, ..., 814, 816, 823],
[821, 818, 816, ..., 823, 820, 816],
[812, 819, 815, ..., 819, 816, 814]], dtype=uint16)
Exercise 1#
Use the statistical tools that you’ve learned to understand this bias data. This will help you to understand the characteristics of the bias frame.
Where do the pixel values concentrate? How to quantify this?
What’s the dispersion of the pixel values? How to quantify this?
Are there any outliers? How to identify them?
What’s the distribution of the pixel values? How to visualize this?
Check that there is no negative pixel value in the bias frame. Why?
Tip: The dispersion in the pixel values of a bias frame is largely contributed by the read-out noise, a noise that is added to each pixel when reading the charge out. We already know the read-out noise from Basic Stats (part 2). Read-out noise may be viewed as a “toll” that must be paid for reading the stored charge. Thanks to the bias voltage, none of the pixels are negative.
## Your answer here
We provide a function show_image to help visualize astronomical images that typically have a large dynamical range. Please feel free to use it.
from utils import show_image
show_image(img, show_colorbar=True, figsize=(12, 12));
The bias frame shows some spatial patterns: some rows have a lower bias value than others. We also see some distinct “bright” pixels. They could be either bad pixels, or caused by Cosmic Rays (CR). As a consequence, when we plot the histogram of pixel values, we will see some outliers from the bulk of pixels.
During the Magellan observation, we took many bias frames. So let’s load another bias frame. What do you expect? Will it be very different from the one above?
show_image(fits.open('./calibration_frames/bias/b0861c1.fits')[0].data);
Yeah, the bright spots changed their positions in the image. So it does seem like they are caused by cosmic rays.
Exercise 2#
What remains the same between different bias frames? What is different?
How might we use all the bias frames to make the best bias correction?
We have provided lots of bias frames in ./calibration_frames/bias. You are encouraged to load them and explore the data.
## Your answer here
1.2 Combine bias frames to make “master” bias#
We have already learned from basic statistics that combining images will reduce the noise level. Belowe we try to combine 40 bias frames. The bias frames are all named bXXXXc1.fits, where XXXX is a four-digit number, from 0860 to 0900.
biases = [fits.open(f'./calibration_frames/bias/b0{ind}c1.fits')[0].data for ind in range(860, 900)]
biases = np.array(biases)
Here I use “Nested List Comprehensions” to construct a list. This is equivalent to the following snippet, but is much shorter.
biases = []
for ind in range(860, 900):
biases.append(fits.open(
f'./calibration_frames/bias/b0{ind}c1.fits'
)[0].data)
biases.shape
(40, 1156, 2176)
# We have introduced how to calculate the mean and median image using the `axis` argument.
# Check `basic_stat.ipynb`
med_bias = np.median(biases, axis=0)
mean_bias = np.mean(biases, axis=0)
fig, [ax1, ax2] = plt.subplots(2, 1, figsize=(12, 12))
ax1 = show_image(med_bias, ax=ax1, fig=fig,)
ax1.text(0.1, 0.8, 'Median-combined', fontsize=25, color='w',
ha='left', transform=ax1.transAxes)
ax2 = show_image(mean_bias, ax=ax2, fig=fig,)
ax2.text(0.1, 0.8, 'Mean-combined', fontsize=25, color='w',
ha='left', transform=ax2.transAxes)
plt.tight_layout()
Wow! The cosmic rays are everywhere in the mean-combined image! The vast difference between the median-combined and mean-combined images highlights that the mean is far less immune to outliers (i.e., cosmic rays) than the median. Because of this, we often use median to combine images.
The histogram below shows again that combining individual bias frames will help reduce the read-out noise, thus boost the SNR.
plt.hist(biases[0].flatten(), histtype='step', density=True,
range=(795, 840), bins=25, lw=1, label='Single bias frames')
for i in range(1, 21):
plt.hist(biases[i].flatten(), histtype='step', density=True,
range=(795, 840), bins=25, alpha=0.1, color='C0')
plt.hist(med_bias.flatten(), histtype='step', density=True,
range=(795, 840), bins=25, lw=2, label='40 biases combined')
plt.xlabel('Pixel value [counts]')
plt.legend(loc='upper center', ncols=2, frameon=False, bbox_to_anchor=(0.5, 1.15));
1.3 Sigma-clipping#
In astronomy, there is yet another way to deal with outliers.
Recall that for a Gaussian distribution, the probability for a data point being \(>3\sigma\) away from the mean \(\mu\) is very very slim. So if we see such a data point that is quite far from the center of the data set, it is probably an outlier and should be removed. This is the idea of sigma clipping.
In the following, we flag out pixels that are \(>3\sigma\) away from the median value (we don’t use the mean value because it is heavily biased by these cosmic rays). You can see the cosmic rays being marked as black in the figure below (they are tiny dots).
hdu = fits.open('./calibration_frames/bias/b0860c1.fits')[0]
img = hdu.data.astype(float)
med, std = np.median(img), np.std(img)
flag = (np.abs(img - med) > 3 * std)
show_image(flag, cmap='Greys', vmin=0., vmax=0.1, show_colorbar=False);
We can replace these outlier pixels with np.nan (not a number, typically used to indicate a bad value). Then we repeat this for all 40 bias frames and combine them.
# This might be slow
biases_clipped = []
for ind in range(860, 900):
img = fits.open(f'./calibration_frames/bias/b0{ind}c1.fits')[0].data.astype(float)
med, std = np.median(img), np.std(img)
flag = (np.abs(img - med) > 3 * std)
img[flag] = np.nan
biases_clipped.append(img)
biases_clipped = np.array(biases_clipped)
# might be slow
med_bias_clipped = np.nanmedian(biases_clipped, axis=0)
mean_bias_clipped = np.nanmean(biases_clipped, axis=0)
fig, [ax1, ax2] = plt.subplots(2, 1, figsize=(12, 12))
ax1 = show_image(med_bias_clipped, ax=ax1, fig=fig, vmin=816.6, vmax=818.4)
ax1.text(0.1, 0.8, 'Median-combined, sigma-clipped', fontsize=25, color='w',
ha='left', transform=ax1.transAxes)
ax2 = show_image(mean_bias_clipped, ax=ax2, fig=fig, vmin=816.6, vmax=818.4)
ax2.text(0.1, 0.8, 'Mean-combined, sigma-clipped', fontsize=25, color='w',
ha='left', transform=ax2.transAxes)
plt.tight_layout()
Yeah! After sigma-clipping, the mean-combined image looks much better and very similar to the median-combined one. This sigma-clipping trick is quite common in astronomy. We used a threshold of 3-sigma and only clipped once. You can repeat the procedure for several times to make sure all outliers are flagged. Actually, astropy already provides us a function sigma_clip to do these jobs. Let’s check it out.
from astropy.stats import sigma_clip
astropy_clipped = sigma_clip(biases, # the bias frames before sigma-clipping
sigma=3, # threshould is 3-sigma
maxiters=5, # do clipping for 5 times at most
cenfunc='median', stdfunc='std',
# use median and standard deviation to calculate clipping criterion
axis=(1, 2), # calc median and std for each bias frame
masked=False)
# might be slow
med_bias_apy = np.nanmedian(astropy_clipped, axis=0)
mean_bias_apy = np.nanmean(astropy_clipped, axis=0)
fig, [ax1, ax2] = plt.subplots(2, 1, figsize=(12, 12))
ax1 = show_image(med_bias_apy, ax=ax1, fig=fig, vmin=816.6, vmax=818.4)
ax1.text(0.1, 0.8, 'Median-combined, astropy sigma-clipped', fontsize=25, color='w',
ha='left', transform=ax1.transAxes)
ax2 = show_image(mean_bias_apy, ax=ax2, fig=fig, vmin=816.6, vmax=818.4)
ax2.text(0.1, 0.8, 'Mean-combined, astropy sigma-clipped', fontsize=25, color='w',
ha='left', transform=ax2.transAxes)
plt.tight_layout()
We can write the median-combined image to a FITS file for future usage.
hdu = fits.PrimaryHDU(med_bias_clipped.astype(np.float32))
hdulist = fits.HDUList([hdu])
hdulist.writeto('./calibration_frames/bias/master_bias.fits', overwrite=True)
2. Dark frames#
Just to give you a flavor of dark frame, here we show the five dark frames that we took during the Magellan observation. They are full of cosmic rays. We construct a master dark frame by combining them using median.
darks = [fits.open(f'./calibration_frames/dark/b{ind}c1.fits')[0].data for ind in range(1212, 1217)]
darks = np.array(darks)
med_dark = np.median(darks, axis=0)
fig, [ax1, ax2] = plt.subplots(2, 1, figsize=(12, 12))
ax1 = show_image(darks[0], ax=ax1, fig=fig)
ax1.text(0.1, 0.8, 'Single dark frame (3600s)', fontsize=25, color='w',
ha='left', transform=ax1.transAxes)
ax2 = show_image(med_dark, ax=ax2, fig=fig)
ax2.text(0.1, 0.8, 'Master dark frame (median-combined)', fontsize=25, color='w',
ha='left', transform=ax2.transAxes)
plt.tight_layout()
The light area in the lower left corner is due to sensor glow, which is light emitted by the electronics of the CCD. The counts in a dark frame scales linearly with time. The exposure time for a dark frame is typically similar to (or longer than) the exposure time of your science frames.
3. Flat frames#
Flat frames are used to correct for pixel-to-pixel sensitivity variations due to dust on the optics, vignetting, CCD quantum efficiency, etc. They are taken by pointing the telescope at a uniformly illuminated screen. The screen can be a white screen, a twilight sky, or a dome screen. Below, we load the flat frames taken during the Magellan observation.
Let’s load a few flat frames of the Magellan IMACS camera and see how they look like. These data are taken with for the i filter using the “screen” in the dome (so-called “dome flat”).
camera = 'ift'
detector = 'c2'
flats = [fits.open(f'./calibration_frames/flat/{camera}/{camera}{ind:04d}{detector}.fits')[0].data.T for ind in range(14, 19)]
flats = np.array(flats)[:, :1024, :2048] # we trim off the overscan region
med_flat = np.median(flats, axis=0)
fig, [ax1, ax2] = plt.subplots(2, 1, figsize=(12, 12))
ax1 = show_image(flats[0], ax=ax1, fig=fig)
ax1.text(0.1, 0.8, f'Single dome flat frame for {camera}-{detector}', fontsize=25, color='w',
ha='left', transform=ax1.transAxes)
ax2 = show_image(med_flat, ax=ax2, fig=fig)
ax2.text(0.1, 0.8, f'Master flat frame for {camera}-{detector} \n (median-combined)', fontsize=25, color='w',
ha='left', transform=ax2.transAxes)
plt.tight_layout()
For the above figure, we see that the flat frames are not flat at all! The vignetting effect (i.e., the illumination is not uniform across the field) is quite obvious. The absolute value of the pixels in the flat frame is not important, but the relative value is. So we normalize the flat frames by dividing them by their median value.
med_flat_norm = med_flat / np.nanmedian(med_flat)
show_image(med_flat_norm, vmin=0.7, vmax=1.1, cmap='inferno',
show_colorbar=True, figsize=(12, 6), colorbar_label='Normalized pixel value');
Note that the above flat frame is only 1/8 of the full field-of-view of the IMACS camera, but the illumination pattern is still quite obvious. Stars at the edge/corner of the field-of-view will look dimmer than those in the center, despite having the same magnitude/intrinsic brightness. In order to correct for this, we need to divide any science frame by the flat frame.
Not all flat frames look this good! Let’s see the flat frames of another camera (also on Magellan)! You’ll see lots of dust spots and artifacts.
camera = 'iff'
detector = 'c3'
flats = [fits.open(f'./calibration_frames/flat/{camera}/{camera}{ind:04d}{detector}.fits')[0].data.T for ind in range(8, 13)]
flats = np.array(flats)[:, :1024, :2048] # we trim off the overscan region
med_flat = np.median(flats, axis=0)
fig, [ax1, ax2] = plt.subplots(2, 1, figsize=(12, 12))
ax1 = show_image(flats[0], ax=ax1, fig=fig)
ax1.text(0.1, 0.8, f'Single dome flat frame for {camera}-{detector}', fontsize=25, color='w',
ha='left', transform=ax1.transAxes)
ax2 = show_image(med_flat, ax=ax2, fig=fig)
ax2.text(0.1, 0.8, f'Master flat frame for {camera}-{detector} \n (median-combined)', fontsize=25, color='w',
ha='left', transform=ax2.transAxes)
plt.tight_layout()
med_flat_norm = med_flat / np.nanmedian(med_flat)
show_image(med_flat_norm, vmin=0.6, vmax=1.1, cmap='inferno',
show_colorbar=True, figsize=(12, 6), colorbar_label='Normalized pixel value');
Oh my gosh… The flat frames are full of dust spots and dead pixels! The dust spots are caused by dust on the optics, while the dead pixels are pixels that are not sensitive to light (you can see their sensitivity is <0.6 from the colorbar). These dust spots and dead pixels are “fixed” in the camera and will show up in every science frame. This highlights the need of flat frames to correct the science exposures!
Exercise 3#
Actually, the median flat we just made is not ready to use for flat-fielding. Do you know why? Does the flat frame also have a bias level and a dark current? How can we correct for this?
Summary#
We have learned how to read and understand bias, dark, and flat frames. We have also learned how to combine multiple calibration frames to make “master” calibration frames. We have also learned how to deal with outliers using sigma-clipping.
Exercise 4#
Does a dark frame also have a bias level? If so, how can we correct this?
Load all dark frames in
./calibration_frames/dark, subtract the master bias frame from each dark frame, and then combine them to make a master dark frame. Try different combining methods (mean, median, sigma-clipping) and compare the results.