4.6. About astropy.fits#
4.6.1. Basic Information about FITS#
FITS: For detailed information about FITS (Flexible Image Transport System), you can refer to the official website and the latest FITS standard document.
§4.4.2 of the document mentioned above describes some frequently encountered standard keywords in FITS files, such as
DATE-OBS,OBJECT,TELESCOP,INSTRUME,OBSERVER,BSCALE,BZERO,BUNIT, etc.Appendix C provides a summary of the standard keywords.
Chap. 8 provides a brief explanation of the World Coordinate System (WCS).
HDU: Header Data Unit. It represents a component of a FITS file that consists of a header and possibly associated data.
HDU List: In Python, an HDU List is a list-like object that contains multiple HDU objects. It allows storing and manipulating multiple HDUs within a single FITS file.
MEF: Multi-Extension FITS, where multiple HDUs are present in a single FITS file. Each HDU can have its own header and data. The simplest FITS file contains only the primary HDU (so it is not an MEF).
In some cases, FITS files consist of multiple HDUs, allowing for the storage of various types of data and associated metadata within a single file. The use of MEF enables the organization and management of complex data structures within FITS files.
4.6.2. Prerequisites#
Before starting this notebook, it is essential to go through the following resources:
MUST READ:
Viewing and manipulating FITS images: This tutorial provides an overview of working with FITS images using the example of the Horsehead nebula image.
Same as above: This tutorial covers similar concepts but may have slightly different implementations, and the outputs shown might be outdated.
Highly recommended to read:
Learn Astropy: It is highly recommended to explore the Learn Astropy website, as it offers a comprehensive set of tutorials and resources covering various aspects of the Astropy library.
Other References
For additional resources or advanced tutorials for further research purposes, you can visit the Learn Astropy website and select the “fits” section from the left sidebar.
from pathlib import Path
import numpy as np
import time
from astropy.io import fits
from astropy.nddata import CCDData
4.6.3. FITS from Astropy#
Making FITS from ndarray#
If you only have ndarray (e.g., numpy array), how can we convert it to FITS?
Making and Saving with CCDData#
The simplest way is to use CCDData:
test_data_01 = np.ones((100, 100))
test_ccd_01 = CCDData(data=test_data_01, header=None, unit='adu')
print(type(test_ccd_01))
print(test_ccd_01.data)
<class 'astropy.nddata.ccddata.CCDData'>
[[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]
...
[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]]
test_ccd_01.header
OrderedDict()
Note that
unitis a mandatory forCCDData(which is, in my humble opinion, very annoying…)The
test_ccd_01.headerdoes not contain anything and returned asOrderedDict().
It is always recommended to attach header to each FITS file, but let me just save it for simple example:
test_ccd_01.write("test_01.fits", overwrite=True)
Indeed, there is a dedicated writer called fits_ccddata_writer in Astropy, which is specifically designed for saving a CCDData object to a FITS file. This writer provides various options for saving the data; refer to the documentation for details of the available saving options and their usage.
Creating and Saving FITS Data with astropy.fits#
When dealing with FITS, astropy.fits is a more general approach compared to astropy.nddata.CCDData class. FITS files can contain various types of data, such as CCD images, star catalogs, simulation data, etc. CCDData class can only deal with CCD-like (intended to be 2-D) image files.
There are at least two ways to save 2D CCD-like data using astropy.fits:
Using
PrimaryHDU(): A simple choice for a single extension.Using
ImageHDU(): A suitable choice for cases involving Multi-Extension FITS (MEF).
In this note, I will demonstrate how to create an HDUList and save it using the straightforward PrimaryHDU() approach. For MEF cases, you can find an example later in this note. You can also refer to the Astropy documentation on HDU to explore other classes of HDUs and their usage.
test_data_01 = np.ones((100, 100))
test_hdu_01 = fits.PrimaryHDU(data=test_data_01, header=None)
print(type(test_hdu_01))
print(test_hdu_01.data)
<class 'astropy.io.fits.hdu.image.PrimaryHDU'>
[[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]
...
[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]]
# If you want, you can make HDUList.
# But it is unnecessary because we have only one extension.
# hdul = fits.HDUList([test_hdu_01])
test_hdu_01.writeto("test_01.fits", overwrite=True, output_verify='fix')
For the options when writing, please look at the official document (here).
Reading FITS#
Reading with CCDData.read()#
Let’s read this:
test_ccd_01_read = CCDData.read("test_01.fits", unit="adu")
print(type(test_ccd_01_read))
print(test_ccd_01_read.data)
<class 'astropy.nddata.ccddata.CCDData'>
[[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]
...
[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]]
test_ccd_01_read.header
SIMPLE = T / conforms to FITS standard
BITPIX = -64 / array data type
NAXIS = 2 / number of array dimensions
NAXIS1 = 100
NAXIS2 = 100
EXTEND = T
Please note that when saving data using astropy.fits, the header is automatically generated for you.
The automatically generated header includes the following mandatory keywords, which follow the FITS standard:
SIMPLE = T: Indicates that the file is in standard FITS format.BITPIX = -64: Indicates that the data is stored as a 64-bit floating-point format.NAXIS = 2: Specifies that the data has 2 axes.NAXIS1 = 100andNAXIS2 = 100: Specify the length of the first axis (X-axis) and the second axis (Y-axis) of the data, respectively. Note that the X-axis corresponds tonumpyaxis 1, and the Y-axis corresponds tonumpyaxis 0.BUNIT = 'adu ': Specifies that the unit of the pixel values is in ADU (analog-to-digital units).
Based on the Astropy documentation, the relationship between the BITPIX value and the corresponding numpy data type is as follows:
BITPIX Numpy Data Type
8 numpy.uint8 (unsigned integer)
16 numpy.int16
32 numpy.int32
-32 numpy.float32
-64 numpy.float64
Reading with fits.open()#
Once it is saved, you can also read the FITS using astropy.fits:
test_ccd_01_hdul = fits.open("test_01.fits")
print(type(test_ccd_01_hdul))
<class 'astropy.io.fits.hdu.hdulist.HDUList'>
Now it’s not CCDData, but HDUList object. It’s like a python list, so you need to specify the index based on the information:
print(test_ccd_01_hdul.info())
Filename: test_01.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 6 (100, 100) float64
None
print(type(test_ccd_01_hdul[0]))
print(type(test_ccd_01_hdul["PRIMARY"]))
<class 'astropy.io.fits.hdu.image.PrimaryHDU'>
<class 'astropy.io.fits.hdu.image.PrimaryHDU'>
As can be seen, you can use index or Name in the .info() result.
print(test_ccd_01_hdul[0].data)
[[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]
...
[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]
[1. 1. 1. ... 1. 1. 1.]]
test_ccd_01_hdul[0].header
SIMPLE = T / conforms to FITS standard
BITPIX = -64 / array data type
NAXIS = 2 / number of array dimensions
NAXIS1 = 100
NAXIS2 = 100
EXTEND = T
Now you have identical result as the above CCDData.read case.
You can also use other methods (see File Handling and Convenience Functions). Here are a few examples:
getheader(): Retrieves only the header information from a FITS file.getdata(): Retrieves only the data from a FITS file.
But if you need both header and data, DON’T use these, but just use .open(). The getheader() and getdata() methods are useful only when you specifically require one of the two separately.
NOTE: It is important to note that using methods like
setval(),delval(), and similar methods to modify the header should be done with caution. These methods can overwrite existing information in the header. If you later discover a mistake or bug in your code, you would need to re-download the FITS files to start over again. Therefore, it is advisable to use these methods only when you have a definite reason and are fully aware of the consequences of modifying the header information.
Adding/Manipulating Header#
There are several ways to add header using Header.
Header.fromstring()#
The most primitive way. Likely this is used only if you want to make a fixed example (such as bug report) or tutorial like this notebook.
test_data_02 = np.ones((10, 20))
header_str = '''OBJECT = 'dark '
GAIN = 1.360000014305115 / [e-/ADU] The electron gain factor.
RDNOISE = 9.0 / [e-] The (Gaussian) read noise.
'''
header = fits.Header.fromstring(header_str, sep='\n')
test_ccd_02 = CCDData(data=test_data_02, header=header, unit='adu')
test_ccd_02.header
OBJECT = 'dark '
GAIN = 1.360000014305115 / [e-/ADU] The electron gain factor.
RDNOISE = 9.0 / [e-] The (Gaussian) read noise.
test_ccd_02.write("test2.fits", overwrite=True)
test_ccd_02_read = CCDData.read("test2.fits")
test_ccd_02_read.header
SIMPLE = T / conforms to FITS standard
BITPIX = -64 / array data type
NAXIS = 2 / number of array dimensions
NAXIS1 = 20
NAXIS2 = 10
OBJECT = 'dark '
GAIN = 1.360000014305115 / [e-/ADU] The electron gain factor.
RDNOISE = 9.0 / [e-] The (Gaussian) read noise.
BUNIT = 'adu '
There are other similar methods, .fromtextfile() (when you have header in text file) and fromkeys() (when you have header in something like dict). See above Header doc.
To set header keyword, value, and comment in a Header object, you have several options available:
header[<key>] = <value>: Assigns a value to a specific keyword in the header.header[<key>] = (<value>, <comment>): Assigns a value and comment to a specific keyword in the header.header.set(<key>, <value>, [<comment>])method: Sets the value and an optional comment for a specific keyword in the header.header.add_comment(<comment_message>)method: Adds a comment using theCOMMENTkeyword. Comments are typically used to provide additional information, references, or explanations that are not crucial as key-value pairs.header.add_history(<history>)method: Adds a history record using theHISTORYkeyword. History records are often used to document the actions or operations performed on the data.
Both COMMENT and HISTORY entries will appear at the end of the header (unless specified), regardless of the order in which they were added.
header = fits.Header.fromstring('') # empty header
# 1) Basic setting
header["object"] = 'dark' # header key will automatically be capitalized
# 2) With comment
header["GAIN"] = (1.36, "[e-/ADU] The electron gain factor.")
# 3) Using .set()
header.set("RDNOISE", 9.0, "[e-] The (Gaussian) read noise.")
header.set("RDNOISE", 10.0, "[e-] The (Gaussian) read noise. Oops, I am adding it again! What will happen?")
# 4) Adding COMMENT line
header.add_comment("This is a testing fits file.")
header.add_comment("This is the second comment.")
header.add_comment("What if the comment is too long? "*6)
# 5) Adding HISTORY line
header.add_history("Bias subtracted 2019-01-01T00:00:01")
header.add_history("Dark corrected 2019-01-01T00:00:02")
header.add_history("Flat corrected 2019-01-01T00:00:03")
header.add_history("Cosmic ray rejected 2019-01-01T00:00:04")
header.add_history("WCS added 2019-01-01T00:00:06")
test_ccd_03 = CCDData(data=test_data_02, header=header, unit='adu')
test_ccd_03.header
WARNING: VerifyWarning: Card is too long, comment will be truncated. [astropy.io.fits.card]
OBJECT = 'dark '
GAIN = 1.36 / [e-/ADU] The electron gain factor.
RDNOISE = 10.0 / [e-] The (Gaussian) read noise. Oops, I am addi
COMMENT This is a testing fits file.
COMMENT This is the second comment.
COMMENT What if the comment is too long? What if the comment is too long? What i
COMMENT f the comment is too long? What if the comment is too long? What if the
COMMENT comment is too long? What if the comment is too long?
HISTORY Bias subtracted 2019-01-01T00:00:01
HISTORY Dark corrected 2019-01-01T00:00:02
HISTORY Flat corrected 2019-01-01T00:00:03
HISTORY Cosmic ray rejected 2019-01-01T00:00:04
HISTORY WCS added 2019-01-01T00:00:06
From the result in the RDNOISE example, you can observe that the header is indeed overwritten by the newer value and comment when using the mentioned methods. If a comment is too long, it will be truncated as indicated by the WARNING message.
Regarding the long comment in the COMMENT example, you can see that a line break occurs when the comment exceeds a certain length. This behavior ensures that the comment fits within the maximum allowed width for each line in the FITS header.
test_ccd_03.write("test3.fits", overwrite=True)
test_ccd_03_read = CCDData.read("test3.fits")
test_ccd_03_read.header
SIMPLE = T / conforms to FITS standard
BITPIX = -64 / array data type
NAXIS = 2 / number of array dimensions
NAXIS1 = 20
NAXIS2 = 10
OBJECT = 'dark '
GAIN = 1.36 / [e-/ADU] The electron gain factor.
RDNOISE = 10.0 / [e-] The (Gaussian) read noise. Oops, I am addi
BUNIT = 'adu '
COMMENT This is a testing fits file.
COMMENT This is the second comment.
COMMENT What if the comment is too long? What if the comment is too long? What i
COMMENT f the comment is too long? What if the comment is too long? What if the
COMMENT comment is too long? What if the comment is too long?
HISTORY Bias subtracted 2019-01-01T00:00:01
HISTORY Dark corrected 2019-01-01T00:00:02
HISTORY Flat corrected 2019-01-01T00:00:03
HISTORY Cosmic ray rejected 2019-01-01T00:00:04
HISTORY WCS added 2019-01-01T00:00:06
Making/Reading MEF (optional)#
If you are not advanced FITS user, it’s generally better NOT to think about making MEF. It’s better to stick to single-extension FITS format.
But for many cases, e.g., HST images, may contain MEF, because there are clear reasons. Thus, I add how to deal with MEF here.
prim = fits.PrimaryHDU(data=None, header=header)
im1 = fits.ImageHDU(data=np.ones((10, 10)), header=None, name="SCI")
im2 = fits.ImageHDU(data=np.zeros((10, 10)), header=None, name="UNCERT")
hdul_mef = fits.HDUList([prim, im1, im2])
hdul_mef.writeto("test_mef.fits", overwrite=True, output_verify='fix')
hdul_mef_read = fits.open("test_mef.fits")
print(hdul_mef_read.info())
Filename: test_mef.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 17 ()
1 SCI 1 ImageHDU 8 (10, 10) float64
2 UNCERT 1 ImageHDU 8 (10, 10) float64
None
Each extension can be accessed by either (1) index (No.) or (2) Name.
np.testing.assert_allclose(hdul_mef_read[1].data, hdul_mef_read["SCI"].data)
The getdata() or getheader(), etc, allows the arguments ext and extname to select only one extension:
hdul_mef_data_1 = fits.getdata("test_mef.fits", ext=1)
hdul_mef_data_2 = fits.getdata("test_mef.fits", extname="UNCERT")
np.testing.assert_allclose(hdul_mef_data_1, hdul_mef_data_1 - hdul_mef_data_2)
Note that, CCDData understands any extension named 'UNCERT' as the standard deviation uncertainty map (or simple 1-sigma error map), and reads it as astropy.nddata.nduncertainty.StdDevUncertainty. Similarly, 'MASK' is understood as mask. This is a reason why HST uses MEF than single-extension.
You can see this in our example:
test_ccd_mef_read = CCDData.read("test_mef.fits", unit='adu')
INFO: first HDU with data is extension 1. [astropy.nddata.ccddata]
test_ccd_mef_read.data
array([[1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
[1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]])
print(type(test_ccd_mef_read.uncertainty))
test_ccd_mef_read.uncertainty
<class 'astropy.nddata.nduncertainty.StdDevUncertainty'>
StdDevUncertainty([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])
tmpfitsfiles = Path('.').glob("test*.fits")
for fpath in tmpfitsfiles:
fpath.unlink()