Lecturer: Anupama Gadiyara
Jupyter Notebook Author: Kishalay De & Cameron Hummels
This is a Jupyter notebook lesson taken from the GROWTH Winter School 2018. For other lessons and their accompanying lectures, please see: http://growth.caltech.edu/growth-astro-school-2018-resources.html
Process raw images from a visible wavelength telescope and make them ready for photometric analysis
See GROWTH school webpage for detailed instructions on how to install these modules and packages. Nominally, you should be able to install the python modules with
pip install <module>. The external astromatic packages are easiest installed using package managers (e.g.,
#import some useful modules from astropy.io import fits import matplotlib.pyplot as plt import glob import os from astropy.stats import sigma_clipped_stats import numpy as np import subprocess import astroscrappy import pyregion import shutil
In order for this jupyter notebook to function correctly, we must have some external software installed, as described above. The following step assures that these are installed properly before getting to the rest of the content of this lesson.
dependencies = ['SWarp', 'sextractor'] def test_dependency(dep): try: subprocess.Popen(dep, stderr=subprocess.PIPE).stderr.read() print("%s is installed properly. OK" % dep) return 1 except ImportError: print("===%s IS NOT YET INSTALLED PROPERLY===" % dep) return 0 i = 0 for dep in dependencies: i += test_dependency(dep) print("\n%i out of %i dependencies installed properly." % (i, len(dependencies))) if i != len(dependencies): print("Please correctly install these programs before continuing.") else: print("You are ready to continue.")
SWarp is installed properly. OK sextractor is installed properly. OK 2 out of 2 dependencies installed properly. You are ready to continue.
The various directories where data live are defined here for use later in this notebook
curpath = os.path.abspath('.') # top level directory dataFolder = os.path.join(curpath, 'data') # data dir biasFolder = os.path.join(dataFolder, 'bias') # biases flatFolder = os.path.join(dataFolder, 'flats') # flats sciFolder = os.path.join(dataFolder, 'science') # science data procFolder = os.path.join(curpath, 'processing')# processing directory if not os.path.isdir(procFolder): os.mkdir(procFolder) else: for f in os.listdir(procFolder): os.remove(os.path.join(procFolder,f)) #clear the processing folder from previous iterations
os.chdir(sciFolder) fileList = glob.glob('*.fits') os.chdir(curpath) procList = [os.path.join(procFolder, file) for file in fileList] sciList = [os.path.join(sciFolder, file) for file in fileList]
The image data received from the telescope is in the form of a FITS file containing a 2D array of counts that represent the brightness distribution on the sky as seen by the detector. We have to remove the effects of the telescope optical chain from this image to get the true brightness distribution on the sky. There are are a few different steps involved here. We will call this image distribution $I(x,y)$ as a function of x, y coordinates.
We'll start this module by creating a master bias frame to be subtracted from each of the science and flat images. The idea is to use a median combination of multiple bias frames to get rid of transient features like cosmic rays in each of the individual bias frames, and create a 'Master bias' frame.
biasList = glob.glob(os.path.join(biasFolder,'*fits')) numBiasFiles = len(biasList) print('Found %d bias files'%numBiasFiles) #Let's load all the files into a numpy array #This is a 3D array with each element representing a 4096 x 4108 pixel array corresponding to each input image biasImages = np.zeros((4108, 4096, numBiasFiles)) for i in range(numBiasFiles): biasImages[:,:,i] = fits.open(biasList[i]).data #Let's plot up a single bias frame plt.figure(figsize=(8,8)) mean, median, std = sigma_clipped_stats(biasImages[:,:,0]) plt.figure(figsize=(8,8)) #set the scale of the image based on its statistics plt.imshow(biasImages[:,:,0], vmin = median - 2*std, vmax = median + 2*std) plt.colorbar() plt.show() #Now creating the master bias frame by doing a median combination for each pixel in the list of images masterBias = np.median(biasImages, axis=2)
Found 3 bias files
<Figure size 576x576 with 0 Axes>
Now plot the master bias frame and check what it looks like
#Lets plot the bias image mean, median, std = sigma_clipped_stats(masterBias) plt.figure(figsize=(8,8)) #set the scale of the image based on its statistics plt.imshow(masterBias, vmin = median - 2*std, vmax = median + 2*std) plt.colorbar() plt.show()