不太熟悉的物体#

在本章中,我们将讨论使用频率较低的FITS数据结构。它们包括ASCII表、可变长度表和随机访问组FITS文件。

ASCII表格#

FITS标准支持二进制和ASCII表。在ASCII表中,所有数据都以人类可读的文本形式存储,因此解析文本以获取数字数据需要更多的空间和额外的处理。根据列的格式,浮点数据也可能会失去精度。

astropy ,ASCII表和二进制表的接口基本相同(即数据在 .data 属性和 field() 方法用于引用列并返回 numpy 阵列)。读表的时候, astropy 将自动检测它是什么样的表。

>>> from astropy.io import fits
>>> filename = fits.util.get_testdata_filepath('ascii.fits')
>>> hdul = fits.open(filename)
>>> hdul[1].data[:1]  
FITS_rec([(10.123, 37)],
         dtype=(numpy.record, {'names':['a','b'], 'formats':['S10','S5'], 'offsets':[0,11], 'itemsize':16}))
>>> hdul[1].data['a']
array([  10.123,    5.2  ,   15.61 ,    0.   ,  345.   ])
>>> hdul[1].data.formats
['E10.4', 'I5']
>>> hdul.close()

注意,记录数组中的格式引用的原始数据是ASCII字符串(因此是“a11”和“a5”),但是 .formats 属性保留了原始格式规范('E10.4'和'I5')。

创建ASCII表#

从头开始创建ASCII表与创建二进制表类似。区别在于列定义。ASCII表中的列/字段比二进制表中的列/字段更受限制。它不允许一个单元格中有多个数值。而且,它只支持二进制表中允许的子集,即字符串、整数和(单精度和双精度)浮点数。不允许使用布尔数和复数。

格式语法(TFORM关键字的值)与二进制表的格式语法不同。他们是:

Aw         Character string
Iw         (Decimal) Integer
Fw.d       Double precision real
Ew.d       Double precision real, in exponential notation
Dw.d       Double precision real, in exponential notation

其中w是宽度,d是小数点后的位数。ASCII表和二进制表之间的语法差异可能会令人困惑。例如,一个由3个字符组成的字段在二进制表中被指定为“3A”,在ASCII表中被指定为“A3”。

另一个区别是在使用 TableHDU.from_columns() 方法,还有那个 Column 应提供 ascii=True 为了不含糊而争论。

备注

尽管二进制表在大多数FITS文件中更常见,但FITS格式的早期版本只支持ASCII表。这就是上课的原因 TableHDU 专门用于表示ASCII表,而 BinTableHDU 更明确的是它表示一个二进制表。这些名字来自于 XTENSION 表标题中的关键字,它是 TABLE 对于ASCII表和 BINTABLE 对于二进制表。

TableHDU.from_columns() 可以这样使用:

>>> import numpy as np

>>> a1 = np.array(['abcd', 'def'])
>>> r1 = np.array([11., 12.])
>>> col1 = fits.Column(name='abc', format='A3', array=a1, ascii=True)
>>> col2 = fits.Column(name='def', format='E', array=r1, bscale=2.3,
...                    bzero=0.6, ascii=True)
>>> col3 = fits.Column(name='t1', format='I', array=[91, 92, 93], ascii=True)
>>> hdu = fits.TableHDU.from_columns([col1, col2, col3])
>>> hdu.data
FITS_rec([('abc', 11.0, 91), ('def', 12.0, 92), ('', 0.0, 93)],
         dtype=(numpy.record, [('abc', 'S3'), ('def', 'S15'), ('t1', 'S10')]))

应该注意的是,当列的格式明确地特定于ASCII表时,不必指定 ascii=TrueColDefs 建造师。在这种情况下 is 模棱两可因为格式代码 'I' 表示二进制表中的16位整数,而在ASCII表中,它在技术上不是有效的格式。ASCII表格格式代码在技术上要求每列有一个字符宽度,例如 'I10' 创建一列,该列可容纳10个字符宽的整数。

然而, astropy 允许在某些情况下省略宽度指定。当它被省略时 'I' 格式列使用精确表示列中所有整数所需的最小宽度。使用此快捷方式的唯一问题是它与二进制表的模糊性 'I' 格式,所以指定 ascii=True 是一个很好的实践(虽然 astropy 在大多数情况下,你仍然会明白你的意思)。

可变长度数组表#

FITS标准还支持可变长度数组表。其基本思想是,有时希望在同一字段(列)中具有相同数据类型但长度/维度不同的单元格的表。与标准的表数据结构相比,变长表在不同单元中具有较大的数据长度动态范围,可以节省存储空间。

可变长度数组表格可以具有一个或多个可变长度的字段(列)。同一个表中的其余字段(列)仍然可以是固定长度的常规字段。表中可变长度数组的数据不存储在主数据表中;它们存储在主数据表之后的补充数据区域堆中。 astropy 在阅读过程中将自动检测它是哪种类型的字段;用户不需要特殊操作。数据类型规范(即TForm关键字的值)使用额外的字母‘P’(或‘Q’),格式为:

rPt(max)

哪里 r 可以是0或1(通常省略,因为它不适用于可变长度数组), t 是基础数据类型(L、B、I、J等)的字母代码之一;目前中的可变长度数组字段不支持X格式 astropy ),以及 max 列中任何数组的最大元素数。因此,对于可变长度的int16字段,对应的格式规范例如是‘pj(100)’。存储在主数据表字段中的是数组描述符。这包括两个32位有符号整数值(在‘P’格式的情况下),或两个64位有符号整数值(在‘Q’格式的情况下):存储的数组的元素数(数组长度),后跟数组第一个元素的零索引字节偏移量,从堆区域的开始测量。

备注

虽然P格式使用32位带符号整数,但FITS标准没有定义负值的含义。P格式从字节0到字节0的索引 \(2^{31} - 1\) 。根据变量数组的格式(整型、浮点型或双精度型)以及行数的不同,可能需要使用Q格式来分配足够的堆空间。

例子#

此示例显示数据类型为int16::

>>> filename = fits.util.get_testdata_filepath('variable_length_table.fits')
>>> hdul = fits.open(filename)
>>> hdul[1].header['tform1']
'PI(3)'
>>> print(hdul[1].data.field(0))
[array([45, 56], dtype=int16) array([11, 12, 13], dtype=int16)]
>>> hdul.close()

在这个字段中,第一行有一个元素,第二行有两个元素,等等。访问可变长度字段与常规字段几乎相同,只是通常不可能同时对整个字段进行操作。用户必须像独立数组一样逐行处理字段。

创建可变长度数组表#

创建可变长度表与创建常规表几乎相同。唯一的区别在于创建可变长度数组的字段定义。首先,数据类型规范需要'P'字母,其次,字段数据必须是objects数组(如 numpy 模块)。

例子#

下面是一个创建包含两个字段的表的示例:一个是常规字段,另一个是可变长度数组:

>>> col1 = fits.Column(
...    name='var', format='PI()',
...    array=np.array([[45, 56], [11, 12, 13]], dtype=np.object_))
>>> col2 = fits.Column(name='xyz', format='2I', array=[[11, 3], [12, 4]])
>>> hdu = fits.BinTableHDU.from_columns([col1, col2])
>>> data = hdu.data
>>> data  
FITS_rec([([45, 56], [11,  3]), ([11, 12, 13], [12,  4])],
         dtype=(numpy.record, [('var', '<i4', (2,)), ('xyz', '<i2', (2,))]))
>>> hdu.writeto('variable_length_table.fits')
>>> with fits.open('variable_length_table.fits') as hdul:
...     print(repr(hdul[1].header))
XTENSION= 'BINTABLE'           / binary table extension
BITPIX  =                    8 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                   12 / length of dimension 1
NAXIS2  =                    2 / length of dimension 2
PCOUNT  =                   10 / number of group parameters
GCOUNT  =                    1 / number of groups
TFIELDS =                    2 / number of table fields
TTYPE1  = 'var     '
TFORM1  = 'PI(3)   '
TTYPE2  = 'xyz     '
TFORM2  = '2I      '

随机访问组#

FITS标准支持的另一个不太熟悉的数据结构是random access group。这个约定是在引入二进制表扩展之前建立的。在大多数情况下,它的使用现在可以被二进制表取代。它主要用于无线电干涉测量。

与主HDU一样,随机访问组HDU总是FITS文件的第一个HDU。它的数据有一个或多个组。每个组可以有任意数量(包括0)的参数,以及一个图像。参数和图像具有相同的数据类型。

同一个HDU中的所有组具有相同的数据结构,即相同的数据类型(由关键字BITPIX指定,如在image HDU中)、相同的参数数目(由PCOUNT指定)、相同的大小和形状(由NAXISn关键字指定)。组数由GCOUNT指定,关键字NAXIS1始终为0。因此,随机访问组HDU的总数据大小为:

|BITPIX| * GCOUNT * (PCOUNT + NAXIS2 * NAXIS3 * ... * NAXISn)

标题和摘要#

访问随机访问组HDU的头与任何其他HDU没有什么不同;可以使用.header属性。

HDU的内容同样可以通过使用 HDUList.info() 方法:

>>> filename = fits.util.get_testdata_filepath('group.fits')
>>> hdul = fits.open(filename)
>>> hdul[0].header['groups']
True
>>> hdul[0].header['gcount']
10
>>> hdul[0].header['pcount']
3
>>> hdul.info()
Filename: ...group.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 GroupsHDU       15   (5, 3, 1, 1)   float32   10 Groups  3 Parameters

数据:组参数#

与其他HDU一样,随机访问组HDU的数据部分位于 .data 属性。它包括参数和图像数组。

实例#

要显示第三组的内容,包括参数和数据:

>>> hdul[0].data[2]  
(2.0999999, 42.0, 42.0, array([[[[30., 31., 32., 33., 34.],
         [35., 36., 37., 38., 39.],
         [40., 41., 42., 43., 44.]]]], dtype=float32))

数据首先列出指定组的所有参数,然后列出图像数组。提醒一下,在Python或C约定中,此文件中的图像数据的形状为(1,1,1,4,3),或IRAF或Fortran惯例中的(3,4,1,1,1)。

要访问参数,首先要找出参数名是什么,使用 .parnames 属性:

>>> hdul[0].data.parnames # get the parameter names
['abc', 'xyz', 'xyz']

组参数可以通过 par() 方法。就像桌子一样 field() 参数可以是:或方法名:

>>> hdul[0].data.par(0)[8]  # Access group parameter by name or by index  
8.1
>>> hdul[0].data.par('abc')[8]  
8.1

请注意,参数名“xyz”出现两次。这是随机访问组中的一个特性,它意味着将值相加。因此:

>>> hdul[0].data.parnames  # get the parameter names
['abc', 'xyz', 'xyz']
>>> hdul[0].data.par(1)[8]  # Duplicate parameter name 'xyz'
42.0
>>> hdul[0].data.par(2)[8]
42.0
>>> # When accessed by name, it adds the values together if the name is
>>> # shared by more than one parameter
>>> hdul[0].data.par('xyz')[8]
84.0

这个 par() 是整个数据对象或一个数据项(组)的方法。因此,有两种可能的方法来获取某个组的组参数,这与表数据中的情况类似(使用 field() 方法:

>>> hdul[0].data.par(0)[8]  
8.1
>>> hdul[0].data[8].par(0)  
8.1

另一方面,要修改组参数,我们可以直接分配新值(如果最后访问行/组号),或者使用 setpar() 方法(如果先访问行/组号)。方法 setpar() 如果参数由多个参数共享,则还需要按名称更新:

>>> # Update group parameter when selecting the row (group) number last
>>> hdul[0].data.par(0)[8] = 99.
>>> # Update group parameter when selecting the row (group) number first
>>> hdul[0].data[8].setpar(0, 99.)  # or:
>>> hdul[0].data[8].setpar('abc', 99.)
>>> # Update group parameter by name when the name is shared by more than
>>> # one parameters, the new value must be a tuple of constants or
>>> # sequences
>>> hdul[0].data[8].setpar('xyz', (2445729., 0.3))
>>> hdul[0].data[8:].par('xyz')  
array([2.44572930e+06, 8.40000000e+01])

数据:图像数据#

可以访问数据部分的图像数组 data 数据对象的属性。A numpy 返回数组::

>>> print(hdul[0].data.data[8])  
[[[[120. 121. 122. 123. 124.]
   [125. 126. 127. 128. 129.]
   [130. 131. 132. 133. 134.]]]]
>>> hdul.close()

创建随机访问组HDU#

要从头开始创建随机访问组HDU,请使用 GroupData 将数据封装到组数据结构中,并使用 GroupsHDU 创建HDU本身。

例子#

要创建随机访问组HDU::

>>> # Create the image arrays. The first dimension is the number of groups.
>>> imdata = np.arange(150.0).reshape(10, 1, 1, 3, 5)
>>> # Next, create the group parameter data, we'll have two parameters.
>>> # Note that the size of each parameter's data is also the number of
>>> # groups.
>>> # A parameter's data can also be a numeric constant.
>>> pdata1 = np.arange(10) + 0.1
>>> pdata2 = 42
>>> # Create the group data object, put parameter names and parameter data
>>> # in lists assigned to their corresponding arguments.
>>> # If the data type (bitpix) is not specified, the data type of the
>>> # image will be used.
>>> x = fits.GroupData(imdata, bitpix=-32,
...                    parnames=['abc', 'xyz', 'xyz'],
...                    pardata=[pdata1, pdata2, pdata2])
>>> # Now, create the GroupsHDU and write to a FITS file.
>>> hdu = fits.GroupsHDU(x)
>>> hdu.writeto('test_group.fits')
>>> hdu.header
SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                  -32 / array data type
NAXIS   =                    5 / number of array dimensions
NAXIS1  =                    0
NAXIS2  =                    5
NAXIS3  =                    3
NAXIS4  =                    1
NAXIS5  =                    1
EXTEND  =                    T
GROUPS  =                    T / has groups
PCOUNT  =                    3 / number of parameters
GCOUNT  =                   10 / number of groups
PTYPE1  = 'abc     '
PTYPE2  = 'xyz     '
PTYPE3  = 'xyz     '
>>> data = hdu.data
>>> hdu.data  
GroupData([ (0.1       , 42., 42., [[[[  0.,   1.,   2.,   3.,   4.], [  5.,   6.,   7.,   8.,   9.], [ 10.,  11.,  12.,  13.,  14.]]]]),
           (1.10000002, 42., 42., [[[[ 15.,  16.,  17.,  18.,  19.], [ 20.,  21.,  22.,  23.,  24.], [ 25.,  26.,  27.,  28.,  29.]]]]),
           (2.0999999 , 42., 42., [[[[ 30.,  31.,  32.,  33.,  34.], [ 35.,  36.,  37.,  38.,  39.], [ 40.,  41.,  42.,  43.,  44.]]]]),
           (3.0999999 , 42., 42., [[[[ 45.,  46.,  47.,  48.,  49.], [ 50.,  51.,  52.,  53.,  54.], [ 55.,  56.,  57.,  58.,  59.]]]]),
           (4.0999999 , 42., 42., [[[[ 60.,  61.,  62.,  63.,  64.], [ 65.,  66.,  67.,  68.,  69.], [ 70.,  71.,  72.,  73.,  74.]]]]),
           (5.0999999 , 42., 42., [[[[ 75.,  76.,  77.,  78.,  79.], [ 80.,  81.,  82.,  83.,  84.], [ 85.,  86.,  87.,  88.,  89.]]]]),
           (6.0999999 , 42., 42., [[[[ 90.,  91.,  92.,  93.,  94.], [ 95.,  96.,  97.,  98.,  99.], [100., 101., 102., 103., 104.]]]]),
           (7.0999999 , 42., 42., [[[[105., 106., 107., 108., 109.], [110., 111., 112., 113., 114.], [115., 116., 117., 118., 119.]]]]),
           (8.10000038, 42., 42., [[[[120., 121., 122., 123., 124.], [125., 126., 127., 128., 129.], [130., 131., 132., 133., 134.]]]]),
           (9.10000038, 42., 42., [[[[135., 136., 137., 138., 139.], [140., 141., 142., 143., 144.], [145., 146., 147., 148., 149.]]]])],
           dtype=(numpy.record, [('abc', '<f4'), ('xyz', '<f4'), ('_xyz', '<f4'), ('DATA', '<f4', (1, 1, 3, 5))]))

压缩图像数据#

开发了一种在FITS二进制表中存储压缩图像数据的通用技术。此约定中使用的原则是首先将n维图像划分为子图像或“平铺”的矩形网格。然后将每个分片压缩为一个连续的数据块,得到的压缩字节流存储在FITS二进制表中可变长度列的一行中。支持压缩图像块的几种常用算法。其中包括Gzip、Rice、IRAF像素列表(PLIO)和Hcompress。

更多关于压缩的建议,请参考图片:

以及“Registered FITS Convention,Tiled Image Compression Convention”:

在中访问压缩图像数据 astropy 可选,使用 astropy.io.fits.compression 包含在C共享库中的模块(压缩.so). 如果在压缩模块不可用时试图访问包含压缩图像数据的HDU,则会通知用户该问题,并将HDU视为标准二进制表HDU。此通知仅在第一次遇到压缩图像数据时发出。这样,就不需要压缩模块 astropy 工作。

标题和摘要#

astropy ,在用户看来,压缩图像HDU的头部与任何图像头一样。存储在FITS文件中的实际头是二进制表HDU的头,其中有一组特殊的关键字,由约定定义,用来描述压缩图像的结构。二进制表HDU头和图像HDU头的转换都是在后台进行的。由于HDU实际上是一个二进制表,所以它可能不会在FITS文件中显示为主HDU。

例子#

可以使用访问HDU报头的内容 .header 属性:

>>> filename = fits.util.get_testdata_filepath('compressed_image.fits')

>>> hdul = fits.open(filename)
>>> hdul[1].header
XTENSION= 'IMAGE   '           / Image extension
BITPIX  =                   16 / data type of original image
NAXIS   =                    2 / dimension of original image
NAXIS1  =                   10 / length of original image axis
NAXIS2  =                   10 / length of original image axis
PCOUNT  =                    0 / number of parameters
GCOUNT  =                    1 / number of groups

相应的二进制表HDU的内容可以使用hidden访问 ._header 属性。但是,所有带有HDU头的用户界面都应该通过图像头( .header 属性):

>>> hdul[1]._header
XTENSION= 'BINTABLE'           / binary table extension
BITPIX  =                    8 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                    8 / width of table in bytes
NAXIS2  =                   10 / number of rows in table
PCOUNT  =                   60 / number of group parameters
GCOUNT  =                    1 / number of groups
TFIELDS =                    1 / number of fields in each row
TTYPE1  = 'COMPRESSED_DATA'    / label for field 1
TFORM1  = '1PB(6)  '           / data format of field: variable length array
ZIMAGE  =                    T / extension contains compressed image
ZTENSION= 'IMAGE   '           / Image extension
ZBITPIX =                   16 / data type of original image
ZNAXIS  =                    2 / dimension of original image
ZNAXIS1 =                   10 / length of original image axis
ZNAXIS2 =                   10 / length of original image axis
ZPCOUNT =                    0 / number of parameters
ZGCOUNT =                    1 / number of groups
ZTILE1  =                   10 / size of tiles to be compressed
ZTILE2  =                    1 / size of tiles to be compressed
ZCMPTYPE= 'RICE_1  '           / compression algorithm
ZNAME1  = 'BLOCKSIZE'          / compression block size
ZVAL1   =                   32 / pixels per block
ZNAME2  = 'BYTEPIX '           / bytes per pixel (1, 2, 4, or 8)
ZVAL2   =                    2 / bytes per pixel (1, 2, 4, or 8)
EXTNAME = 'COMPRESSED_IMAGE'   / name of this binary table extension

HDU的内容可以通过使用 info() 方便函数或方法:

>>> fits.info(filename)
Filename: ...compressed_image.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU       4   ()
  1  COMPRESSED_IMAGE    1 CompImageHDU      7   (10, 10)   int16

>>> hdul.info()
Filename: ...compressed_image.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU       4   ()
  1  COMPRESSED_IMAGE    1 CompImageHDU      7   (10, 10)   int16

数据#

与头部一样,压缩图像HDU的数据在用户看来是标准的未压缩图像数据。实际数据以二进制表数据的形式存储在FITS文件中,该表至少包含一列(压缩的数据)。此可变长度列的每一行都包含压缩相应图像块所生成的字节流。也可能出现几个可选列。这些包括未压缩的_数据,用于保存无法压缩的图像块的未压缩像素值,ZSCALE和ZZERO用于保存线性比例因子和零点偏移,这可能需要将原始未压缩值转换回原始图像像素值,和ZBLANK保存用于表示图像中未定义像素(如果有)的整数值。

例子#

未压缩的HDU数据的内容可以使用 .data 属性:

>>> hdul[1].data
array([[ 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]], dtype=int16)
>>> hdul.close()

压缩数据可以通过 .compressed_data 属性,但很少需要直接访问。它可能有助于直接复制压缩数据,而无需首先解压缩。

创建压缩图像HDU#

要从头开始创建压缩图像HDU,请构造 CompImageHDU 对象从未压缩的图像数据数组及其关联的图像标头。从那里,HDU可以像任何其他图像HDU一样被处理。

例子#

要创建压缩图像HDU::

>>> imageData = np.arange(100).astype('i2').reshape(10, 10)
>>> imageHeader = fits.Header()
>>> hdu = fits.CompImageHDU(imageData, imageHeader)
>>> hdu.writeto('compressed_image.fits')

的API文档 CompImageHDU 初始值设定项方法描述了构造 CompImageHDU 对象。