重新绑定时保留FITS文件的WCS信息

Question

Aim：重新绑定现有图像（FITS文件），并将新条目写入新的重新绑定图像（也是FITS文件）

问题：重新固定的FITS文件和原始FITS文件似乎坐标不匹配（问题后面显示的图）

过程：我将简要描述我的过程，以提供更多信息。第一步是读取现有的fits文件并定义numpy数组

from math import *
import numpy as np
import matplotlib.pyplot as plt 
from astropy.visualization import astropy_mpl_style
from astropy.io import fits 
import matplotlib.pyplot as plt
%matplotlib notebook
import aplpy
from aplpy import FITSFigure



file = 'F0621_HA_POL_0700471_HAWDHWPD_PMP_070-199Augy20.fits'
hawc = fits.open(file)

stokes_i = np.array(hawc[0].data)
stokes_i_rebinned = congrid(stokes_i,newdim,method="neighbour", centre=False, minusone=False)

这里的“conglid”是一个函数，我用它在近一小时内重新绑定，将原始数组重新绑定到“newdim”给定的新维度。现在的目标是将这个重新绑定的数组作为新文件写入FITS文件格式。我还有几个这样的数组，但为了简单起见，我只举了一个数组作为例子。为了保持头信息相同，我从现有FITS文件中读取该数组的头信息，并使用该信息将新数组写入新的FITS文件。写入后，可以像原始文件一样读取重新绑定的文件：-

header_0= hawc[0].header
fits.writeto("CasA_HAWC+_rebinned_congrid.fits", rebinned_stokes_i, header_0, overwrite=True)

rebinned_file = 'CasA_HAWC+_rebinned_congrid.fits'
hawc_rebinned= fits.open(rebinned_file)

为了检查重新固定的图像现在看起来如何，我绘制了它们

cmap = 'rainbow'

stokes_i = hawc[0]
stokes_i_rebinned = hawc_rebinned[0]

axi = FITSFigure(stokes_i, subplot=(1,2,1))  # generate FITSFigure as subplot to have two axes together
axi.show_colorscale(cmap=cmap)              # show I


axi_rebinned = FITSFigure(stokes_i_rebinned, subplot=(1,2,2),figure=plt.gcf())
axi_rebinned.show_colorscale(cmap=cmap)              # show I rebinned

# FORMATTING
axi.set_title('Stokes I (146 x 146)')
axi_rebinned.set_title('Rebinned Stokes I (50 x 50)')
axi_rebinned.axis_labels.set_yposition('right')
axi_rebinned.tick_labels.set_yposition('right')
axi.tick_labels.set_font(size='small')
axi.axis_labels.set_font(size='small')
axi_rebinned.tick_labels.set_font(size='small')
axi_rebinned.axis_labels.set_font(size='small')

正如您看到的原始图像和重新固定的图像，X、Y坐标似乎不匹配，我的最佳猜测是原始FITS文件的WCS（世界坐标系）没有正确复制到新FITS文件，从而导致任何不匹配。那么我该如何调整这些坐标呢

任何帮助都将不胜感激！谢谢

Answer 1

正如@astrochun已经说过的，您的重新装箱功能不会调整重新装箱图像的WCS。除了reproject和Montage之外，astropy.wcs.WCS对象还有slice()方法。您可以尝试使用它“重新存储”WCS，如下所示：

from astropy.wcs import WCS
import numpy as np
wcs = WCS(hawc[0].header, hawc)
wcs_rebinned = wcs.slice((np.s_[::2], np.s_[::2]))
wcs_hdr = wcs_rebinned.to_header()
header_0.update(wcs_hdr)  # but watch out for CD->PC conversion

您还应该在header_0= hawc[0].header中创建hawc[0].header的“真实”副本，例如作为header_0= hawc[0].header.copy()或者header_0.update(wcs_hdr)也将修改hawc[0].header

Answer 2

我将我的答案发布在一个非常松散的频道，如果这对其他人有用的话

congrid将不起作用，因为它不包含有关WCS的信息。例如，CD矩阵绑定到原始图像，而不是重新组合的集。有多种方法可以使用适当的WCS重新存储数据。您可能会考虑reproject，尽管这常常需要另一个WCS头重新返回到/

Montage（虽然不是Python工具，但有Python包装器）可能是另一种方式