Astropy版本： 5.0.1

FITS(Flexible Image Transport System)，为天文学界中通用的数据传输格式，由IAU指定了统一的标准。而著名的天文工具包Astropy中就提供了读取与处理FITS文件的工具，可将读取的FIT转化为ndarray在程序中进行后续处理。

astropy.io.fits的文档请参见FITS File Handling (astropy.io.fits) — Astropy v5.0.1

本文主要关注具体读取FITS的方法，FITS文件的其他介绍可参照其他资料。

1. FITS文件基本操作

FITS文件打开：

1
2
3
4
from astropy.io import fits

IO = './...'
hdulist = fits.open(IO +'155_0.fits')

此时即可查阅FITS文件中的基本信息。

1
2
3
4
5
6
7
8
# 查看文件基本信息
print(hdulist.info())

Filename: ./155_0.fits 
No.   Name Ver    Type     Cards Dimensions      Format
0   PRIMARY 1   PrimaryHDU   4       () 
1           1   BinTableHDU  21    40R x 6C [J, J, J, J, J, 41943E] 
None

可以看到这个FITS文件中中有两组数据，可进行提取查看，一般第一组（标号为0）数据只有头文件，没有数据：

1
2
3
4
5
6
7
>>>hdulist[0].header

SIMPLE =                T / conforms to FITS standard 
BITPIX =                8 / array data type 
NAXIS =                 0 / number of array dimensions 
EXTEND =                T 
END

这是标题记录文件，记录该FITS文件的基本情况。各行释义如下：

  • SIMPLE=T：T为逻辑值True，意为该文件的数据记录与标题记录都符合FITS文件格式约定。
  • BITPIX=8：指明每一个值在内存中所占位数为8。
  • NAXIS=0：数据维数，该处说明原数据只有一维。若数据为多维数组，还会有NAXIS1、NAXIS2等说明各维度的数据大小。
    不同的FITS文件的标题记录文件中会有不同信息，可查看相应注释了解。
    查看第二组数据可得到数据的相关信息，其中具体的观测数据也储存在第二组数据中。
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    >>>hdulist[1].header

    XTENSION= 'BINTABLE' / binary table extension 
    BITPIX = 8 / array data type 
    NAXIS = 2 / number of array dimensions 
    NAXIS1 = 167792 / length of dimension 1 
    NAXIS2 = 40 / length of dimension 2 
    PCOUNT = 0 / number of group parameters 
    GCOUNT = 1 / number of groups 
    TFIELDS = 6 / number of table fields 
    TTYPE1 = 'EXPOSURE' #参数名
    TFORM1 = 'J ' #参数类型
    TUNIT1 = 's ' #参数单位
    TTYPE2 = 'NCHAN ' 
    TFORM2 = 'J ' 
    TTYPE3 = 'CHAN_BW ' 
    TFORM3 = 'J ' 
    TTYPE4 = 'SPAN ' 
    TFORM4 = 'J ' 
    TTYPE5 = 'TIMESPAN' 
    TFORM5 = 'J ' 
    TTYPE6 = 'DATA ' 
    TFORM6 = '41943E '

头文件中各项含义在上文与注释中介绍。可根据头文件中的信息读取数据与各个参数：

1
2
3
4
5
6
data1 = hdulist[1].data #数据
header1 = hdulist[1].header

nline = header1['NAXIS2']#行数
dt = data1['EXPOSURE'][0] #采样时间
...

data1中的数组可化为ndarray数组以便进行后续处理。在实际使用过程中，因为FITS文件比较大，若在打开后不及时关闭，在后续处理中会导致电脑卡死（笔者血泪经验），所以在进行FITS文件操作时一定要养成随手关文件的好习惯：

1
hdulist.close()#关闭fit文件

在批量处理大量FITS文件时，可用下列方法打开文件，更安全：

1
2
3
4
5
6
7
8
9
10
11
for i in tqdm(range(1,41)):
	if i < 10:
		with fits.open(IO_data + filename + '0' + str(i) + '.fit') as hdul:
			data = (np.array(hdul[0].data) - bia_mean) * np.mean(flat_mean - bia_mean) / (flat_mean - bia_mean)
			time_fit.append(hdul[0].header['DATE-OBS'][11:])#记录时间
			header = hdul[0].header
			wcs = WCS(header)
			#提取要使用的数据
			hdul.close()#关闭文件
	else:
....

即在with...后提取出后续要用到的数据中马上关闭该fit文件，节省内存。

2. FITS文件画图

有些FITS文件数据里会自带wcs坐标，在使用该数据画图时可利用wcs坐标使图像附带坐标信息：

1
2
3
4
5
6
7
8
9
from astropy import wcs
w = wcs.WCS(hdulist[0].header)

plt.figure(figsize=(8,6))
plt.subplot(projection=w)
im = plt.imshow(data,cmap='CMRmap')
plt.colorbar(im)
#plt.xlabel('Right Ascension')
plt.ylabel('Declination')

2022/3/21更新
若要画多张子图，可在设置子图时添加以下参数：

1
fig,ax = plt.subplots(1,2,figsize=(10,8),subplot_kw={'projection': wcs})

即可为每一张子图添加wcs坐标。
（感谢python - How do I change matplotlib’s subplot projection of an existing axis? - Stack Overflow）

对于FITS文件处理，笔者也只是刚刚入门，此后会根据学习进度继续更新本文。