In [ ]:
%pylab inline
#This is an ipython magic to import numpy and matplotlib (to produce graphs)
In [ ]:
import astropy.io.fits as fits #This library handles FITS files (There are other libraries to do so in python as well)
from skimage.morphology import disk
from skimage.filters import rank
from astropy.visualization import make_lupton_rgb
from reproject import reproject_interp, reproject_exact
from photutils import CircularAperture
from photutils import aperture_photometry
Let's take a look at the original images that you took
In [ ]:
science_image_path_g = 'data/seo_m66_g-band_180s_apagul_1.fits' #Type the path to your image
sci_g = fits.open(science_image_path_g)
sci_im_g = fits.open(science_image_path_g)[0].data
plt.imshow(sci_im_g,cmap='gray', vmax=1800, norm=matplotlib.colors.LogNorm())
plt.colorbar()
This image is not science-ready yet...
Dark image: If you take a shot with the shutter closed (i.e., no light/photons incoming in the camera) you still have a non-zero image.
In [ ]:
dark_image_path='data/dark.fits' #Type the path to your dark image
drk_im = fits.open(dark_image_path)[0].data
plt.imshow(drk_im,cmap='gray', vmax=2000)
plt.colorbar()
In [ ]:
bias_image_path = 'data/bias.fits' #Type the path to your bias image
bias_image = fits.open(bias_image_path)[0].data
plt.imshow(bias_image, cmap='gray')
plt.colorbar()
In [ ]:
plt.hist(drk_im.flatten());
plt.yscale('log')
plt.xlabel('Output counts')
plt.ylabel('Number of pixels')
Why is this?
Another interesting feature of CCD cameras is that the chips do not respond equally to the same light intensity. For example if you illuminate the camera with uniform light (this is called flat image).
In [ ]:
flat_image_path = 'data/FLAT_g-band_2016-10-06_bin1_id5908.fits' #Type the path to your flat image here
flat_image = fits.open(flat_image_path)[0].data
#You can try cmap='hot' or cmap='jet' to see how it changes
plt.imshow(flat_image, cmap='gray')
plt.colorbar()
In [ ]:
plt.hist(flat_image.flatten())
In [ ]:
def reduce_image(sci_im,drk_im,flat_im, bias_im, filter_dark=True):
from scipy.stats import mode
dkr_im = drk_im - bias_im
#First part: We take "zero" the image
#The next part is optional and averages the dark image in a 10 pixel radius
#to get rid of salt/pepper noise
if(filter_dark):
selem = disk(10) #We are going to perform averages in 10 pixel radius disks
selem2 = disk(4)
drk_im = rank.mean(drk_im, selem=selem) #We perform an average to remove salt-pepper noise
flat_im = rank.mean(flat_im, selem=selem2)
#Second part: Make every part have the same sensitivity
#flat_im = (flat_im - drk_im)/mode(flat_im-drk_im,axis=None)[0] #most common pixel value will equal 1
flat_im = (flat_im - drk_im)/np.median(flat_im-drk_im)
#Lower than 1 where the CCD is less sensitive and more than 1 where it's more sensitive
sci_im = (sci_im -drk_im)/flat_im
#Error image
return sci_im
Let's create a better image!
In [ ]:
new_sci_image_g = reduce_image(sci_im_g,drk_im,flat_image,bias_image, filter_dark=False)
In [ ]:
plt.imshow(new_sci_image_g, cmap='gray', vmax=4000, vmin=50, norm=matplotlib.colors.LogNorm())
plt.colorbar()
Compare to the original!
In [ ]:
fig, ax = plt.subplots(nrows=1,ncols=3,figsize=(10,8))
ax[0].imshow(sci_im_g,cmap='gray',vmax=1800, norm=matplotlib.colors.LogNorm())
ax[0].set_title('Before reduction')
ax[1].imshow(new_sci_image_g,cmap='gray',vmax=2000, vmin=50, norm=matplotlib.colors.LogNorm())
ax[1].set_title('After reduction')
ax[2].imshow(sci_im_g-new_sci_image_g,cmap='gray', vmax=1050, vmin=1000)
ax[2].set_title('Difference')
In [ ]:
science_image_path_r = 'data/seo_m66_r_180s_apagul_1.fits'
sci_im_r = fits.open(science_image_path_r)[0].data
science_image_path_i = 'data/seo_m66_i-band_180s_apagul_1.fits'
sci_im_i = fits.open(science_image_path_i)[0].data
In [ ]:
flat_r = fits.open('data/FLAT_r-band_2016-10-06_bin1_id5906.fits')[0].data
flat_i = fits.open('data/FLAT_i-band_2016-10-06_bin1_id5907.fits')[0].data
Reduce the rest of images (in principle we should take a different bias image for each filter) because the CCD has different sensitivity at different wavelengths
In [ ]:
new_sci_image_r = reduce_image(sci_im_r,drk_im,flat_r,bias_image)
new_sci_image_i = reduce_image(sci_im_i,drk_im,flat_i,bias_image)
An example from SDSS:
In [ ]:
# Read in the three images downloaded from here:
# g: http://dr13.sdss.org/sas/dr13/eboss/photoObj/frames/301/1737/5/frame-g-001737-5-0039.fits.bz2
# r: http://dr13.sdss.org/sas/dr13/eboss/photoObj/frames/301/1737/5/frame-r-001737-5-0039.fits.bz2
# i: http://dr13.sdss.org/sas/dr13/eboss/photoObj/frames/301/1737/5/frame-i-001737-5-0039.fits.bz2
g = fits.open('data/frame-g-001737-5-0039.fits.bz2')[0]
r = fits.open('data/frame-r-001737-5-0039.fits.bz2')[0]
i = fits.open('data/frame-i-001737-5-0039.fits.bz2')[0]
# remap r and i onto g
r_new, r_mask = reproject_interp(r, g.header)
i_new, i_mask = reproject_interp(i, g.header)
# zero out the unmapped values
i_new[np.logical_not(i_mask)] = 0
r_new[np.logical_not(r_mask)] = 0
# red=i, green=r, blue=g
# make a file with the default scaling
rgb_default = make_lupton_rgb(i_new, r_new, g.data, filename="ngc6976-default.jpeg")
# this scaling is very similar to the one used in Lupton et al. (2004)
rgb = make_lupton_rgb(i_new, r_new, g.data, Q=10, stretch=0.5, filename="ngc6976.jpeg")
In [ ]:
plt.imshow(rgb)
If you want to know more about Jupyter: https://github.com/fjaviersanchez/JupyterTutorial/blob/master/TutorialJupyter.ipynb
Astronomers use the magnitude scale to characterize the bright of an object. With the magnitude scale you quantify the brightness of an object by comparing it with other objects. Astronomers have agreed to use "Vega" as the zero magnitude point (like the freezing point for water is the zero-point for the Celsius temperature scale). The magnitude scale goes "backwards" in the sense that brighter objects have smaller magnitude. For example the Sun has magnitude -27, the full Moon -13, and Venus -5.
How can we measure magnitudes from an image?
A first approach is to use an object which magnitude we know, called "standard" and refer the rest of the objects in an image to it.
But what do you use to count the total brightness of an object?
In [ ]:
positions = [(550., 600.), (450., 500.)] #Change it and include the position of an object in your image
apertures = CircularAperture(positions, r=20.)
phot_table = aperture_photometry(new_sci_image_g, apertures)
In [ ]:
print phot_table
Astronomers use both disks and complicated shapes to measure the brightness and they then refer to a known object!
In [ ]: