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
/Users/erik/miniconda3/envs/astropy-tutorials/lib/python3.5/site-packages/matplotlib/__init__.py:913: UserWarning: axes.color_cycle is deprecated and replaced with axes.prop_cycle; please use the latter.
  warnings.warn(self.msg_depr % (key, alt_key))

The following line is needed to download the example FITS files used here.

In [2]:
from astropy.utils.data import download_file

Viewing and manipulating FITS images

In [3]:
from astropy.io import fits
In [4]:
image_file = download_file('http://data.astropy.org/tutorials/FITS-images/HorseHead.fits', cache=True )

Opening FITS files and loading the image data

I will open the FITS file and find out what it contains.

In [5]:
hdu_list = fits.open(image_file)
hdu_list.info()
Filename: /Users/erik/.astropy/cache/download/py3/2c9202ae878ecfcb60878ceb63837f5f
No.    Name         Type      Cards   Dimensions   Format
  0  PRIMARY     PrimaryHDU     161   (891, 893)   int16   
  1  er.mask     TableHDU        25   1600R x 4C   [F6.2, F6.2, F6.2, F6.2]   

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)
<class 'numpy.ndarray'>
(893, 891)

At this point, we can just close the FITS file. We have stored everything we wanted to a variable.

In [8]:
hdu_list.close()

SHORTCUT

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)
<class 'numpy.ndarray'>
(893, 891)

Viewing the image data and getting basic statistics

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]:
<matplotlib.colorbar.Colorbar at 0x1114a4e10>

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))
Min: 3759
Max: 22918
Mean: 9831.48167629
Stdev: 3032.3927542

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 ndarray.flatten() to return a 1-D numpy array.

In [12]:
print(type(image_data.flatten()))
<class 'numpy.ndarray'>
In [13]:
NBINS = 1000
histogram = plt.hist(image_data.flatten(), NBINS)

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.

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]:
[<matplotlib.text.Text at 0x1134c1a20>,
 <matplotlib.text.Text at 0x1134287f0>,
 <matplotlib.text.Text at 0x1139417f0>]

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.

In [16]:
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 [17]:
# 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 [18]:
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 [19]:
plt.imshow(final_image, cmap='gray', vmin=2.e3, vmax=3.e3)
plt.colorbar()
Out[19]:
<matplotlib.colorbar.Colorbar at 0x1166f0940>

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.

In [20]:
outfile = 'stacked_M13_blue.fits'

hdu = fits.PrimaryHDU(final_image)
hdu.writeto(outfile, clobber=True)
WARNING: AstropyDeprecationWarning: "clobber" was deprecated in version 1.3 and will be removed in a future version. Use argument "overwrite" instead. [astropy.utils.decorators]

Exercises

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 [21]:
 

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 [21]:
 

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 [21]: