从头开始创建非常大的FITS文件

这个例子演示了如何使用 astropy.io.fits .

作者：Erik Bray

许可证：BSD

通常要创建一个图像拟合文件，可以执行以下操作：

import os
import numpy as np
from astropy.io import fits
data = np.zeros((40000, 40000), dtype=np.float64)
hdu = fits.PrimaryHDU(data=data)

然后使用 astropy.io.fits.writeto() 方法将新文件写入磁盘

hdu.writeto('large.fits')

然而，一个40000×40000的双倍数组几乎是12G！大多数系统无法在内存中创建它，仅仅是为了写到磁盘上。为了高效地创建如此大的文件，需要一些额外的工作和一些假设。

首先，可以预测标题中的标题有多大（如，有多少关键字）。FITS头必须用2880字节的块编写，足够每个块36个关键字（包括最后一个块中的END关键字）。典型的头有1到4个块之间，尽管有时更多。

因为我们写入FITS文件的第一件事就是头文件，所以我们需要写足够多的头块，这样就有足够的填充空间来添加新的关键字，而不必调整整个文件的大小。假设您希望标头默认使用4个块。然后，排除Astropy将自动添加的结束卡，创建标题并将其填充到36*4张卡。

创建一个存根数组来初始化HDU；它的确切大小无关紧要，只要它具有所需的维数

data = np.zeros((100, 100), dtype=np.float64)
hdu = fits.PrimaryHDU(data=data)
header = hdu.header
while len(header) < (36 * 4 - 1):
    header.append()  # Adds a blank card to the end

现在将NAXISn关键字调整为所需的数组大小，并只将头写到文件中。使用 hdu.writeto() 方法将使astropy“有用地”重置NAXISn关键字以匹配虚拟数组的大小。这是因为它努力确保只写入有效的FITS文件。相反，我们可以使用 astropy.io.fits.Header.tofile 方法：

header['NAXIS1'] = 40000
header['NAXIS2'] = 40000
header.tofile('large.fits')

最后，扩展文件的末尾以匹配数据的长度（加上头的长度）。通过查找文件的末尾并编写单个字节，这样可以在大多数系统上非常有效地完成这一点：

with open('large.fits', 'rb+') as fobj:
    # Seek past the length of the header, plus the length of the
    # Data we want to write.
    # 8 is the number of bytes per value, i.e. abs(header['BITPIX'])/8
    # (this example is assuming a 64-bit float)
    # The -1 is to account for the final byte that we are about to
    # write:
    fobj.seek(len(header.tostring()) + (40000 * 40000 * 8) - 1)
    fobj.write(b'\0')

一般来说，可以这样写：

shape = tuple(header[f'NAXIS{ii}'] for ii in range(1, header['NAXIS']+1))
with open('large.fits', 'rb+') as fobj:
    fobj.seek(len(header.tostring()) + (np.product(shape) * np.abs(header['BITPIX']//8)) - 1)
    fobj.write(b'\0')

在现代操作系统上，这将导致文件（超过头）被0填充到大约12GB的空间，以容纳40000 x 40000个映像。在支持稀疏文件创建的文件系统上（大多数Linux文件系统，但不是大多数mac使用的HFS+文件系统），这是一个非常快速、高效的操作。在其他系统中，您的里程数可能会有所不同。

这不是建立大型文件的唯一方法，但可能是最安全的方法之一。此方法还可以用于创建大型多扩展名FITS文件，只需稍加小心。

最后，我们将删除创建的文件：

os.remove('large.fits')