<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 The following line is needed to download the example FITS files used here.
Step1:  Viewing and manipulating FITS images
Step2:  Opening FITS files and loading the image data
I will open the FITS file and find out what it contains.
Step3:  Generally the image information is located in the <code>PRIMARY</code> block. The blocks are numbered and can be accessed by indexing <code>hdu_list</code>.
Step4:  You data is now stored as a 2-D numpy array.  Want to know the dimensions of the image?  Just look at the shape of the array.
Step5:  At this point, we can just close the FITS file.  We have stored everything we wanted to a variable.
Step6:  SHORTCUT
If you don't need to examine the FITS header, you can call fits.getdata to bypass the previous steps.
Step7:  Viewing the image data and getting basic statistics
Step8:  Let's get some basic statistics about our image
Step9:  Plotting a histogram
To make a histogram with matplotlib.pyplot.hist(), I need to cast the data from a 2-D to array to something one dimensional.
In this case, I am using the iterable python object img_data.flat.
Step10:  Displaying the image with a logarithmic scale
Want to use a logarithmic color scale? To do so we need to load the LogNorm object from matplotlib.
Step11:  Basic image math
Step12:  Now I'll stack the images by summing my concatenated list.
Step13:  I'm going to show the image, but I want to decide on the best stretch. To do so I'll plot a histogram of the data.
Step14:  I'll use the keywords vmin and vmax to set limits on the color scaling for imshow.
Step15:  Writing image data to a FITS file
This is easy to do with the writeto() method.
You will receive an error if the file you are trying to write already exists.  That's why I've set clobber=True.
<ASSISTANT_TASK:>
Python Code:
from astropy.utils.data import download_file
Explanation: The following line is needed to download the example FITS files used here.
End of explanation
from astropy.io import fits
image_file = download_file('http://data.astropy.org/tutorials/FITS-images/HorseHead.fits', cache=True )
Explanation: Viewing and manipulating FITS images
End of explanation
hdu_list = fits.open(image_file)
hdu_list.info()
Explanation: Opening FITS files and loading the image data
I will open the FITS file and find out what it contains.
End of explanation
image_data = hdu_list[0].data
Explanation: Generally the image information is located in the <code>PRIMARY</code> block. The blocks are numbered and can be accessed by indexing <code>hdu_list</code>.
End of explanation
print(type(image_data))
print(image_data.shape)
Explanation: You data is now stored as a 2-D numpy array.  Want to know the dimensions of the image?  Just look at the shape of the array.
End of explanation
hdu_list.close()
Explanation: At this point, we can just close the FITS file.  We have stored everything we wanted to a variable.
End of explanation
image_data = fits.getdata(image_file)
print(type(image_data))
print(image_data.shape)
Explanation: SHORTCUT
If you don't need to examine the FITS header, you can call fits.getdata to bypass the previous steps.
End of explanation
plt.imshow(image_data, cmap='gray')
plt.colorbar()
# To see more color maps
# http://wiki.scipy.org/Cookbook/Matplotlib/Show_colormaps
Explanation: Viewing the image data and getting basic statistics
End of explanation
print('Min:', np.min(image_data))
print('Max:', np.max(image_data))
print('Mean:', np.mean(image_data))
print('Stdev:', np.std(image_data))
Explanation: Let's get some basic statistics about our image
End of explanation
print(type(image_data.flat))
NBINS = 1000
histogram = plt.hist(image_data.flat, NBINS)
Explanation: Plotting a histogram
To make a histogram with matplotlib.pyplot.hist(), I need to cast the data from a 2-D to array to something one dimensional.
In this case, I am using the iterable python object img_data.flat.
End of explanation
from matplotlib.colors import LogNorm
plt.imshow(image_data, cmap='gray', norm=LogNorm())
# I chose the tick marks based on the histogram above
cbar = plt.colorbar(ticks=[5.e3,1.e4,2.e4])
cbar.ax.set_yticklabels(['5,000','10,000','20,000'])
Explanation: Displaying the image with a logarithmic scale
Want to use a logarithmic color scale? To do so we need to load the LogNorm object from matplotlib.
End of explanation
image_list = [ download_file('http://data.astropy.org/tutorials/FITS-images/M13_blue_000'+n+'.fits', cache=True ) \
              for n in ['1','2','3','4','5'] ]
# The long way
image_concat = []
for image in image_list:
    image_concat.append(fits.getdata(image))
    
# The short way
#image_concat = [ fits.getdata(image) for image in IMAGE_LIST ]
Explanation: Basic image math: image stacking
You can perform math with the image data like any other numpy array.  In this particular example, I will stack several images of M13 taken with a ~10'' telescope.
I open a series of FITS files and store the data in a list, which I've named image_concat.
End of explanation
# The long way
final_image = np.zeros(shape=image_concat[0].shape)
for image in image_concat:
    final_image += image
# The short way
#final_image = np.sum(image_concat, axis=0)
Explanation: Now I'll stack the images by summing my concatenated list.
End of explanation
image_hist = plt.hist(final_image.flat, 1000)
Explanation: I'm going to show the image, but I want to decide on the best stretch. To do so I'll plot a histogram of the data.
End of explanation
plt.imshow(final_image, cmap='gray', vmin=2.e3, vmax=3.e3)
plt.colorbar()
Explanation: I'll use the keywords vmin and vmax to set limits on the color scaling for imshow.
End of explanation
outfile = 'stacked_M13_blue.fits'
hdu = fits.PrimaryHDU(final_image)
hdu.writeto(outfile, clobber=True)
Explanation: Writing image data to a FITS file
This is easy to do with the writeto() method.
You will receive an error if the file you are trying to write already exists.  That's why I've set clobber=True.
End of explanation
<END_TASK>