<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Let's take a look at the original images that you took
Step1:  This image is not science-ready yet...
Dark image
Step2:  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).
Step3:  Let's create a better image!
Step4:  Compare to the original!
Step5:  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
Step6:  An example from SDSS
Step7:  If you want to know more about Jupyter
<ASSISTANT_TASK:>
Python Code:
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()
Explanation: Let's take a look at the original images that you took
End of explanation
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()
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()
plt.hist(drk_im.flatten());
plt.yscale('log')
plt.xlabel('Output counts')
plt.ylabel('Number of pixels')
Explanation: 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.
End of explanation
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()
plt.hist(flat_image.flatten())
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
Explanation: 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).
End of explanation
new_sci_image_g = reduce_image(sci_im_g,drk_im,flat_image,bias_image, filter_dark=False)
plt.imshow(new_sci_image_g, cmap='gray', vmax=4000, vmin=50, norm=matplotlib.colors.LogNorm())
plt.colorbar()
Explanation: Let's create a better image!
End of explanation
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')
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
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
Explanation: Compare to the original!
End of explanation
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)
Explanation: 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
End of explanation
# 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")
plt.imshow(rgb)
Explanation: An example from SDSS:
End of explanation
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)
print phot_table
Explanation: If you want to know more about Jupyter:
https://github.com/fjaviersanchez/JupyterTutorial/blob/master/TutorialJupyter.ipynb
Aperture photometry
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?
Use the brightest pixel?
Add the brightness in a certain radius?
Count only the pixels which belong to each object?
End of explanation
<END_TASK>