Calibration frames#
Prerequisites
Basic statistics (covered in this notebook)
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
I suppose Jenny already introduced bias, dark, flat frames. Here I will demonstrate what real bias frames look like, and why we need to stack them.
%load_ext autoreload
%autoreload 2
import os, sys
import numpy as np
import scipy
import matplotlib
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
try:
import google.colab
IN_COLAB = True
except:
IN_COLAB = False
if not IN_COLAB:
if not os.path.exists("../_static/ObsAstroData/calibration_frames/bias/"):
os.makedirs("../_static/ObsAstroData/calibration_frames/bias/")
os.chdir('../_static/ObsAstroData/calibration_frames/')
if os.path.isfile('b0860c1.fits'):
print('Data is already there')
else:
print('Download data')
# !wget https://tigress-web.princeton.edu/~jiaxuanl/cutoutRGB_NIRCAM.pkl
Download data
1. Bias frames#
As introduced in class, a voltage is applied to every pixel to make sure the pixel value is always positive. This voltage introduces a “bias” to pixel values. However, this bias is slightly different for each pixel and might change over time. Therefore, we need to take some zero-time exposures (such that dark current and real photons don’t 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.
Read a bias frame#
from astropy.io import fits
It’s our first time reading a FITS file. 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('./bias/b0860c1.fits')
[<astropy.io.fits.hdu.image.PrimaryHDU object at 0x7f91c2a327c0>]
fits.open('./bias/b0860c1.fits').info()
Filename: ./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('./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)
np.median(img), np.std(img)
(817.0, 10.545205968480271)
The pixel values in this bias frame are centered around ~817, and have a dispersion of ~10.5. This dispersion is largely contributed by the read-out noise, a noise that added to each pixel when reading the charge out. Read-out noise may be viewed as a “toll” that must be paid for reading the stored charge.
np.sum(img < 0)
0
Thanks to the bias voltage, none of the pixels are negative. We can further visualize the bias frame as follows.
from utils import show_image
show_image(img, show_colorbar=True);
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.
fig, [ax1, ax2] = plt.subplots(1, 2, figsize=(12, 4))
ax1.hist(img.flatten(), bins=30)
ax1.set_xlabel('Pixel value [counts]')
ax1.set_ylabel('Frequency')
ax2.hist(img.flatten(), bins=30)
ax2.set_yscale('log')
ax2.set_xlabel('Pixel value [counts]')
ax2.set_ylabel('log Frequency')
plt.tight_layout()
During the Magellan observation, we took many bias frames (why??). So let’s load another bias frame. What do you expect? Will it be very different from the one above?
show_image(fits.open('./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.
Combine bias frames to make “master” bias#
We have already learned from this notebook that combining images will reduce the noise level. Belowe we try to combine 40 bias frames.
biases = [fits.open(f'./bias/b0{ind}c1.fits')[0].data for ind in range(860, 900)]
biases = np.array(biases)
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));
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 red in the figure below (they are tiny dots).
hdu = fits.open('./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='Reds', vmin=0., vmax=0.2, 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'./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('./bias/master_bias.fits', overwrite=True)
2. Dark frame#
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'./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.