4.4 Working with FITS Data in Python

A comprehensive demonstration of all Python tools available for HST data analysis is outside the scope of this document. Users are encouraged to consult the documentation for the specific tools they are interested in using for examples and syntax, as well as any additional information such as limitations and known bugs.

Python arrays are zero-indexed, and the convention is to specify indices of two-dimensional data as [row, column] or [y, x]. This indexing order is different than other languages and tools familiar to astronomers such as IDL and IRAF. When specifying ranges of indices, the first element is inclusive, while the second is exclusive. For example, the indices specified by [:15] include all elements 0 through 14 with the zeroth index being implicit by omission. One can also use negative indices to count backwards from the end of an array. For example, we can select the first 2,000 rows and all columns of a 2048 row x 4096 column array using either [0:2001, :] or [:-48, :]. Notice that by not specifying any indices, we return all indices.

4.4.1 Opening and Updating FITS Files

As discussed in Section 3.2, FITS files are comprised of multiple HDUs, which are a combination of an ASCII header and a data array. In Python, FITS files can be read into a HDUList object using the astropy.io.fits.open() function (see http://docs.astropy.org/en/stable/io/fits/ for more information). These HDUList objects can be indexed like arrays to access the desired HDU, or alternatively indexed using the EXTNAME and EXTVER discussed in Section 3.2.

The header and data array of an HDU can be accessed using the ".data" and ".header" attributes. For example:

from astropy.io import fits

with fits.open('jdl803dxq_flt.fits') as hdu:

	science_data1 = hdu[1].data
	science_header1 = hdu[1].header

	error_data2 = hdu['err',2].data
	error_header2 = hdu['err',2].header

One can easily update the header and data array of an HDU as follows:

from astropy.io import fits
import numpy as np

with fits.open('jdl803dxq_flt.fits') as hdu:

	hdu[0].header['FILTER1'] = 'NEW_FILTER_NAME'
	hdu[1].data = np.zeros(2048, 4096)

In this example, we have updated the keyword FILTER1 with the new value "NEW_FILTER_NAME" and we have replaced the first science extension with an 2048 row x 4096 column array of all zeros. Image data stored within FITS files is returned as a numpy ndarray object. Note that our changes will not be written out to the FITS file unless we call the astropy.io.fits.open() function with the optional argument mode='update'. Also notice that we have opened the FITS file inside of a "with block." Unindenting after this block will automatically close the file. We can repeat the above example, saving our changes to the FITS file, and explicitly close the file without using the "with block" syntax:

from astropy.io import fits 
import numpy as np 

hdu = fits.open('jdl803dxq_flt.fits', mode='update')

hdu[0].header['FILTER1'] = 'NEW_FILTER_NAME'
hdu[1].data = np.zeros(2048, 4096)


Use of the "with block" is often considered more readable and safer for file input/output, and one may encounter it in many examples.

4.4.2 FITS Table Data

Binary FITS tables are stored just like image arrays inside the HDUList object, and we can also access these with the ".data" attribute. By simply returning the ".data" attribute of the HDUList object, we obtain a FITS_rec class. We can more easily access the binary table data using the astropy.table.Table() class (see http://docs.astropy.org/en/stable/table/ for more information). We can use this to access, for example, the table information contained within an HST association file.

from astropy.io import fits
from astropy.table import Table

with fits.open('icdm0a070_asn.fits') as hdu:
	asn_table = Table(hdu[1].data)


This will yield the resulting table:

-------------- -------------- --------
ICDM0AEIQ      EXP-DTH            True
ICDM0AEKQ      EXP-DTH            True
ICDM0AEMQ      EXP-DTH            True
ICDM0A070      PROD-DTH           True

Astropy tables are indexed just like other array-like objects, but we can also use the column names to our advantage. For example, we can select the rootnames (the MEMNAME column) of the first three entries in the table (the individual exposures in the observation; the last row represents the combined product of the association) using asn_table['MEMNAME'][:-1].data. By including the ".data" attribute here, we convert the indexed result into a numpy array.