astropy.io.fits 教程
astropy.io.fits 教程
**astropy.io.fits**包提供了对fits文件的访问。FITS(灵活图像传输系统)是一个可移植的文件标准，广泛用于天文存储图像和表格。
读取fits文件
Note
get_testdata_filepath()函数在这里的示例中使用，用于访问由astropy附带的数据。要使用您自己的数据，请使用astropy.io.fit .open()，它采用相对或绝对路径。
如果已经安装了astropy模块，就可以打开一个已经存在的fits文件:
from astropy.io import fits
fits_image_filename = fits.util.get_testdata_filepath('test0.fits')
hdul = fits.open(fits_image_filename)
open()函数有几个可选参数，我们将在后面的章节中讨论。在上面的例子中，默认模式是“readonly”。open函数返回一个名为HDUList的对象，它是一个类似列表的HDU对象集合。HDU(头数据单元)是FITS文件结构的最高级别组件，由一个头和(通常)一个数据数组或表组成。
HDUList有一个有用的方法HDUList.info()，它总结了打开的FITS中的内容:
from astropy.io import fits
fits_image_filename = fits.util.get_testdata_filepath('test0.fits')
hdul = fits.open(fits_image_filename)
hdu_info = hdul.info()
print(hdu_info)
当对fits文件所有处理都完成后，可以用**HDUList.close()**方法关闭它:
hdul.close()
也有一种方法可以避免特定语句关闭fits文件：
with fits.open(fits_image_filename) as hdul:
hdu_info = hdul.info()
print(hdu_info)
退出with范围后，文件将自动关闭。这(通常)是Python中打开文件的首选方式，因为即使发生异常，它也会关闭文件。
如果文件是用lazy_load_hdus=False打开的，那么在关闭HDUList之后，仍然可以访问所有的头文件。标题和数据可能或可能不可访问取决于是否数据被触摸和如果他们是内存映射;详情请参阅后面的章节。
处理较大的fits文件
open()函数支持memmap=True参数，该参数允许使用mmap访问每个HDU的数组数据，而不是一次全部读入内存。这对于处理不能完全装入物理内存的非常大的数组特别有用。模块中默认memmap=True，该值是从配置项astropy.io.fit . conf .use_memmap获得的。
这对较小的文件的影响也很小，尽管一些操作(比如顺序读取数组数据)可能会带来一些额外的开销。在32位系统上，大于2到3 GB的数组不能被mmap(这很好，因为到那个时候无论如何都有可能耗尽物理内存)，但是64位系统在这方面的限制要小得多。
当memmap=True打开文件时，由于mmap的工作方式，这意味着当HDU数据被访问时(例如，hdul[0].data)，另一个FITS文件的句柄被mmap打开。这意味着，即使在调用了hdul.close()之后，mmap仍然持有一个开放的数据句柄，这样，在假定.data属性拥有内存中的所有数据的情况下构建的不小心的程序仍然可以访问该数据。
为了强制关闭mmap，要么等待包含的HDUList对象超出作用域，要么手动调用del hdul[0].data。(只要没有对数据数组的其他引用，这种方法就可以工作。
无符号整数 （Unsigned integers）
由于FITS格式的Fortran起源，FITS在图像或表中不天生支持无符号整数数据。但是，有一种常见的约定是将无符号整数存储为有符号整数，同时使用一个移位指令(一个值为2 ** (BITPIX - 1)的BZERO关键字)将所有有符号整数上移为无符号整数。例如，当将值0写成无符号的32位整数时，它将在FITS文件中存储为-32768，以及头关键字BZERO = 32768。
默认情况下，astropy识别并应用这个约定，因此，所有看起来应该被解释为无符号整数的数据都将自动转换(这适用于图像和表)。在v1.1.0之前的astropy版本中，这不会自动应用，因此必须将参数uint=True传递给open()。在v1.1.0或更高版本中，这是默认值。
即使使用uint=False，仍然应用了BZERO移位，但是返回的数组是“float64”类型。要完全禁用缩放/移动，使用do_not_scale_image_data=True(看看为什么包含整数数据的图像被意外地转换为浮点数?)详情请参阅常见问题解答)。
压缩后的fits文件
**open()**函数将无缝地打开使用gzip、bzip2或pkzip压缩的FITS文件。注意，在本文中，我们讨论的是用这些工具之一压缩的FITS文件(例如.fit .gz文件)。
使用压缩文件时有一些限制。例如，对于包含多个压缩文件的Zip文件，只有第一个文件是可访问的。另外，bzip2不支持追加或更新访问模式。
当写入一个文件(例如，使用**writeto()**函数)时，压缩将根据给定的文件名扩展名或正在写入的已存在文件中使用的压缩来确定。
处理非标准的fits文件
当astropy.io.fits将读取一个不符合fits标准的fits文件，它将试图对不符合标准的字段做出有意义的解释。这可能不会每次都成功，在访问头文件时可能会触发警告，在写入文件时可能会触发异常。可以使用**open()**的output_verify参数来控制写入输出文件的字段的验证。打开读取的文件可以用hdulist方法进行验证和修复。这个方法通常在打开文件后调用，但在访问任何头文件或数据之前:
>>> with fits.open(fits_image_filename) as hdul:
... hdul.verify('fix')
... data = hdul[1].data
在上面的示例中，调用astropy.io中的hdul.verify(“fix”),适合修复不兼容的字段并打印信息消息。在FITS验证中描述了除“fix”外的其他选项
修改fits文件
处理FITS 的 HEADER (表头)
如前所述，HDUList的每个元素都是一个带有.header和.data属性的HDU对象，可用于访问HDU的header和数据部分。
对于FITS的HERDER，它们由一个80字节的“卡片”列表组成，其中卡片包含一个关键字、一个值和一个注释。关键字和注释都必须是字符串，而值可以是字符串或整数、浮点数、复数或True/False。关键字在标题中通常是唯一的，除了少数特殊情况。
header属性是一个header实例，另一个astropy对象。要获得与一个头关键字相关的值，需要按照Python字典操作方法:
>>> hdul = fits.open(fits_image_filename)
>>> hdul[0].header['DATE']
'01/04/99'
按照上述代码获取到关键字“DATE”的值，它是一个字符串’ 01/04/99 '。
尽管FITS文件中的关键字名称总是用大写，但为了方便用户，使用astropy指定关键字名称是不区分大小写的。如果指定的关键字名称不存在，它将引发KeyError异常。
我们也可以通过索引(Python列表)得到关键字值:
>>> hdul[0].header[7]
32768.0
这个示例返回第8个(与Python列表一样，它是0索引的)关键字的值—一个浮点数—32768.0。
同样，在astropy中也可以通过关键字名称或索引更新关键字的值:
>>> hdr = hdul[0].header
>>> hdr['targname'] = 'NGC121-a'
>>> hdr[27] = 99
请注意，几乎所有的应用程序代码都应该通过关键字名称而不是位置索引来更新标题值。这是因为大多数FITS关键字可能出现在表头的任何位置。
也有可能更新与关键字相关联的值和注释，将它们分配为元组:
>>> hdr = hdul[0].header
>>> hdr['targname'] = ('NGC121-a', 'the observation target')
>>> hdr['targname']
'NGC121-a'
>>> hdr.comments['targname']
'the observation target'
与dict一样，您也可以使用上述语法来添加一个新的keyword/value对(以及可选的注释)。在这种情况下，新卡片被追加到头的末尾(除非它是评论关键字，比如COMMENT或HISTORY，在这种情况下，它被追加到带有该关键字的最后一张卡片之后)。
更新现有卡或追加新卡的另一种方法是使用**Header.set()**方法:
>>> hdr.set('observer', 'Edwin Hubble')
Comment或history像普通卡片一样被添加，尽管在它们的情况下总是创建一个新的卡片，而不是更新一个现有的历史或评论卡片:
>>> hdr['history'] = 'I updated this file 2/26/09'
>>> hdr['comment'] = 'Edwin Hubble really knew his stuff'
>>> hdr['comment'] = 'I like using HST observations'
>>> hdr['history']
I updated this file 2/26/09
>>> hdr['comment']
Edwin Hubble really knew his stuff
I like using HST observations
注意:注意不要将注释卡与普通卡片的注释值混淆。
要更新现有的Comment或history，请按索引参考:
>>> hdr['history'][0] = 'I updated this file on 2/27/09'
>>> hdr['history']
I updated this file on 2/27/09
>>> hdr['comment'][1] = 'I like using JWST observations'
>>> hdr['comment']
Edwin Hubble really knew his stuff
I like using JWST observations
要查看FITS文件中出现的整个头(结束卡和填充被剥离)，请输入头对象本身，或print(repr(hdr)):
>>> hdr
SIMPLE = T / file does conform to FITS standard
BITPIX = 16 / number of bits per data pixel
NAXIS = 0 / number of data axes
...
>>> print(repr(hdr))
SIMPLE = T / file does conform to FITS standard
BITPIX = 16 / number of bits per data pixel
NAXIS = 0 / number of data axes
...
只输入打印(hdr)也可以，但是在大多数显示中可能不是很清楚，因为这样显示的头文件是写在FITS文件本身中的，这意味着卡片之间没有换行符。这是新用户经常感到困惑的一个原因。
它也可以查看一个切片的头:
>>> hdr[:2]
SIMPLE = T / file does conform to FITS standard
BITPIX = 16 / number of bits per data pixel
上面只显示了前两张卡片。
要获得所有关键字的列表，使用**Header.keys()**方法，就像使用dict一样:
>>> list(hdr.keys())
['SIMPLE', 'BITPIX', 'NAXIS', ...]
结构关键字 (Structural Keywords)
FITS关键字将元数据和解析文件所需的文件结构的关键信息混合在一起。这些结构关键字由astropy.io.fits内部管理。一般来说，使用者不应触摸适合的部位。相反，应该使用astropy.io.fits的相关属性。适合类(参见下面的示例)。
FITS标准使用的特定结构关键字集随HDU类型的不同而不同。下表列出了与每种HDU类型相关的关键字:
| HDU Type | Structural Keywords |
|---|---|
| All | SIMPLE, BITPIX, NAXIS |
PrimaryHDU | EXTEND |
ImageHDU, TableHDU, BinTableHDU | PCOUNT, GCOUNT |
GroupsHDU | NAXIS1, GCOUNT, PCOUNT, GROUPS |
TableHDU, BinTableHDU | TFIELDS, TFORM, TBCOL |
还有许多其他保留的关键字，例如用于数据缩放的关键字，或用于表的列属性的关键字，如FITS标准中所述。其中大多数都可以通过Column或HDU对象的属性访问，例如用HDU .name设置EXTVER的 EXTNAME或HDU.ver。结构关键字将作为通用操作的结果进行检查和/或更新。例如,当:
-
设置数据。关键字
NAXIS*从数据形状(.data.shape)中设置，BITPIX从数据类型(.data.dtype)中设置。 -
设置页眉。它的关键字根据数据属性进行更新(如上所述)。
-
写一个文件。所有必要的关键字都被删除、更新或添加到标题中。
-
调用HDU的验证方法(例如PrimaryHDU.verify())。一些关键字可以自动修复。
在这种情况下，用户可能为这些关键字分配的任何值都将被覆盖。
处理数据 DATA
fits中的数据一般分为两种，图像数据(Image Data)和表格数据(Table Data)，对于脉冲星观测产生的的数据为表格数据，所以在此并不包括对图像数据的内容。
本节描述直接使用FITS包以FITS格式读取和写入表数据。然而，在某些情况下，高级的统一文件Unified File Read/Write Interface通常就足够了，而且在某种程度上使用起来更方便。有关详细信息，请参阅Unified I/O FITS 部分。
FITS表扩展的数据部分也在.data属性中:
>>> fits_table_filename = fits.util.get_testdata_filepath('tb.fits')
>>> hdul = fits.open(fits_table_filename)
>>> data = hdul[1].data # assuming the first extension is a table
如果您熟悉numpy recarray(数组)对象，您会发现表数据基本上是一个带有一些额外属性的数组。但是熟悉数组并不是本教程的先决条件。
要查看表的第一行:
>>> print(data[0])
(1, 'abc', 3.7000000715255736, False)
表中的每一行都是一个FITS_record对象，它看起来像一个(Python)元组，包含异构数据类型的元素。在这个例子中:一个整数、一个字符串、一个浮点数和一个布尔值。因此，表数据只是这些记录的数组。更常见的情况是，用户可能以按列的方式访问数据。这是通过使用field()方法完成的。要获取表的第一列(或“字段”在NumPy的说法-它在这里与“列”互换使用)，使用:
>>> data.field(0)
array([1, 2]...)
返回具有指定字段的数据类型的numpy对象。
像标题关键字，一个列可以通过索引，如上，或通过名称:
>>> data.field('c1')
array([1, 2]...)
当通过名称访问一个列时，类似dict的访问也是可能的(甚至更好):
>>> data['c1']
array([1, 2]...)
在大多数情况下，最好通过列的名称访问列，因为列名完全独立于它在表中的物理顺序。与标题关键字一样，列名不区分大小写。
但是我们怎么知道表中有哪些列呢?首先，我们将介绍表HDU的另一个属性:columns属性
>>> cols = hdul[1].columns
这个属性是一个ColDefs(列定义)对象。如果我们使用交互式提示符中的ColDefs.info()方法:
>>> 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']
start:
['', '', '', '']
dim:
['', '', '', '']
coord_type:
['', '', '', '']
coord_unit:
['', '', '', '']
coord_ref_point:
['', '', '', '']
coord_ref_value:
['', '', '', '']
coord_inc:
['', '', '', '']
time_ref_pos:
['', '', '', '']
它将显示表中所有列的属性，比如它们的名称、格式、bscales、bzeros等等。一个类似的输出，将显示列名及其格式，可以打印从脚本:
>>> hdul[1].columns
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'
)
修改columns中某个keyword的值
columns = fits.ColDefs(hdul[1].columns)
clumns1.change_attrib(col_name='DATA',attrib='format',new_value='8388608B')
我们也可以分别得到这些列的名称;例如:
>>> cols.names
['c1', 'c2', 'c3', 'c4']
返回字段名称的(Python)列表。
因为每个字段都是一个numpy对象，所以我们将拥有整个numpy工具库来使用。我们可以重新分配(更新)值:
>>> data['c4'][:] = 0
取一列的平均值:
>>> data['c3'].mean()
5.19999989271164
保存对文件的修改
如前所述，在用户打开文件，对头文件或数据做一些更改后，用户可以使用HDUList.writeto()保存更改。这将获取内存中的头文件和数据，并将它们写入磁盘上的一个新的FITS文件。后续操作可以执行的数据在内存和写入到另一个不同的文件，所有不重新使用原始数据到(更多)内存:
hdul.writeto('newtable.fits')
将当前hdulist的内容写入一个新的磁盘文件newfile.fits。如果以更新模式打开文件，还可以使用**HDUList.flush()方法将自open()**以来所做的所有更改写回原始文件。**close()**方法将对以更新模式打开的FITS文件执行相同的操作:
with fits.open('original.fits', mode='update') as hdul:
# Change something in hdul.
hdul.flush() # changes are written back to original.fits
# closing the file will also flush any changes and prevent further writing
物理常量：astropy.constants
astropy.constants包含许多在天文学中有用的物理常量，其中包括：
import astropy.constants as const
const.M_earth ## Earth mass
const.Ryd ## Rydberg constant
const.hbar ## Reduced Planck constant
const.mu0 ## Vacuum magnetic permeability
const.M_jup ## Jupiter mass
const.a0 ## Bohr radius
const.muB ## Bohr magneton
const.G ## Gravitational constant
const.M_sun ## Solar mass
const.pc ## Parsec , 1pc = 3.26 light year
const.GM_earth ## Nominal Earth mass parameter
const.N_A ## Avogadro's number
const.atm ## Standard atmosphere
const.k_B ## Boltzmann constant
const.GM_jup ## Nominal Jupiter mass parameter
const.R ## Gas constant
const.au ## Astronomical Unit
const.e ## Electron charge
const.kpc ## Kiloparsec
const.sigma_T ## Thomson scattering cross-section
const.GM_sun ## Nominal solar mass parameter
const.R_earth ## Nominal Earth equatorial radius
const.b_wien ## Wien wavelength displacement law constant
const.eps0 ## Vacuum electric permittivity
const.m_e ## Electron mass
const.sigma_sb ## Stefan-Boltzmann constant
const.L_bol0 ## Luminosity for absolute bolometric magnitude 0
const.R_jup ## Nominal Jupiter equatorial radius
const.c ## Speed of light in vacuum
const.g0 ## Standard acceleration of gravity
const.m_n ## Neutron mass
const.u ## Atomic mass
const.L_sun ## Nominal solar luminosity
const.R_sun ## Nominal solar radius
const.h ## Planck constant
const.m_p ## Proton mass