Handling FITS files¶
%matplotlib inline
import matplotlib.pyplot as plt
Documentation
For more information about the features presented below, you can read the astropy.io.fits docs.
Data
The data used in this page is a model of the gamma-ray sky background used for the LAT instrument on the Fermi telescope, as well as the Fermi/LAT point source catalog.
Reading FITS files and accessing data
Opening a FITS file is relatively straightforward. We can open the background model included in the tutorial files:
from astropy.io import fits
hdulist = fits.open('data/gll_iem_v02_P6_V11_DIFFUSE.fit')
The returned object, hdulist, (an instance of the HDUList class) behaves like a Python list, and each element
maps to a Header-Data Unit (HDU) in the FITS file. You can view more
information about the FITS file with:
hdulist.info()
As we can see, this file contains two HDUs. The first contains an image cube, the second a table with one column of double-precision floating point values. To access the primary HDU, which contains the main data, you can then do:
hdu = hdulist[0]
The hdu object then has two important attributes: data, which behaves
like a Numpy array, can be used to access the data, and header, which
behaves like a dictionary, can be used to access the header information.
First, we can take a look at the data:
hdu.data.shape
This tells us that it is a 3-d cube. We can now take a peak at the header:
hdu.header
which shows that this is a plate carrée (-CAR) projection in galactic
coordinates, and the third axis is photon energy. We can access individual
header keywords using standard item notation:
hdu.header['TELESCOP']
hdu.header['INSTRUME']
Provided that we started up ipython with the --pylab flag, we can plot
one of the slices in photon energy:
plt.imshow(hdu.data[0, ...], origin='lower')
Note that this is just a plot of an array, so the coordinates are just pixel coordinates at this stage. The data is stored with longitude increasing to the right (the opposite of the normal convention), but the Level 3 problem at the bottom of this page shows how to correctly flip the image.
Modifying data or header information in a FITS file object is easy. We can update existing header keywords:
hdu.header['TELESCOP'] = "Fermi Gamma-ray Space Telescope"
or add new ones:
hdu.header['MODIFIED'] = '2014-07-06' # adds a new keyword
and we can also change the data, for example extracting only the first slice in photon energy:
hdu.data = hdu.data[0,:,:]
Note that this does not change the original FITS file, simply the FITS file object in memory. Since the data is now 2-dimensional, we can remove the WCS keywords for the third dimension:
hdu.header.remove('CRPIX3')
hdu.header.remove('CRVAL3')
hdu.header.remove('CDELT3')
hdu.header.remove('CUNIT3')
hdu.header.remove('CTYPE3')
You can write the FITS file object to a file with:
hdu.writeto('lat_background_model_slice.fits', clobber=True)
if you want to simply write out this HDU to a file, or:
hdulist.writeto('lat_background_model_slice_allhdus.fits', clobber=True)
if you want to write out all of the original HDUs, including the modified one, to a file.
Creating a FITS file from scratch
If you want to create a FITS file from scratch, you need to start off by creating an HDU object:
hdu = fits.PrimaryHDU()
and you can then populate the data and header attributes with whatever information you like:
import numpy as np
hdu.data = np.random.random((128, 128))
Note that setting the data automatically populates the header with basic information:
hdu.header
and you should never have to set header keywords such as NAXIS, NAXIS1, and so on manually. We can then set additional header keywords:
hdu.header['telescop'] = 'Python Observatory' # Note the keyword is case-insensitive; it will be written all caps
and we can then write out the FITS file to disk:
hdu.writeto('random_array.fits', clobber=True)
Creating a multi-extension FITS file
Many observatories format their FITS files such that no data is stored in the primary HDU. Instead, multiple image extensions are used to store each of the main image data, the data quality array, the error values, etc. This is the case for all modern HST observations, for example. The primary HDU stores no data, but does use the primary header to store metadata common to the observation.
As mentioned before, the HDUList object is similar to a Python list, and can be manipulated like one to store multiple HDUs as well as reorder them in a FITS file. For example we can create a primary HDU containing only metadata:
pri_hdu = fits.PrimaryHDU()
pri_hdu.header['telescop'] = 'Python Observatory'
and then create an image extension HDU using the ImageHDU class (this is no different from PrimaryHDU except that PrimaryHDU must come first in an HDUList; in fact, in most cases Astropy will ensure that before writing out a multi-extension file the first HDU is converted to a PrimaryHDU if possible.):
img_hdu = fits.ImageHDU(data=np.random.random((128, 128)))
Finally, add both HDUs to an HDUList and write it to disk:
hdul = fits.HDUList([pri_hdu, img_hdu])
hdul.writeto('random_array2.fits', clobber=True)
We can check that the output file is in the format we expected with the fits.info() convenience function. This is a shortcut for the pattern we saw earlier of hdul = fits.open(...); hdul.info():
fits.info('random_array2.fits')
Because an HDUList is like a Python list, new HDUs can also be appended or inserted:
new_hdu = fits.ImageHDU(data=np.random.random((256, 256)))
hdul.insert(1, new_hdu)
hdul.writeto('random_array3.fits', clobber=True)
This writes a new FITS file similar to the previous one, but with the new HDU inserted:
fits.info('random_array3.fits')
Inplace updates
In all the examples we've seen thus far, modifications have been made to parts of a FITS file and saved by writing out an entirely new FITS file to a different filename. For large files this can be an expensive and time-consuming operation, especially when only making small changes.
Both the header and data of an existing FITS file can be updated inplace by opening the file with the mode='update' option:
hdu_orig = fits.open('data/gll_iem_v02_P6_V11_DIFFUSE.fit')
tile = hdu_orig[0].data[0, 52:52 + 256, 232:232 + 256]
hdul = fits.open('random_array3.fits', mode='update') # Update the file from the previous example
hdul[0].header['OBSERVER'] = 'Vera Rubin'
hdul[1].data = tile # Replace some or all of the data values
hdul.close() # Close the file to see that all changes were saved
Confirm that the file was updated:
hdul = fits.open('random_array3.fits')
hdul.info()
hdul[0].header['OBSERVER']
plt.imshow(hdul[1].data, origin='lower')
One caveat about this is that if the file needs to be resized as the result of an update (many header cards were added or removed, or the size of the data was changed) this may still result in the file being rewritten. This happens automatically in the background, but may be more time consuming than modifications that don't significantly change the size of the file.
If you're familiar with Python's with statement, this can also be used when updating FITS files. This is a good idea to ensure that the file is closed/saved once all updates are complete:
with fits.open('random_array3.fits') as hdul:
hdul[0].header['DATE-OBS'] = '2014-07-06'
Convenience functions
Like fits.info, a few other shortcut "convenience" functions are provided. For example, in cases where you just want to access the data or header in a specific HDU, you can use the following convenience functions:
data = fits.getdata('data/gll_iem_v02_P6_V11_DIFFUSE.fit')
header = fits.getheader('data/gll_iem_v02_P6_V11_DIFFUSE.fit')
To get the data or header for an HDU other than the first, you can specify the
extension name or index. The second HDU is called energies, so we can do:
data = fits.getdata('data/gll_iem_v02_P6_V11_DIFFUSE.fit', extname='energies')
or:
data = fits.getdata('data/gll_iem_v02_P6_V11_DIFFUSE.fit', ext=1)
and similarly for getheader. The documentation provides a list of other such functions.
Note: While these functions are useful and convenient for interactive sessions, directly manipulating the data structures returned by fits.open as described above is currently more performant and should be preferred for automated scripts.
Accessing Tabular Data
Tabular data behaves very similarly to image data such as that shown above, but the data array is a structured Numpy array which requires column access via the item notation:
from astropy.io import fits
hdulist = fits.open('data/gll_psc_v08.fit')
hdulist[1].name
Display information on the first five columns (or all columns without the slice notation, but there are many in this table):
hdulist[1].columns[:5]
hdulist[1].data['RAJ2000']
hdulist[1].data['DEJ2000']
FITS tables can also be accessed via Astropy's Table interface like we saw earlier:
from astropy.table import Table
t = Table.read('data/gll_psc_v08.fit')
Note that if a FITS file contains a table, the first table in the file is returned by default. The hdu= keyword argument to Table.read can be used to specify and HDU to read from.
t[t.colnames[:5]] # Display the first 5 columns
In a future version of Astropy, FITS tables will always be returned in this format. At present, however, there remain some issues with supporting all FITS tables with the Table class.
Practical Exercises
Level 1
Try and read in one of your own FITS files using astropy.io.fits, and
see if you can also plot the array values in Matplotlib. Also, examine the
header, and try and extract individual values. You can even try and modify
the data/header and write the data back out - but take care not to write
over the original file!
Level 2
Read in the LAT Point Source Catalog and make a scatter plot of the Galactic Coordinates of the sources (complete with axis labels). Bonus points if you can make the plot go between -180 and 180 instead of 0 and 360 degrees. Note that the Point Source Catalog contains the Galactic Coordinates, so no need to convert them.
Level 3
Using Matplotlib, make an all-sky plot of the LAT Background Model in the
Plate Carrée projection showing the LAT Point Source Catalog overlaid with
markers, and with the correct coordinates on the axes. You should do this
using only astropy.io.fits, Numpy, and Matplotlib (no WCS or
coordinate conversion library). Hint: the -CAR projection is such that the
x pixel position is proportional to longitude, and the y pixel position to
latitude. Bonus points for a pretty colormap.