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 = _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 = _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的每个元素都是⼀个带有.heade r和.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'
>>> hdrments['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 file2/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
EXTEND
, ,PCOUNT, GCOUNT
NAXIS1, GCOUNT, PCOUNT, GROUPS
,TFIELDS, TFORM, TBCOL
还有许多其他保留的关键字,例如⽤于数据缩放的关键字,或⽤于表的列属性的关键字,如FITS标准中所述。其中⼤多数都可以通过Column或HDU对象的属性访问,例如⽤HDU .name设置EXTVER的 EXTNAME或HDU.ver。结构关键字将作为通⽤操作的结果进⾏检查和/或更新。例如,当:
1. 设置数据。关键字NAXIS*从数据形状(.data.shape)中设置,BITPIX从数据类型(.data.dtype)中设置。
2. 设置页眉。它的关键字根据数据属性进⾏更新(如上所述)。
3. 写⼀个⽂件。所有必要的关键字都被删除、更新或添加到标题中。
4. 调⽤HDU的验证⽅法(例如PrimaryHDU.verify())。⼀些关键字可以⾃动修复。
在这种情况下,⽤户可能为这些关键字分配的任何值都将被覆盖。
处理数据 DATA
fits中的数据⼀般分为两种,图像数据(Image Data)和表格数据(Table Data),对于脉冲星观测产⽣的的数据为表格数据,所以在此并不包括对图像数据的内容。
本节描述直接使⽤FITS包以FITS格式读取和写⼊表数据。然⽽,在某些情况下,⾼级的统⼀⽂件通常就⾜够了,⽽且在某种程度上使⽤起来更⽅便。有关详细信息,请参阅 部分。
FITS表扩展的数据部分也在.data属性中:
>>> fits_table_filename = _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对象。
python怎么读取py文件
像标题关键字,⼀个列可以通过索引,如上,或通过名称:
>>> 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等等。⼀个类似的输出,将显⽰列名及其格式,可以打印从脚本:

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。