不太熟悉的物体#
在本章中,我们将讨论使用频率较低的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=True
在 ColDefs
建造师。在这种情况下 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
对象。