使用astropy.fits公司和numpy对SWIFT-fits-imag应用一致性校正

Question

这个问题可能有点专业，但希望有人能帮上忙。我通常使用IDL，但对于开发管道，我希望使用python来提高运行时间。你知道吗

我的文件处理设置如下：

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或^{}对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)

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

但是你想填满你的两张照片。请小心，因为校正较小（我们输入了第一行和最后一行/列），所以您必须在校正图像中获取相同的区域：

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或意外的结果。你知道吗