<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 FITS file (Flexible Image Transport System) - Standard Astro Data Format
FITS (Flexible Image Transport System) is the data format most widely used within astronomy for transporting, analyzing, and archiving scientific data files. FITS is much more than just another image format (such as JPG or GIF) and is primarily designed to store scientific data sets consisting of multidimensional arrays (images) and 2-dimensional tables organized into rows and columns of information.
Step1:  FITS files consist of at least two parts - Header and Data
A FITS file is comprised of segments called Header/Data Units (HDUs), where the first HDU is called the Primary HDU.
The primary data array can contain a 1-D spectrum, a 2-D image, or a 3-D data cube.
Every HDU consists of an ASCII formatted Header Unit. 
Each header unit contains a sequence of fixed-length 80-character keyword (Cards)
Step2:  FITS format preserves the full dynamic range of data
Step3:  You can use masks on images
Step4:  You can use specific bins for histograms
Step5:  You can add and subtract images
Step6:  The two images above may look the same but they are not!
Subtracting the two images reveals the truth.
Step7:  FITS Tables - An astronomical example
Stellar spectra data from the ESO Library of Stellar Spectra
Step8:  Stellar spectral classes
Step9:  FITS Images - An astronomical example
Step10:  Notice the origin is in the upper left corner (the image is upside down)
Step11:  Better, the origin is in the lower left corner
World Coordinate System wcs
Step12:  RGB - Pseudocolor Images
Step13:  Sigma Clipping
Step14:  Rejected Data
Step15:  Accepted Data
Step16:  Sigma clipping an image
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import astropy.io.fits as fits
plt.style.use('ggplot')
plt.rc('axes', grid=False)   # turn off the background grid for images
Explanation: FITS file (Flexible Image Transport System) - Standard Astro Data Format
FITS (Flexible Image Transport System) is the data format most widely used within astronomy for transporting, analyzing, and archiving scientific data files. FITS is much more than just another image format (such as JPG or GIF) and is primarily designed to store scientific data sets consisting of multidimensional arrays (images) and 2-dimensional tables organized into rows and columns of information.
End of explanation
data_file = "./MyData/bsg01.fits"
my_fits_file = fits.open(data_file)
my_fits_file.info()
image_data = my_fits_file[0].data
image_header = my_fits_file[0].header
image_header
print("The image has a shape [height,width] of {0}".format(image_data.shape))
print("The image is made up of data of type {0}".format(image_data.dtype))
Explanation: FITS files consist of at least two parts - Header and Data
A FITS file is comprised of segments called Header/Data Units (HDUs), where the first HDU is called the Primary HDU.
The primary data array can contain a 1-D spectrum, a 2-D image, or a 3-D data cube.
Every HDU consists of an ASCII formatted Header Unit. 
Each header unit contains a sequence of fixed-length 80-character keyword (Cards)
End of explanation
print("The image has a maximum value of {0}".format(image_data.max()))
print("The image has a minimum value of {0}".format(image_data.min()))
fig, ax = plt.subplots(1,2)
fig.set_size_inches(12,6)
fig.tight_layout()
ax[0].imshow(image_data,cmap=plt.cm.gray)
ax[1].hist(image_data.flatten(),bins=20);
Explanation: FITS format preserves the full dynamic range of data
End of explanation
copy_data = np.copy(image_data)        # make a copy of the data to work with
cut_off = 40
mask = np.where(copy_data > cut_off)
copy_data[mask] = 60                   # You can not just throw data away, you have to set it to something.
print("The cropped image has a maximum value of {0}".format(copy_data.max()))
print("The cropped image has a minimum value of {0}".format(copy_data.min()))
Explanation: You can use masks on images
End of explanation
my_bins = np.arange(0,110,5)
my_bins
fig, ax = plt.subplots(1,2)
fig.set_size_inches(12,6)
fig.tight_layout()
ax[0].imshow(copy_data,cmap=plt.cm.gray)
ax[1].hist(image_data.flatten(),bins=my_bins)
ax[1].hist(copy_data.flatten(),bins=my_bins, alpha = 0.40);
Explanation: You can use specific bins for histograms
End of explanation
another_image_file = "./MyData/bsg02.fits"
another_image_data = fits.getdata(another_image_file)     # a quick way to just get the data
fig, ax = plt.subplots(1,2)
fig.set_size_inches(12,6)
fig.tight_layout()
ax[0].set_title("Original Image")
ax[1].set_title("Another Image")
ax[0].imshow(image_data, cmap=plt.cm.gray)
ax[1].imshow(another_image_data, cmap=plt.cm.gray);
Explanation: You can add and subtract images
End of explanation
real_image = image_data - another_image_data                 # Subtract the images pixel by pixel
fig, ax = plt.subplots(1,3)
fig.set_size_inches(12,6)
fig.tight_layout()
ax[0].set_title("Original Image")
ax[1].set_title("Another Image")
ax[2].set_title("Real Image!")
ax[0].imshow(image_data, cmap=plt.cm.gray)
ax[1].imshow(another_image_data, cmap=plt.cm.gray);
ax[2].imshow(real_image, cmap=plt.cm.gray);
print("The real image has a maximum value of {0}".format(real_image.max()))
print("The real image has a minimum value of {0}".format(real_image.min()))
Explanation: The two images above may look the same but they are not!
Subtracting the two images reveals the truth.
End of explanation
my_spectra_file = './MyData/Star_G5.fits'
my_spectra_fits = fits.open(my_spectra_file)
my_spectra_fits.info()
spectra_data = my_spectra_fits[0].data
spectra_header = my_spectra_fits[0].header
spectra_header[0:15]
# The FITS header has the information to make an array of wavelengths
start_wavelength = spectra_header['CRVAL1']
num_of_points = spectra_header['NAXIS1']
width_of_points  = spectra_header['CDELT1']
end_wavelength = start_wavelength + (num_of_points * width_of_points)
wavelength_range = np.arange(start_wavelength,end_wavelength,width_of_points)
fig, ax = plt.subplots(2,1)
fig.set_size_inches(11,8.5)
fig.tight_layout()
# Full spectra
ax[0].plot(wavelength_range, spectra_data, color='b')
ax[0].set_ylabel("Flux")
ax[0].set_xlabel("Wavelength [angstroms]")
# Just the visible range with the hydrogen Balmer lines
ax[1].set_xlim(4000,7000)
ax[1].set_ylim(0.6,1.2)
ax[1].plot(wavelength_range, spectra_data, color='b')
ax[1].set_ylabel("Flux")
ax[1].set_xlabel("Wavelength [angstroms]")
H_Balmer = [6563,4861,4341,4102,3970,3889,3835,3646]
ax[1].vlines(H_Balmer,0,2, color='r', linewidth=3, alpha = 0.25);
Explanation: FITS Tables - An astronomical example
Stellar spectra data from the ESO Library of Stellar Spectra
End of explanation
import glob
star_list = glob.glob('./MyData/Star_*5.fits')
star_list
fig, ax = plt.subplots(1,1)
fig.set_size_inches(9,5)
fig.tight_layout()
# Full spectra
ax.set_ylabel("Flux")
ax.set_xlabel("Wavelength [angstroms]")
ax.set_ylim(0.0, 1.05)
for file in star_list:
        
    spectra = fits.getdata(file)
    
    spectra_normalized = spectra / spectra.max()
    
    ax.plot(wavelength_range, spectra_normalized, label=file)
ax.legend(loc=0,shadow=True);   
Explanation: Stellar spectral classes
End of explanation
star_file = "./MyData/star_field.fits"
star_fits_file = fits.open(star_file)
star_fits_file.info()
star_data = star_fits_file[0].data
star_header = star_fits_file[0].header
plt.imshow(star_data, cmap=plt.cm.gray)
Explanation: FITS Images - An astronomical example
End of explanation
plt.imshow(star_data, origin='lower', cmap=plt.cm.gray)
Explanation: Notice the origin is in the upper left corner (the image is upside down)
End of explanation
from astropy.wcs import WCS
my_wcs = WCS(star_header)
my_wcs
fig = plt.figure()
ax = fig.add_subplot(111, projection=my_wcs)
fig.set_size_inches(6,6)
fig.tight_layout()
ax.grid(color='y', ls='-')
ax.set_xlabel('Right Ascension (J2000)')
ax.set_ylabel('Declination (J2000)')
ax.imshow(star_data, origin='lower', cmap=plt.cm.gray);
Explanation: Better, the origin is in the lower left corner
World Coordinate System wcs
End of explanation
from astropy.visualization import make_lupton_rgb
red_img = fits.getdata("./MyData/m51_IR.fits").astype(float)
green_img = fits.getdata("./MyData/m51_red.fits").astype(float)
blue_img= fits.getdata("./MyData/m51_blue.fits").astype(float)
clean_r =  red_img - np.median(red_img)
clean_g =  green_img - np.median(green_img)
clean_b =  blue_img - np.median(blue_img)
fig, ax = plt.subplots(1,3)
fig.set_size_inches(12,6)
fig.tight_layout()
ax[0].set_title("R")
ax[1].set_title("G")
ax[2].set_title("B")
ax[0].imshow(clean_r, origin='lower', cmap=plt.cm.gray)
ax[1].imshow(clean_g, origin='lower', cmap=plt.cm.gray);
ax[2].imshow(clean_b, origin='lower', cmap=plt.cm.gray);
image = make_lupton_rgb(clean_r, clean_g, clean_b, stretch = 5000)
fig, ax = plt.subplots(1,1)
fig.set_size_inches(6,6)
fig.tight_layout()
ax.imshow(image, origin='lower', cmap=plt.cm.gray);
Explanation: RGB - Pseudocolor Images
End of explanation
from astropy import stats
my_data = np.array([1, 5, 6, 8, 100, 5, 3, 2, 4, 5])
my_data.mean(), my_data.std()
filtered_data = stats.sigma_clip(my_data, sigma=2, iters=5) 
filtered_data
Explanation: Sigma Clipping
End of explanation
my_data[filtered_data.mask]
Explanation: Rejected Data
End of explanation
my_data[~filtered_data.mask]
Explanation: Accepted Data
End of explanation
star_data.mean(), star_data.std(), star_data.max(), star_data.min()
clip_star = stats.sigma_clip(star_data, sigma=8, iters=5) 
copy_data = np.copy(star_data)        # make a copy of the data to work with
copy_data[~clip_star.mask] = star_data.min()
fig, ax = plt.subplots(1,2)
fig.set_size_inches(12,6)
fig.tight_layout()
ax[0].imshow(star_data, origin='lower', cmap=plt.cm.gray)
ax[1].imshow(copy_data, origin='lower', cmap=plt.cm.gray);
Explanation: Sigma clipping an image
End of explanation
<END_TASK>