In [1]:
import numpy as np
# Set up matplotlib and use a nicer set of plot parameters
%config InlineBackend.rc = {}
import matplotlib
matplotlib.rc_file("../templates/matplotlibrc")
import matplotlib.pyplot as plt
%matplotlib inline
The following line is needed to download the example FITS files used here.
In [2]:
from astropy.utils.data import download_file
In [3]:
from astropy.io import fits
In [4]:
image_file = download_file('http://data.astropy.org/tutorials/FITS-images/HorseHead.fits', cache=True )
I will open the FITS file and find out what it contains.
In [5]:
hdu_list = fits.open(image_file)
hdu_list.info()
Generally the image information is located in the PRIMARY block. The blocks are numbered and can be accessed by indexing hdu_list.
In [6]:
image_data = hdu_list[0].data
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.
In [7]:
print(type(image_data))
print(image_data.shape)
At this point, we can just close the FITS file. We have stored everything we wanted to a variable.
In [8]:
hdu_list.close()
If you don't need to examine the FITS header, you can call fits.getdata to bypass the previous steps.
In [9]:
image_data = fits.getdata(image_file)
print(type(image_data))
print(image_data.shape)
In [10]:
plt.imshow(image_data, cmap='gray')
plt.colorbar()
# To see more color maps
# http://wiki.scipy.org/Cookbook/Matplotlib/Show_colormaps
Out[10]:
Let's get some basic statistics about our image
In [11]:
print('Min:', np.min(image_data))
print('Max:', np.max(image_data))
print('Mean:', np.mean(image_data))
print('Stdev:', np.std(image_data))
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 ndarray.flatten() to return a 1-D numpy array.
In [12]:
print(type(image_data.flatten()))
In [13]:
NBINS = 1000
histogram = plt.hist(image_data.flatten(), NBINS)
Want to use a logarithmic color scale? To do so we need to load the LogNorm object from matplotlib.
In [14]:
from matplotlib.colors import LogNorm
In [15]:
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'])
Out[15]:
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.
In [ ]:
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 ]
Now I'll stack the images by summing my concatenated list.
In [ ]:
# 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)
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.
In [ ]:
image_hist = plt.hist(final_image.flatten(), 1000)
I'll use the keywords vmin and vmax to set limits on the color scaling for imshow.
In [ ]:
plt.imshow(final_image, cmap='gray', vmin=2.e3, vmax=3.e3)
plt.colorbar()
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.
In [ ]:
outfile = 'stacked_M13_blue.fits'
hdu = fits.PrimaryHDU(final_image)
hdu.writeto(outfile, clobber=True)
Determine the mean, median, and standard deviation of a part of the stacked M13 image where there is not light from M13. Use those statistics with a sum over the part of the image that includes M13 to estimate the total light in this image from M13.
In [ ]:
Show the image of the Horsehead Nebula, but in to units of surface brightness (magnitudes per square arcsecond). (Hint: the physical size of the image is 15x15 arcminutes.)
In [ ]:
Now write out the image you just created, preserving the header the original image had, but add a keyword 'UNITS' with the value 'mag per sq arcsec'. (Hint: you may need to read the astropy.io.fits documentation if you're not sure how to include both the header and the data)
In [ ]: