使用 astropy 从旧数据写入新的 FITS 文件

Question

我对 FITS 文件执行了一个非常简单的操作（数据是 numpy 数组格式），但我无法将其保存为新文件或覆盖现有文件。

我正在重写一些使用 numpy pyfits 模块处理天文 FITS 文件的旧代码 - 我想更新它以使用 astropy.io fits 模块。具体来说，我使用的一些数据是 3D 的，有些是 4D 的。4D 的东西只是一个约定 - 第 4 轴不包含有用的信息（可以在此处找到数据示例：http : //www.mpia.de/THINGS/Data_files/NGC_628_NA_CUBE_THINGS.FITS）。所以我更喜欢删除额外的轴，然后我的其余代码可以在没有任何特殊要求的情况下继续进行。

这是我使用的基于 pyfits 的旧代码，效果很好：

import numpy
import pyfits

filename = 'NGC628.fits'
outfile  = 'NGC628_reshaped.fits'

# Get the shape of the file
fitsfile=pyfits.open(filename)
image = fitsfile[0].data
header =fitsfile[0].header
z = image.shape[1]  # No. channels                  
y = image.shape[2]  # No. x pixels
x = image.shape[3]  # No. y pixels

newimage = numpy.reshape(image,[z,y,x])

pyfits.core.writeto(outfile,newimage,header, clobber=True)

没有什么复杂的，它只是重塑一个数组并将其保存到一个新文件中。奇妙。现在我想用 astropy fits 模块替换它。如果我做 ：

import numpy
from astropy.io import fits as pyfits

fitsfile = pyfits.open('NGC628.fits', mode='update')
image  = fitsfile[0].data
header = fitsfile[0].header
z = image.shape[1]                  
y = image.shape[2]
x = image.shape[3]
image = numpy.reshape(image,[z,y,x])

......那么到目前为止，很好。正如 image.shape 所确认的那样，“图像”数组被正确地重塑。但是我一生都无法弄清楚如何将其保存到新的（或旧的）FITS 文件中。使用旧语法：

fitsfile.writeto('NGC628_2.fits',image,header)

...给出错误消息，“AttributeError: 'numpy.ndarray' object has no attribute 'lower'。如果相反，根据 astropy 文档，我只是省略图像和标题并尝试保存到原始文件：

fitsfile.writeto('NGC628.fits')

然后我收到一个错误，该文件已经存在。如果我提供关键字“overwrite=True”，那么它会抱怨“WinError 32：该进程无法访问该文件，因为它正被另一个进程使用：NGCC628.fits”。该文件绝对不会在任何其他程序中打开。

如果我指定新文件名 NGC628_2.fits，则 Python 会崩溃（将我送回命令提示符）而不会出现错误。写入一个非常小的文件，其中仅包含标题数据，不包含任何图像数据。编辑：如果我使用正确的命令使用图像和标题数据编写新的 FITS 文件，则会发生完全相同的事情，例如 pyfits.writeto('NGC628_2.fits',image,header)。

只是为了让事情更加混乱，如果我做一个稍微简单的操作，比如说，将所有图像数据设置为一个常量值，然后关闭文件：

import numpy
from astropy.io import fits as pyfits
fitsfile = pyfits.open('NGC628.fits', mode='update')
image  = fitsfile[0].data
header = fitsfile[0].header
image[:,:,:,:] = 5.0
fitsfile.close()

然后这有效 - 原始文件现在是一个数组，其中每个值都等于 5。我从 astropy 文档中收集到，只需在更新模式下打开文件并关闭它就足够了，在这种情况下就是这样。但这同样的伎俩并不能重塑形象，当工作- FITS文件是没有改变的。

那我到底做错了什么？更新原始文件或保存到新文件都可以（最好是后者），但我无法使任一操作正常工作。

编辑：我有 Python 版本 3.5.3、numpy 版本 1.17.3、astropy 版本 3.2.3，我运行的是 Windows 10。

Answer 1

好吧，我认为你很早就犯了一个小错误，你一直只是忽视，这让你陷入了疯狂的追逐，寻找不同的解决方案。除了一个小而微妙的事情之外，您的其他无效尝试也可能会奏效。我先把哪里出错了，然后在后面的一些案例中说明哪里出错了，所以这有点长......

首先，在您的原始代码中，您有：

pyfits.core.writeto(outfile,newimage,header, clobber=True)

这可以写得更简单：

pyfits.writeto(outfile, newimage, header, clobber=True)

该pyfits.core模块是一个实现细节，几乎没有理由直接引用它。你可以直接调用这个函数作为pyfits.writeto().

如果你有这样的，那么你的现有代码会工作正是通过改变进口为前

from astropy.io import fits as pyfits

唯一需要注意的是，该clobber参数已重命名overwrite，尽管我认为clobber它仍然有效，但带有弃用警告。

第一个错误

您将其更改为但似乎忽略了的是：

fitsfile.writeto('NGC628_2.fits',image,header)

这不是一回事。在第一种情况下，它是高级别的writeto()“便利功能”。但是现在您正在调用fitsfile.writeto. 这里fitsfile不是astropy.io.fits模块，而是文件对象本身：

>>> type(fitsfile)
<class 'astropy.io.fits.hdu.hdulist.HDUList'>

的HDUList类也有方法HDUList.writeto，其执行类似的功能，但它需要不同的参数。您可以通过多种方式进行检查，例如：

>>> help(fitsfile.writeto)
Help on method writeto in module astropy.io.fits.hdu.hdulist:

writeto(fileobj, output_verify='exception', overwrite=False, checksum=False) method of astropy.io.fits.hdu.hdulist.HDUList instance
    Write the `HDUList` to a new file.

它唯一需要的参数是将文件写入的文件名。与另一个不同，writeto它不接受数据或标头参数，因为HDUList它已经是一个已经具有关联数据和标头的 HDU 的集合。事实上，如果您只想对现有的 FITS 文件进行一些更改并将这些更改写到一个新文件中，我不会writeto()像您最初那样使用高级别的- 如果您正在创建一个数组，它最有用从头开始，只想快速将其写入新的 FITS 文件。所以你的原始代码也可以这样写：

import numpy as np
import astropy.io.fits as pyfits

filename = 'NGC628.fits'
outfile  = 'NGC628_reshaped.fits'

# Get the shape of the file
fitsfile = pyfits.open(filename)
image = fitsfile[0].data
z = image.shape[1]  # No. channels                  
y = image.shape[2]  # No. x pixels
x = image.shape[3]  # No. y pixels

# Replace the primary HDU's data in-place; this just manipulates the
# in-memory HDU data structure; it does not change anything on disk
fitsfile[0].data = np.reshape(image, [z,y,x])

# Write the new HDU structure to outfile
fitsfile.writeto(outfile, overwrite=True)

更新模式

接下来，您尝试通过使用 来打开文件来修改文件mode='update'，但这实际上并不是要使用更新模式的方式。 mode='update'更类似于在写入模式下打开纯文本文件，例如open('foo.txt', 'w')：它旨在就地修改磁盘上的现有文件。当文件关闭时，您所做的任何修改都会刷新到磁盘，并且不需要使用writeto()或类似的东西。

你写了：

使用旧语法：

fitsfile.writeto('NGC628_2.fits',image,header)

但是正如我之前写的，这里没有“旧语法”，您只是尝试使用HDUList.writeto()method 而不是writeto() function。

如果我提供关键字“overwrite=True”，那么它会抱怨“WinError 32：该进程无法访问该文件，因为它正被另一个进程使用：NGCC628.fits”。该文件绝对不会在任何其他程序中打开。

我看你是在Windows上-这是Windows的一个特定的限制，一般的：Windows不允许写入，移动或删除一个文件，如果已经有任何打开的句柄到该文件。我认为这里的错误消息具有误导性（此消息直接来自 Windows 本身）：“它正在被另一个进程使用”。它应该类似于“它正在被这个进程或另一个进程使用”。当你这样做时：

fitsfile = pyfits.open('NGC628.fits', mode='update')

你现在有一个打开的文件句柄，所以 Windows 不会让你在没有先关闭它的情况下覆盖它（例如fitsfile.close()）。

如果我指定新文件名 NGC628_2.fits，则 Python 会崩溃（将我送回命令提示符）而不会出现错误。写入一个非常小的文件，其中仅包含标题数据，不包含任何图像数据。

现在这听起来像是一个真正的错误。听起来您的 Python 解释器出现了段错误。我无法解释。但是我尝试按照您描述的相同步骤顺序复制它（在 Windows 上），但我不能。最后的fitsfile.writeto('NGC628_2.fits')工作没有问题。我能想象的一件事是，当你打开一个文件时mode='update'，内部状态管理可能会变得非常复杂，因为它必须跟踪您更改的有关 FITS 数据结构的所有内容，以便它可以在现有文件中的磁盘上正确移动内容。有可能在尝试一些温和的病理操作（例如尝试在文件更新过程中覆盖文件）的过程中，某些东西进入了未定义的状态。我不知道如何做同样的事情。

mode='update'就地修改文件的“正确”用法可能类似于：

with pyfits.open(filename, mode='update') as fitsfile:
    image = fitsfile[0].data
    z, y, x = image.shape[1:]
    # Replace the original HDU data in-place
    fitsfile[0].data = image.reshape((z, y, x))

就是这样！我会习惯使用该with语句：当您with执行任何需要在with语句主体内打开文件的操作时，当with退出该语句时，它将确保您所做的所有更改（如果有的话） ) 刷新到磁盘，并且文件已关闭。即使只是在只读模式下打开文件，我也会使用它。这在 Windows 上尤其重要，因为正如您所发现的，Windows 对您是否正确关闭文件特别挑剔。

同样，您的最后一个示例可以重写：

with pyfits.open('NGC628.fits', mode='update') as fitsfile:
    image  = fitsfile[0].data
    header = fitsfile[0].header
    image[:,:,:,:] = 5.0

最后，你写道：

我从 astropy 文档中收集到，只需在更新模式下打开文件并关闭它就足够了，在这种情况下就是这样。但这同样的伎俩并不能重塑形象，当工作- FITS文件是没有改变的。

但我不确定你在这里尝试了什么；你没有表现出这种尝试。如果我不得不猜测你写的是这样的：

    image = np.reshape(image, (z, y, x))

但所有这些都是用一个新数组替换image变量指向的值。它不是更新原始fitsfile[0].data. 前面的示例image[:] = 0.5有效，因为这里image指向fitsfile[0].data与您相同的数组，并且您正在就地修改其内容。但np.reshape(...)不是就地操作；它返回一个新数组。请参阅我之前的示例以了解正确的方法。