使用python将fits文件中的低信号像素归零

Question

我正在尝试编写一个简单的脚本，该脚本将读取适合的图像文件，将低于某个阈值（我将根据背景统计数据输入）的所有像素设置为零或 nan。我在下面包含了我正在使用的代码。这将写入一个类似于原始文件的新拟合文件，但像素值似乎被重新缩放（即，我希望新文件中的最小值是我指定的最小值，但事实并非如此）。我认为这是我如何将标题写入新文件的问题，但不知道如何修复它，并且欢迎任何做过类似事情的天文学家的建议或 sn-ps。谢谢！

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

f=fits.open('filename.fits')
data = f[0].data
header = f[0].header
data[data<noise_cutoff]=np.nan
fits.writeto('outfilename.fits', data, header)
f.close()

Answer 1

在 python 2.7 中，我一直是这样做的：

img[numpy.where(img<min_flux)]=0
img[numpy.where(img>0)]+=add_flux

看看区别：

这会返回满足条件的字段的实际值：

>>> data2[data2>noise_cutoff]
array([ 0.15600586,  0.15576172,  0.16992188, ...,  0.15063477,
        0.21899414,  0.15722656], dtype=float32)

这将返回条件返回True的索引：

>>> np.where(data2>noise_cutoff)
(array([   0,    0,    0, ..., 1488, 1488, 1488]), array([  88,  789, 1065, ..., 1683, 1833, 1872]))

您想将索引处的值设置为 nan 或零，我不知道第一个是如何解释的，但根据经验我知道它不起作用。

编辑 我对处理 nan 值的 FITS 观众也有不好的体验。通常人们会写一个适合字段，即RANGE，它应该描述你所有可能的值。我经常看到[0, 65535]，但这主要是一种协议和数据类型（毕竟它是 FITS 格式）。

我也对发送索引和值之间的区别感到好奇，这就是我所做的：

先显示一些默认值：

>>> data
array([[ 0.01800537,  0.00421143, -0.01644897, ..., -0.03686523,
         0.05981445, -0.00924683],
       [-0.00267029, -0.02334595,  0.03179932, ...,  0.09436035,
         0.05981445,  0.00457001],
       [-0.13354492, -0.0302124 , -0.00266266, ...,  0.05291748,
        -0.06445312,  0.09436035],
       ..., 
       [ 0.04669189, -0.02218628, -0.06347656, ..., -0.01507568,
         0.10229492,  0.02636719],
       [ 0.00536346, -0.00842285,  0.04669189, ..., -0.00816345,
         0.00565338, -0.02886963],
       [-0.07043457, -0.00840759, -0.09106445, ...,  0.06787109,
        -0.11865234, -0.05645752]], dtype=float32)

为了不编辑我需要的原始文件，我制作了一个副本，然后我使用data2[data2>noise_cutoff] 和np.where(data2>noise_cutoff) 操作并将结果复制到[ ] 运算符中，并进行了一些编辑，因此它是一个有效的表达式.

>>> data2 = data.copy

# result of data2[data2>noise_cutoff]
>>> data2[[0.21899414,  0.15722656]]
array([[ 0.01800537,  0.00421143, -0.01644897, ..., -0.03686523,
         0.05981445, -0.00924683],
       [ 0.01800537,  0.00421143, -0.01644897, ..., -0.03686523,
         0.05981445, -0.00924683]], dtype=float32)

#result of np.where(data2>noise_cutoff)
>>> data2[[0,0, 0, 1488, 1488, 1488], [88,789,1065,1683,1833,1872]]
array([ 0.15600586,  0.15576172,  0.16992188,  0.15063477,  0.21899414,
        0.15722656], dtype=float32)

所以 numpy 似乎喜欢解压发送到[ ] 的参数。例如，发送data2[[1, 2]] 将返回data2 的第一行和第二行。由于我已经发送了浮点值0.21899414 和0.15722656，它们显然被转换为ints，向下舍入到0，并且第一行返回了两次。

发送data2[[1],[1]] 但会返回data[1,1] 处的浮点数。发送这两个的列表会返回一系列写入这些索引的值，即：data2[[1,2], [1,2]] 获取元素 [1,1]、[2,2]

>>> data2[1,1] == data2[[1],[1]]
array([ True], dtype=bool)