使用 astropy 从旧数据写入新的 FITS 文件
Writing a new FITS file from old data with astropy
我要对 FITS 文件执行一个非常简单的操作（数据是 numpy 数组格式），但我无法使用 astropy 将其另存为新文件或覆盖现有文件。
我是 re-writing 一些使用 numpy pyfits 模块来处理天文 FITS 文件的旧代码 - 我想更新它以使用 astropy.io 适合模块。具体来说，我处理的一些数据是 3D 的，一些是 4D 的。 4D 的东西只是一个约定——第 4 轴不包含任何有用的信息（可以在此处找到数据示例：http://www.mpia.de/THINGS/Data_files/NGC_628_NA_CUBE_THINGS.FITS）。所以我更愿意删除多余的轴，然后我的其余代码可以继续进行而无需任何特殊要求。
这是我使用的旧 pyfits-based 代码，它运行良好：
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 所确认的，"image" 数组已正确重塑。但我终其一生都无法弄清楚如何将其保存到新的（或旧的）FITS 文件中。使用旧语法：
fitsfile.writeto('NGC628_2.fits',image,header)
... 给出错误消息，“AttributeError: 'numpy.ndarray' object 没有属性 'lower'。
如果相反，根据 astropy 文档，我只是省略图像和 header 并尝试保存到原始文件：
fitsfile.writeto('NGC628.fits')
然后我得到一个文件已经存在的错误。如果我提供关键字 "overwrite=True"，那么它会抱怨 "WinError 32 : The process cannot access the file because it is being used by another process : NGCC628.fits"。该文件绝对不能在任何其他程序中打开。
如果我指定新文件名 NGC628_2.fits，那么 Python 会崩溃（将我送回命令提示符）并且没有错误。写入了一个非常小的文件，其中仅包含 header 数据和图像数据的 none。编辑：如果我使用正确的命令使用图像和 header 数据编写一个新的 FITS 文件，则会发生完全相同的事情，例如pyfits.writeto('NGC628_2.fits',图像,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.
好吧，我想你很早就犯了一个小错误，你一直忽略了这个小错误，这让你陷入了寻找不同解决方案的困境。除了一件小而微妙的事情外，您其他无效的尝试也可能会奏效。我会先回顾一下哪里出了问题，然后在后面的一些案例中解释出了什么问题，所以这篇有点长...
首先，在您的原始代码中：
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)
这不是同一件事。在第一种情况下，它是 high-level writeto() "convenience function"。但现在您正在呼叫 fitsfile.writeto。这里 fitsfile 不是 astropy.io.fits 模块，而是文件 object 本身：
>>> type(fitsfile)
<class 'astropy.io.fits.hdu.hdulist.HDUList'>
HDUList class 也有方法 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 不同，它不接受数据或 header 参数，因为 HDUList 已经是一个 collection HDU，它已经具有关联数据和 headers。事实上，如果您只想对现有的 FITS 文件进行一些更改并将这些更改写入新文件，我不会像您最初那样使用 high-level 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' 更类似于在 write-mode 中打开纯文本文件，如 open('foo.txt', 'w')：它旨在修改磁盘上的现有文件 in-place。当文件关闭时，您所做的任何修改都会刷新到磁盘，并且无需使用 writeto() 或类似的东西。
您写道：
Using the old syntax :
fitsfile.writeto('NGC628_2.fits',image,header)
但正如我之前所写，这里没有 "old syntax"，您只是想使用 HDUList.writeto() 方法而不是 writeto() 函数 .
If I provide the keyword "overwrite=True", then it complains about "WinError 32 : The process cannot access the file because it is being used by another process : NGCC628.fits". The file is definitely NOT open in any other programs.
我看到你在 Windows——这是 Windows 的特殊限制 通常 ：Windows 没有如果该文件已经有任何打开的句柄，则允许写入、移动或删除该文件。我认为此处的错误消息具有误导性（此消息直接来自 Windows 本身）："it is being used by another process"。它应该类似于 "it is being used by this process or another process"。当你这样做时：
fitsfile = pyfits.open('NGC628.fits', mode='update')
您现在有一个打开的文件句柄，因此 Windows 不会让您在不先关闭该文件的情况下覆盖它（例如 fitsfile.close()）。
If I specify the new filename NGC628_2.fits, then Python crashes (sends me back to the command prompt) with no errors. A very small file is written that contains only the header data and none of the image data.
现在 这 听起来像是一个真正的错误。听起来您的 Python 解释器出现了段错误。我无法解释。但是我尝试按照您描述的相同步骤顺序重现它（在 Windows 上），但我做不到。最后的 fitsfile.writeto('NGC628_2.fits') 没有问题。我可以想象的一件事是，当您使用 mode='update' 打开文件时，内部状态管理会变得相当复杂，因为它必须跟踪您对 FITS 数据结构所做的所有更改，以便它可以正确地移动内容在现有文件中的磁盘上。有可能在尝试一些轻度病态操作的过程中（比如在更新过程中尝试覆盖文件）某些东西进入了未定义的状态。不过我不知道如何做同样的事情。
"correct" 使用 mode='update' 修改文件 in-place 可能类似于：
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 statement：当你使用 with 时，你会执行任何需要你在 with 语句的 body 中打开文件的操作，并且然后当 with 语句退出时，它将确保您所做的所有更改（如果有）都刷新到磁盘，并且文件已关闭。即使只是在 read-only 模式下打开文件，我也会使用它。这 在 Windows 上尤为重要，因为如您所见，Windows 对您是否已正确关闭文件特别挑剔。
同样，您的最后一个示例可以重写：
with pyfits.open('NGC628.fits', mode='update') as fitsfile:
image = fitsfile[0].data
header = fitsfile[0].header
image[:,:,:,:] = 5.0
最后，您写道：
I gather from the astropy documentation that simply opening the file in update mode and closing it should be enough, and in this case it is. But this same trick does not work when reshaping the image - the FITS file is unaltered.
但我不确定你在这里尝试了什么；你没有表现出那种尝试。如果我不得不猜测你写了类似的东西：
image = np.reshape(image, (z, y, x))
但这只是用新数组替换变量指向的值image。它没有更新原始 fitsfile[0].data。前面的 image[:] = 0.5 示例有效，因为此处 image 指向与 fitsfile[0].data 相同的数组，并且您正在修改其内容 in-place。但是 np.reshape(...) 不是 in-place 操作；它 returns 一个新数组。请参阅我之前的示例以了解正确的方法。
我要对 FITS 文件执行一个非常简单的操作（数据是 numpy 数组格式），但我无法使用 astropy 将其另存为新文件或覆盖现有文件。
我是 re-writing 一些使用 numpy pyfits 模块来处理天文 FITS 文件的旧代码 - 我想更新它以使用 astropy.io 适合模块。具体来说，我处理的一些数据是 3D 的，一些是 4D 的。 4D 的东西只是一个约定——第 4 轴不包含任何有用的信息（可以在此处找到数据示例：http://www.mpia.de/THINGS/Data_files/NGC_628_NA_CUBE_THINGS.FITS）。所以我更愿意删除多余的轴，然后我的其余代码可以继续进行而无需任何特殊要求。
这是我使用的旧 pyfits-based 代码，它运行良好：
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 所确认的，"image" 数组已正确重塑。但我终其一生都无法弄清楚如何将其保存到新的（或旧的）FITS 文件中。使用旧语法：
fitsfile.writeto('NGC628_2.fits',image,header)
... 给出错误消息，“AttributeError: 'numpy.ndarray' object 没有属性 'lower'。 如果相反，根据 astropy 文档，我只是省略图像和 header 并尝试保存到原始文件：
fitsfile.writeto('NGC628.fits')
然后我得到一个文件已经存在的错误。如果我提供关键字 "overwrite=True"，那么它会抱怨 "WinError 32 : The process cannot access the file because it is being used by another process : NGCC628.fits"。该文件绝对不能在任何其他程序中打开。
如果我指定新文件名 NGC628_2.fits，那么 Python 会崩溃（将我送回命令提示符）并且没有错误。写入了一个非常小的文件，其中仅包含 header 数据和图像数据的 none。编辑：如果我使用正确的命令使用图像和 header 数据编写一个新的 FITS 文件，则会发生完全相同的事情，例如pyfits.writeto('NGC628_2.fits',图像,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.
好吧，我想你很早就犯了一个小错误，你一直忽略了这个小错误，这让你陷入了寻找不同解决方案的困境。除了一件小而微妙的事情外，您其他无效的尝试也可能会奏效。我会先回顾一下哪里出了问题，然后在后面的一些案例中解释出了什么问题，所以这篇有点长...
首先，在您的原始代码中：
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)
这不是同一件事。在第一种情况下，它是 high-level writeto() "convenience function"。但现在您正在呼叫 fitsfile.writeto。这里 fitsfile 不是 astropy.io.fits 模块，而是文件 object 本身：
>>> type(fitsfile)
<class 'astropy.io.fits.hdu.hdulist.HDUList'>
HDUList class 也有方法 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 不同，它不接受数据或 header 参数，因为 HDUList 已经是一个 collection HDU，它已经具有关联数据和 headers。事实上，如果您只想对现有的 FITS 文件进行一些更改并将这些更改写入新文件，我不会像您最初那样使用 high-level 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' 更类似于在 write-mode 中打开纯文本文件，如 open('foo.txt', 'w')：它旨在修改磁盘上的现有文件 in-place。当文件关闭时，您所做的任何修改都会刷新到磁盘，并且无需使用 writeto() 或类似的东西。
您写道：
Using the old syntax :
fitsfile.writeto('NGC628_2.fits',image,header)
但正如我之前所写，这里没有 "old syntax"，您只是想使用 HDUList.writeto() 方法而不是 writeto() 函数 .
If I provide the keyword "overwrite=True", then it complains about "WinError 32 : The process cannot access the file because it is being used by another process : NGCC628.fits". The file is definitely NOT open in any other programs.
我看到你在 Windows——这是 Windows 的特殊限制 通常 ：Windows 没有如果该文件已经有任何打开的句柄，则允许写入、移动或删除该文件。我认为此处的错误消息具有误导性（此消息直接来自 Windows 本身）："it is being used by another process"。它应该类似于 "it is being used by this process or another process"。当你这样做时：
fitsfile = pyfits.open('NGC628.fits', mode='update')
您现在有一个打开的文件句柄，因此 Windows 不会让您在不先关闭该文件的情况下覆盖它（例如 fitsfile.close()）。
If I specify the new filename NGC628_2.fits, then Python crashes (sends me back to the command prompt) with no errors. A very small file is written that contains only the header data and none of the image data.
现在 这 听起来像是一个真正的错误。听起来您的 Python 解释器出现了段错误。我无法解释。但是我尝试按照您描述的相同步骤顺序重现它（在 Windows 上），但我做不到。最后的 fitsfile.writeto('NGC628_2.fits') 没有问题。我可以想象的一件事是，当您使用 mode='update' 打开文件时，内部状态管理会变得相当复杂，因为它必须跟踪您对 FITS 数据结构所做的所有更改，以便它可以正确地移动内容在现有文件中的磁盘上。有可能在尝试一些轻度病态操作的过程中（比如在更新过程中尝试覆盖文件）某些东西进入了未定义的状态。不过我不知道如何做同样的事情。
"correct" 使用 mode='update' 修改文件 in-place 可能类似于：
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 statement：当你使用 with 时，你会执行任何需要你在 with 语句的 body 中打开文件的操作，并且然后当 with 语句退出时，它将确保您所做的所有更改（如果有）都刷新到磁盘，并且文件已关闭。即使只是在 read-only 模式下打开文件，我也会使用它。这 在 Windows 上尤为重要，因为如您所见，Windows 对您是否已正确关闭文件特别挑剔。
同样，您的最后一个示例可以重写：
with pyfits.open('NGC628.fits', mode='update') as fitsfile:
image = fitsfile[0].data
header = fitsfile[0].header
image[:,:,:,:] = 5.0
最后，您写道：
I gather from the astropy documentation that simply opening the file in update mode and closing it should be enough, and in this case it is. But this same trick does not work when reshaping the image - the FITS file is unaltered.
但我不确定你在这里尝试了什么；你没有表现出那种尝试。如果我不得不猜测你写了类似的东西：
image = np.reshape(image, (z, y, x))
但这只是用新数组替换变量指向的值image。它没有更新原始 fitsfile[0].data。前面的 image[:] = 0.5 示例有效，因为此处 image 指向与 fitsfile[0].data 相同的数组，并且您正在修改其内容 in-place。但是 np.reshape(...) 不是 in-place 操作；它 returns 一个新数组。请参阅我之前的示例以了解正确的方法。