Lab 4 - Using Python to Read and Display .fits Files


In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
%matplotlib inline

The cell above imports all of the .fits file handling functions from the astropy python library. To call any function from the library, type fits.functionname. To see a list of available functions, click in the cell below, move your cursor to just after the period, and hit tab. You'll see a dropdown list of available functions.


In [ ]:
## don't execute this cell. It's incomplete and won't do anything.
fits.

To start, let's read in a .fits file with fits.open


In [ ]:
data = fits.open('M13-001_I.proc.fit')

We've created a python object called data above, but data is a special type of python object called an HDUlist, which is more or less just a list of python objects. The .fits file format can have many extensions besides just a data array and a header (other data arrays, other text files, etc.). In most cases in astronomy, and in probably all the cases for this class, we are only concerned with the first object in the list, which contains the main data array and header. This first object has index 0. One useful method for an HDUlist is .info, which will tell you about how many objects are in the list. If you execute the cell below, you can see that there's only one, called PRIMARY, and .info also lists its dimensions and the type of numbers stored in the array (64-bit floats)


In [ ]:
data.info()

More useful may be displaying the header in python, just like we did in DS9. The header is associated with the Primary HDU, so you need to index data with [0] in order to see it.


In [ ]:
data[0].header

To display a .fits image, we can use the function imshow, as in the cell below, but first we need to make an ordinary 2D array from the image data, which is what the first line is for.


In [ ]:
image = data[0].data
plt.imshow(image)

You can probably tell that the image is of a globular cluster, but obviously the choice of scale, etc. is not ideal. If you want to know what the colorscale and its min and max look like, add the line plt.colorbar(), as below. At the same time, let's increase the size of the image so that we can see it better.


In [ ]:
plt.figure(figsize=(15,7.5))
plt.imshow(image)
plt.colorbar()

Exercise 1

Spend no more than 5 minutes playing with the imshow function's keywords cmap (for "color map", options below) and vmin and vmax, the minimum and maximum values for the colorbar. Stop when you think you can better see the stars in the image.


In [ ]:
#plotting code goes here

You can also zoom in on a region of the image by using the indices of the pixels you want. Note though, that the indices are in the form [ymin:ymax,xmin:xmax] AND that, if you look at the image above, python images display with the pixel 0,0 in the upper left corner rather than the lower left, as is more typical.

Since most astronomical images are oriented for a lower left pixel origin (so that, for example, North is up and East is left in the image), you should generally use the option origin="lower" when displaying. For example, to zoom in on the center of the cluster:


In [ ]:
plt.figure(figsize=(8,8))
plt.imshow(image[450:600,700:850], origin="lower")
plt.colorbar()

Exercise 2

Create a nicely scaled image of the globular cluster that shows most of its stars and not a lot of empty space.


In [ ]:
#code for image here