FITS（Flexible Image Transport System）[1]是一种用来存储、传输和操作天文科学数据和相关信息的通用格式，这些数据可以是一维的光谱、二维的图像、三维或更高维的数据块，也可以是ASCII或二进制的表格。FITS同时会包含数据的观测信息和校准信息，以方便后续的数据处理。
一个FITS文件可能包括一组或多组HDU（Header Data Unit），每一组HDU都包括一个头文件和一个数据单元。第一个HDU被称为主HDU（primary HDU），其它的HDU则被称为FITS扩展（FITS extension），FITS文件的主HDU也有可能是空的，只有扩展HDU中包含数据。
FITS包括三种基本的数据类型：
图像扩展（主HDU中也可以是图像数据）：N维数据块，使用统一的数据类型（一共有六种数据类型可供使用：8位无符号整型，16、32、64位有符号整型，32、64位浮点数）。
ASCII表格扩展：可以用来储存星表或天文信息列表，每一行都是固定的长度，并被分成固定长度的不同部分，可以存放ASCII字符串，或者字符串格式的数字。
二进制表格扩展：比ASCII表格提供了更多的数字存储特性，即可以存储可变长度的数据，并且占用的空间更小。
FITS的头文件或者数据单元的大小必须是2880字节的整数倍，末尾未使用的空间要用0或者空格来补齐。头文件由很多长度为80字符的关键字记录组成，通用格式为关键字 = 数值 / 注释（KEYNAME = value / comment）。头文件会记录对应数据单元的大小和格式这样的必要信息，以及观测的日期和时间等信息，关键字COMMENT和HISTORY也经常被用来进一步记录数据的内容。
FITS图像数据的直接查看和简单分析可以通过DS9[2]来完成。DS9中显示的FITS文件的数据部分和头文件信息分别如图所示（这个FITS来自HASH，显示了一个Halpha波段的行星状星云）：
本教程所需全部Python包：
# 本教程所需的全部Python包
from astropy.table import Table
from astropy.io import fits
import numpy as np1. 读取FITS：
[3]在Python中，我们使用astropy.io.fits来进行FITS文件的操作。FITS的读取主要有两种形式：
1.1 astropy.io.fits.open()
可以使用fits.open()直接打开FITS文件，并将HDU数据储存在变量hdu中。然后通过使用hdu[0]来索引主HDU（如果有扩展HDU，依次为1,2,3……），并用.header和.data来获取其中的头文件和数据。得到的头文件和数据与hdu.info()输出的信息一致。
读取结束后需要使用hdu.close()关闭文件。（好像很少留意到有人这样做，但还是要尽量符合代码规范。）
读取FITS并输出HDU信息（如果已经import所需的包，不需要重复操作）：
from astropy.io import fits
# 这个函数仅用来搜索astropy自带的FITS示例，这一行可以给定自己的路径和文件名fits_image_filename = fits.util.get_testdata_filepath('test0.fits')
# fits_image_filename = '/path/to/your/file.fits'
# 打开FITS文件
hdu = fits.open(fits_image_filename)
# 输出HDU信息
hdu.info()输出：
Filename: /opt/anaconda3/lib/python3.9/site-packages/astropy/io/fits/tests/data/test0.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 138 ()
1 SCI 1 ImageHDU 61 (40, 40) int16
2 SCI 2 ImageHDU 61 (40, 40) int16
3 SCI 3 ImageHDU 61 (40, 40) int16
4 SCI 4 ImageHDU 61 (40, 40) int16读取头文件和数据单元，并关闭文件：
# 读出头文件
hd = hdu[1].header
print(len(hd),type(hd))
# 读出数据
data = hdu[1].data
print(len(data),type(data))
# 关闭文件
hdu.close()输出：
61 <class 'astropy.io.fits.header.Header'>
40 <class 'numpy.ndarray'>如果不想使用hdu.close()来关闭文件，可以使用如下方法，退出with时文件会被自动关闭：
# 如果不想手动关闭文件，可以使用如下方式
with fits.open(fits_image_filename) as hdu:
hdu.info()输出：
Filename: /opt/anaconda3/lib/python3.9/site-packages/astropy/io/fits/tests/data/test0.fits
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 138 ()
1 SCI 1 ImageHDU 61 (40, 40) int16
2 SCI 2 ImageHDU 61 (40, 40) int16
3 SCI 3 ImageHDU 61 (40, 40) int16
4 SCI 4 ImageHDU 61 (40, 40) int16 1.2 一些便捷的函数
astropy.io.fits也提供了一些方便调用的函数，这样可以直接获得所需的数据，而且不需要关闭文件。但是需要注意的是，这些函数的效率比较低，不适合在复杂的分析脚本和应用程序中使用。近期看到很多程序都直接使用了类似的方法，在应用到大数据量的时候要特别注意。
以下通过几个例子来简单介绍用来获取整个头文件的getheader()、直接获取头文件关键字的getval()、获取数据信息的getdata()的使用：
from astropy.io import fits
# 获取第二个'SCI'扩展的头文件
hd = fits.getheader(fits_image_filename, extname='sci', extver=2)
print(len(hd),type(hd))
# 获取主HDU的观测时间
flt = fits.getval(fits_image_filename, 'date', 0)
print(flt)
# 获取第三个'SCI'扩展的数据
data = fits.getdata(fits_image_filename, 'sci', 3)
print(len(data),type(data))
# 获取第一个扩展HDU的数据和头文件
data, hd = fits.getdata(fits_image_filename, 1, header=True)
print(len(data),type(data))
print(len(hd),type(hd))输出：
61 <class 'astropy.io.fits.header.Header'>
01/04/99
40 <class 'numpy.ndarray'>
40 <class 'numpy.ndarray'>
61 <class 'astropy.io.fits.header.Header'>2. 处理FITS头文件：
头文件包括很多行，每一行都是80个字节，被分为三部分：关键字、数值和注释。其中关键字和注释必须是字符串，数值可以是字符串、整数、浮点数、复数或者布尔值（True/False）。
2.1 读取头文件信息
头文件中记录了数据拍摄和校准的重要信息，除了COMMENT和HISTORY等特殊例子，关键字一般是不会重复的，所以可以根据关键字来调用数值，另外也可以通过索引/行数来调用。这里的索引是按Python的惯例，从0开始计数的。
一般来说，使用关键字来调用比较合理，因为每个关键字所在的行数可能是不固定的，使用行数来调用可能会出现问题，尤其是在修改数值的时候要小心。
头文件中的关键字一般是大写的，但是调用的时候对大小写不敏感。如果这个头文件中不含有所调用的关键字，就会产生报错。
通过关键字和索引来读取头文件信息，以及读出头文件信息的注释：
from astropy.io import fits
# 打开FITS文件
hdu = fits.open(fits_image_filename)
# 通过关键字读取主HDU的日期
date = hdu[0].header['DATE']
# 通过索引来读取主HDU的日期
date1 = hdu[0].header[9]
# 读出头文件注释
hd = hdu[0].header
cmt = hd.comments['DATE']
print(date,date1,cmt)输出：
01/04/99 01/04/99 Date FITS file was generated2.2 修改头文件信息
可以用直接赋值的形式来修改头文件信息，如果关键字不存在，这个操作会自动在头文件末尾添加一行，如果是COMMENT和HISTORY的话，则会在对应的关键字末尾添加一行。
另一种修改或新增头文件信息的方法是header.set()。
修改或新增头文件信息：
# 读出头文件
hd = hdu[0].header
# 用关键字修改头文件信息，如果关键字不存在，就会在头文件末尾添加
hd['targname'] = 'NGC121-a'
# 以元组的形式，同时修改数值和注释
hd['targname'] = ('NGC121-a', 'the observation target')
# 用索引修改头文件信息
hd[27] = 99
# 修改或新增关键字的另一种方法
hd.set('observer', 'Edwin Hubble')
print(hd['targname'],hd[27],hd['observer'])输出：
NGC121-a 99 Edwin Hubble如果对COMMENT和HISTORY关键字进行同样的操作，只会增加一行，而不会修改已有的信息。如果想修改已有的信息，需要使用索引。
另外，不要混淆COMMENT关键字和每一行的注释（英文也是comment）。
添加和修改HISTORY关键字（如果重复运行这几行代码，会输出四行，读者可以思考一下为什么会这样）：
# 添加history关键字
hd['history'] = 'I updated this file on 2/26/09'
print(hd['history'])
# 修改history关键字
hd['history'][0] = 'I updated this file on 2/27/09'
print(hd['history'])输出：
I updated this file on 2/26/09
I updated this file on 2/27/092.3 显示头文件信息
由于很多头文件是不标准的，在调用前有时需要查看所需的关键字，或打印整个头文件来查看。除了用hdu.info()输出HDU信息，以及用关键字和索引来引用单行的头文件信息，还有更多的显示头文件的方法，这里再介绍三种：
使用list(hd.keys())可以获取头文件的关键字，快速找到需要使用的关键字：
# 获取头文件关键字
print(list(hd.keys()))输出：
['SIMPLE', 'BITPIX', 'NAXIS', 'EXTEND', 'GROUPS', 'NEXTEND', 'BSCALE', 'BZERO', 'ORIGIN', 'DATE', 'IRAF-TLM', '', '', '', '', '', '', '', '', '', '', '', '', 'INSTRUME', 'ROOTNAME', 'FILETYPE', '', '', 'MODE', 'SERIALS', '', '', 'IMAGETYP', 'CDBSFILE', 'PKTFMT', '', '', 'FILTNAM1', 'FILTNAM2', 'FILTER1', 'FILTER2', 'FILTROT', 'LRFWAVE', '', '', 'UCH1CJTM', 'UCH2CJTM', 'UCH3CJTM', 'UCH4CJTM', 'UBAY3TMP', 'KSPOTS', 'SHUTTER', 'ATODGAIN', '', '', 'MASKCORR', 'ATODCORR', 'BLEVCORR', 'BIASCORR', 'DARKCORR', 'FLATCORR', 'SHADCORR', 'DOSATMAP', 'DOPHOTOM', 'DOHISTOS', 'OUTDTYPE', '', '', 'MASKFILE', 'ATODFILE', 'BLEVFILE', 'BLEVDFIL', 'BIASFILE', 'BIASDFIL', 'DARKFILE', 'DARKDFIL', 'FLATFILE', 'FLATDFIL', 'SHADFILE', 'PHOTTAB', 'GRAPHTAB', 'COMPTAB', '', '', 'SATURATE', 'USCALE', 'UZERO', '', '', 'READTIME', '', '', 'PA_V3', 'RA_SUN', 'DEC_SUN', 'EQNX_SUN', 'MTFLAG', 'EQRADTRG', 'FLATNTRG', 'NPDECTRG', 'NPRATRG', 'ROTRTTRG', 'LONGPMER', 'EPLONGPM', 'SURFLATD', 'SURFLONG', 'SURFALTD', '', '', 'PODPSFF', 'STDCFFF', 'STDCFFP', 'RSDPFILL', '', '', 'UEXPODUR', 'NSHUTA17', 'DARKTIME', 'UEXPOTIM', 'PSTRTIME', 'PSTPTIME', '', '', '', '', 'ORIENTAT', 'SUNANGLE', 'MOONANGL', 'SUN_ALT', 'FGSLOCK', '', 'DATE-OBS', 'TIME-OBS', 'EXPSTART', 'EXPEND', 'EXPTIME', 'EXPFLAG', 'FILENAME', 'TARGNAME', 'OBSERVER', 'HISTORY']使用切片来查看多行头文件（可以参考Python的切片的使用方法）：
# 显示头文件前两行，不分行，这里的0可以省略
print(hd[0:2])
# 显示头文件前两行，分行
print(repr(hd[0:2]))输出（这里调小了字号以显示格式的区别）：
SIMPLE = T / file does conform to FITS standard BITPIX = 16 / number of bits per data pixel END
SIMPLE = T / file does conform to FITS standard
BITPIX = 16 / number of bits per data pixel 用合适的格式查看整个头文件（如果直接打印整个头文件，格式会比较混乱，加上repr()之后，能比较清楚地分行显示）：
# 显示整个头文件，不加repr()也可以显示，但是不会分行
print(repr(hd))部分输出：
SIMPLE = T / file does conform to FITS standard
BITPIX = 16 / number of bits per data pixel
NAXIS = 0 / number of data axes
EXTEND = T / FITS dataset may contain extensions
GROUPS = F / data has groups
NEXTEND = 4 / Number of standard extensions
BSCALE = 1.000000E0 / REAL = TAPE*BSCALE + BZERO
BZERO = 3.276800E4 /
ORIGIN = 'NOAO-IRAF FITS Image Kernel Aug 1 1997' / FITS file originator
DATE = '01/04/99 ' / Date FITS file was generated
IRAF-TLM= 'xxx ' / Time of last modification
......3. 处理FITS数据：
FITS的数据单元是FITS中重点处理的部分，后面的文章还会进行更多的介绍，这里只介绍一些基本操作。
前文已经介绍过如何简单地用索引来读取数据单元，还可以使用拓展名（EXTNAME）来引用，如果有多个同名的拓展，第一个拓展仍然可以只使用拓展名来引用，但是需要同时使用索引（EXTVER）来引用后面的拓展。（NAME和VER见hdu.info()的输出。）
在这个例子里，主HDU是hdu[0]，第一个拓展HDU是hdu[1]或hdu['sci']或hdu['sci',1]，第二个拓展HDU是hdu[2]或hdu['sci',2]，以此类推。这里的'SCI'也是大小写不敏感的。
另外，要特别注意FITS的x和y与Python中是相反的，可以自己检查一下。
3.1 处理图像数据
读出图像数据并获得数据形状、数据类型、像素值的信息：
# 使用拓展名和索引来读出数据，如果不引起歧义，可以省略其中之一
data = hdu['SCI',2].data
# 获得数据的形状、数据类型信息
print(data.shape,data.dtype.name)
# 查看x=5，y=2的像素值
print(data[1, 4])
# 查看x=11到20，y=31到40的像素值
print(data[30:40, 10:20])输出：
(40, 40) int16
348
[[350 349 349 348 349 348 349 347 350 348]
[348 348 348 349 348 349 347 348 348 349]
[348 348 347 349 348 348 349 349 349 349]
[349 348 349 349 350 349 349 347 348 348]
[348 348 348 348 349 348 350 349 348 349]
[348 347 349 349 350 348 349 348 349 347]
[347 348 347 348 349 349 350 349 348 348]
[349 349 350 348 350 347 349 349 349 348]
[349 348 348 348 348 348 349 347 349 348]
[349 349 349 348 350 349 349 350 348 350]]3.2 处理表格数据
我们经常用Pandas来处理表格数据，所以这里也只是简单介绍一下astropy.io.fits中的一些基本的处理方法。
表格的每一行都是一个FITS_record对象，由不同数据类型的元素组成，所以表格就是由这些记录组成的数组。可以通过索引来读取表格中的行，或者通过列名来读出表格中的列，也可以使用field()方法，通过索引或者列名来获得表格的每一列：
from astropy.io import fits
# 读取表格数据
fits_table_filename = fits.util.get_testdata_filepath('tb.fits')
hdu = fits.open(fits_table_filename)
data = hdu[1].data
# 使用索引来读出表格的第一行
print(data[0])
# 使用列名来读出表格的第一列
print(data['c1'])
# field()方法，使用索引读出表格的第一列
print(data.field(0))
# field()方法，用列名读出表格的第一列
print(data.field('c1'))
输出：
(1, 'abc', 3.7000000715255736, False)
[1 2]
[1 2]
[1 2]我们更推荐使用列名来获得表格中的列，这样不会受到顺序的影响，另外，列名也是大小写不敏感的。如何知道表格中有哪些列名，以及获得其它信息：
# 获得表格中列的信息
cols = hdu[1].columns
print(cols)
# 获得某一项列信息，比如列名
print(cols.names)输出：
ColDefs(
name = 'c1'; format = '1J'; null = -2147483647; disp = 'I11'
name = 'c2'; format = '3A'; disp = 'A3'
name = 'c3'; format = '1E'; bscale = 3; bzero = 0.4; disp = 'G15.7'
name = 'c4'; format = '1L'; disp = 'L6'
)
['c1', 'c2', 'c3', 'c4']另一种获得表格中的列信息的方法：
# 另一种获得表格中的列信息的方法
print(cols.info())部分输出：
name:
['c1', 'c2', 'c3', 'c4']
format:
['1J', '3A', '1E', '1L']
unit:
['', '', '', '']
null:
[-2147483647, '', '', '']
bscale:
['', '', 3, '']
bzero:
['', '', 0.4, '']
disp:
['I11', 'A3', 'G15.7', 'L6']
......因为这些列都是Numpy对象，所以可以对其进行各种Numpy的操作，比如修改表格，以及对表格信息进行统计：
# 修改表格，比如将第四列赋值为0（变为False）
print(data['c4'])
data['c4'][:] = 0
print(data['c4'])
# 求第三列的平均值
print(data['c3'].mean())输出：
[False False]
[False False]
5.199999892711644.写入FITS：
4.1 写入修改的结果
在对FITS进行操作之后，我们可以对结果进行保存，仍然写入FITS中。
使用hdu.writeto()可以将内存中的数据和头文件版本写入一个新的FITS文件：
# 将内存中的数据和头文件写入新的FITS，如果文件已存在则覆盖原文件hdu.writeto('newtable.fits', overwrite=True)也可以通过添加overwrite=True，来写入原文件，但是为了保留备份，还是比较推荐写入新的文件。另外，设置这个参数也可以方便重复运行这个命令，不然会因为文件已存在而报错。
4.2 建立一个新的FITS文件
我们经常需要从头开始建立一个FITS文件，并写入需要保存的数据，这里我们用一个新建的Numpy数组来作为数据单元：
from astropy.io import fits
import numpy as np
# 建立一个从0到99的浮点数数组作为数据单元
n = np.arange(100.0)
print(n)
# 建立一个PrimaryHDU对象
hdu = fits.PrimaryHDU(n)
# 建立一个HDU列表并写入一个新文件
hdul = fits.HDUList([hdu])
hdul.writeto('new1.fits', overwrite=True)
# 也可以直接将HDU对象写入文件，等同于以上两行
hdu.writeto('new2.fits', overwrite=True)输出：
[ 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17.
18. 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35.
36. 37. 38. 39. 40. 41. 42. 43. 44. 45. 46. 47. 48. 49. 50. 51. 52. 53.
54. 55. 56. 57. 58. 59. 60. 61. 62. 63. 64. 65. 66. 67. 68. 69. 70. 71.
72. 73. 74. 75. 76. 77. 78. 79. 80. 81. 82. 83. 84. 85. 86. 87. 88. 89.
90. 91. 92. 93. 94. 95. 96. 97. 98. 99.]我们也可以新建一个新的表格文件，如果只想新建一个二进制FITS表格，不需要其它扩展，可以使用astropy.table中的Table包，然后写入FITS，这个方法会比astropy.io.fits的方法更简洁：
from astropy.table import Table
# 使用Table包，新建一个二进制表格，并写入FITS
# 建立二进制表格
t = Table([[1, 2], [4, 5], [7, 8]], names=('a', 'b', 'c'))
print(t)
# 写入FITS
t.write('table1.fits', format='fits', overwrite=True)输出：
a b c
--- --- ---
1 4 7
2 5 8使用astropy.io.fits的方法如下：
from astropy.io import fits
import numpy as np
# 新建一个表格并写入FITS
# 分别建立每一列，输入列名、数据和格式
c1 = fits.Column(name='a', array=np.array([1, 2]), format='K')
c2 = fits.Column(name='b', array=np.array([4, 5]), format='K')
c3 = fits.Column(name='c', array=np.array([7, 8]), format='K')
# 建立二进制表格HDU
t = fits.BinTableHDU.from_columns([c1, c2, c3])
print([c1, c2, c3])
print(t)
# 写入FITS
t.writeto('table2.fits', overwrite=True)输出：
[name = 'a'; format = 'K', name = 'b'; format = 'K', name = 'c'; format = 'K']
<astropy.io.fits.hdu.table.BinTableHDU object at 0x7fd980c7bbe0>注意表格只能作为扩展HDU，而不能作为主HDU，所以新建的FITS的主HDU中不会包含数据，读者可以自行用本文开头介绍的hdu.info()方法来查看。
4.3 建立多扩展的FITS
同一个FITS中可以同时包括不同类型、不同大小的主HDU和扩展HDU：
from astropy.io import fits
import numpy as np
# 新建一个多扩展的FITS
# 新建图像数组
n = np.ones((3, 3))
n2 = np.ones((100, 100))
# 新建头文件，并添加一些信息
hdr = fits.Header()
hdr['OBSERVER'] = 'Edwin Hubble'
hdr['COMMENT'] = "Here's some commentary about this FITS file."
# 建立HDU对象（主HDU和图像扩展HDU）
primary_hdu = fits.PrimaryHDU(n, header=hdr)
image_hdu = fits.ImageHDU(n2)
# 建立二进制表格HDU对象（同上）
c1 = fits.Column(name='a', array=np.array([1, 2]), format='K')
c2 = fits.Column(name='b', array=np.array([4, 5]), format='K')
c3 = fits.Column(name='c', array=np.array([7, 8]), format='K')
table_hdu = fits.BinTableHDU.from_columns([c1, c2, c3])
# 建立HDU列表
hdul = fits.HDUList([primary_hdu, image_hdu, table_hdu])
# 也可以追加更多的HDU
n3 = np.ones((10, 10, 10))
image_hdu2 = fits.ImageHDU(n3)
hdul.append(image_hdu2)
print(hdul.info())
# 写入FITS
hdul.writeto('new3.fits', overwrite=True)输出：
Filename: (No file associated with this HDUList)
No. Name Ver Type Cards Dimensions Format
0 PRIMARY 1 PrimaryHDU 8 (3, 3) float64
1 1 ImageHDU 7 (100, 100) float64
2 1 BinTableHDU 14 2R x 3C ['K', 'K', 'K']
3 1 ImageHDU 8 (10, 10, 10) float64
None有时候我们也需要建立一个空的主HDU，只包含一些基本的头文件信息，这时候primary_hdu = fits.PrimaryHDU(header=hdr)只输入头文件即可。在这一步，函数会自动建立一个包含基本信息（比如SIMPLE和NAXIS等关键字）的头文件，读者只需要处理观测相关的关键字即可。
4.4 一些便捷的函数
除了之前介绍的一些获取信息的便捷函数，还有一些可以用来修改和写入FITS文件。因为不太常见，所以这里只做简单介绍：
from astropy.io import fits
# 写入FITS，这里的data和hdr是之前读出的，可以改成自己的数据fits.writeto('out.fits', data, hdr, overwrite=True)
# 追加写FITS，如果不存在，则新建，每次运行都会追加写一次
fits.append('out1.fits', data, hdr)
# 比较两个FITS版本的差别
fits.printdiff('out.fits', 'out1.fits')
# 可以比较指定的HDU的差别
fits.printdiff('out.fits', 'out1.fits', ext=0)第一次运行输出：
fitsdiff: 5.0.4
a: out.fits
b: out1.fits
Maximum number of different data values to be reported: 10
Relative tolerance: 0.0, Absolute tolerance: 0.0
No differences found.
No differences found.第二次运行输出：
fitsdiff: 5.0.4
a: out.fits
b: out1.fits
Maximum number of different data values to be reported: 10
Relative tolerance: 0.0, Absolute tolerance: 0.0
Files contain different numbers of HDUs:
a: 2
b: 3
No differences found between common HDUs.
No differences found.这里out.fits每次运行都会被重复写一次，而out1.fits每次运行都会追加写一次，于是使用printdiff可以查看HDU数量的差别，而HDU的内容并没有差别。
扩展阅读：
[3]FITS File Handling (astropy.io.fits): https://docs.astropy.org/en/stable/io/fits/index.html