使用 astropy.fits 和 numpy 对 SWIFT 拟合图像应用巧合校正

Question

这个问题可能有点专业，但希望有人能提供帮助。我通常使用 IDL，但为了开发管道，我希望使用 python 来提高运行倍。

我的 fits 文件处理设置如下：

import numpy as numpy
from astropy.io import fits

#Directory: /Users/UCL_Astronomy/Documents/UCL/PHASG199/M33_UVOT_sum/UVOTIMSUM/M33_sum_epoch1_um2_norm.img
with fits.open('...') as ima_norm_um2:
    #Open UVOTIMSUM file once and close it after extracting the relevant values:
    ima_norm_um2_hdr  = ima_norm_um2[0].header
    ima_norm_um2_data = ima_norm_um2[0].data
    #Individual dimensions for number of x pixels and number of y pixels:
    nxpix_um2_ext1 = ima_norm_um2_hdr['NAXIS1']
    nypix_um2_ext1 = ima_norm_um2_hdr['NAXIS2']
    #Compute the size of the images (you can also do this manually rather than calling these keywords from the header):
    #Call the header and data from the UVOTIMSUM file with the relevant keyword extensions:
    corrfact_um2_ext1 = numpy.zeros((ima_norm_um2_hdr['NAXIS2'], ima_norm_um2_hdr['NAXIS1']))
    coincorr_um2_ext1 = numpy.zeros((ima_norm_um2_hdr['NAXIS2'], ima_norm_um2_hdr['NAXIS1']))

#Check that the dimensions are all the same:
print(corrfact_um2_ext1.shape)
print(coincorr_um2_ext1.shape)
print(ima_norm_um2_data.shape)

# Make a new image file to save the correction factors:
hdu_corrfact = fits.PrimaryHDU(corrfact_um2_ext1, header=ima_norm_um2_hdr)
fits.HDUList([hdu_corrfact]).writeto('.../M33_sum_epoch1_um2_corrfact.img')

# Make a new image file to save the corrected image to:
hdu_coincorr = fits.PrimaryHDU(coincorr_um2_ext1, header=ima_norm_um2_hdr)
fits.HDUList([hdu_coincorr]).writeto('.../M33_sum_epoch1_um2_coincorr.img')

我希望应用以下更正：

    # Define the variables from Poole et al. (2008) "Photometric calibration of the Swift ultraviolet/optical telescope":

alpha =  0.9842000
ft    =  0.0110329
a1    =  0.0658568
a2    = -0.0907142
a3    =  0.0285951
a4    =  0.0308063

for i in range(nxpix_um2_ext1 - 1): #do begin
    for j in range(nypix_um2_ext1 - 1): #do begin
        if (numpy.less_equal(i, 4) | numpy.greater_equal(i, nxpix_um2_ext1-4) | numpy.less_equal(j, 4) | numpy.greater_equal(j, nxpix_um2_ext1-4)): #then begin
            #UVM2
            corrfact_um2_ext1[i,j] == 0
            coincorr_um2_ext1[i,j] == 0
        else:
            xpixmin = i-4
            xpixmax = i+4
            ypixmin = j-4
            ypixmax = j+4
            #UVM2
            ima_UVM2sum = total(ima_norm_um2[xpixmin:xpixmax,ypixmin:ypixmax])
            xvec_UVM2 = ft*ima_UVM2sum
            fxvec_UVM2 = 1 + (a1*xvec_UVM2) + (a2*xvec_UVM2*xvec_UVM2) + (a3*xvec_UVM2*xvec_UVM2*xvec_UVM2) + (a4*xvec_UVM2*xvec_UVM2*xvec_UVM2*xvec_UVM2)
            Ctheory_UVM2 = - alog(1-(alpha*ima_UVM2sum*ft))/(alpha*ft)
            corrfact_um2_ext1[i,j] = Ctheory_UVM2*(fxvec_UVM2/ima_UVM2sum)
            coincorr_um2_ext1[i,j] = corrfact_um2_ext1[i,j]*ima_sk_um2[i,j]

上面的代码片段是搞砸的地方，因为我混合了 IDL 语法和 python 语法。我只是不确定如何将 IDL 的某些方面转换为 python。例如，ima_UVM2sum = total(ima_norm_um2[xpixmin:xpixmax,ypixmin:ypixmax])我不太确定如何处理。

我还遗漏了更新校正因子和巧合校正图像文件的部分，我会说。如果有人有耐心用细齿梳检查一下并提出我需要的必要更改，那就太好了。

原始归一化图像可以在这里下载：Replace ... in above code with this file

Answer 1

关于 numpy 的一个非常重要的事情是它在元素的基础上执行每个数学或比较函数。所以你可能不需要遍历数组。

因此，也许可以从 convolve 图像 sum-filter 开始。这可以通过 astropy.convolution.convolve 或 scipy.ndimage.filters.uniform_filter

对 2D 图像完成

我不确定你想要什么，但我认为你想要一个可以通过

实现的 9x9 求和过滤器

from scipy.ndimage.filters import uniform_filter
ima_UVM2sum = uniform_filter(ima_norm_um2_data, size=9)

因为您想丢弃边界处的任何像素（4 像素），您可以简单地 slice 将它们移除：

ima_UVM2sum_valid = ima_UVM2sum[4:-4,4:-4]

这里忽略了前4行和后4列（最后是通过使停止值为负来实现的）

现在您要计算修正值：

xvec_UVM2 = ft*ima_UVM2sum_valid
fxvec_UVM2 = 1 + (a1*xvec_UVM2) + (a2*xvec_UVM2**2) + (a3*xvec_UVM2**3) + (a4*xvec_UVM2**4)
Ctheory_UVM2 = - np.alog(1-(alpha*ima_UVM2sum_valid*ft))/(alpha*ft)

这些都是数组，所以你仍然不需要循环。

但是你想要填充你的两个图像。要小心，因为校正较小（我们忽略了第一个和最后一个rows/columns）所以你必须在校正图像中采用相同的区域：

corrfact_um2_ext1[4:-4,4:-4] = Ctheory_UVM2*(fxvec_UVM2/ima_UVM2sum_valid)
coincorr_um2_ext1[4:-4,4:-4]  = corrfact_um2_ext1[4:-4,4:-4] *ima_sk_um2

仍然没有循环，只使用 numpys 数学函数。这意味着它要快得多（快得多！）并且做同样的事情。

也许我忘记了一些切片，这会产生 Not broadcastable error 如果是这样请报告。

请注意您的循环：Python 的第一个轴是 FITS 中的第二个轴，第二个轴是第一个 FITS 轴。因此，如果您需要在轴上循环，请记住这一点，这样您就不会得到 IndexErrors 或意外结果。

使用 astropy.fits 和 numpy 对 SWIFT 拟合图像应用巧合校正

Using astropy.fits and numpy to apply coincidence corrections to SWIFT fits image

python

numpy

idl-programming-language

astropy