python - 在 astropy 表列中插入一个数组

Question

我有一个辐射源的多波段目录(来自 SourceExtractor，如果你想知道的话)，我已将其读入以下形式的天体表:

Source # | FLUX_APER_BAND1 | FLUXERR_APER_BAND1  ...  FLUX_APER_BANDN | FLUXERR_APER_BANDN
1           np.array(...)      np.array(...)     ...   np.array(...)      np.array(...)
...

FLUX_APER_BAND1、FLUXERR_APER_BAND1 等中的数组各有 14 个元素，给出给定波段内给定源的光子数，在 14 个不同的范围内距光源中心的距离(孔径光度法)。我有孔径阵列(2、3、4、6、8、10、14、20、28、40、60、80、100 和 160 像素)，我想将 14 个样本插值到一个(假设)在其他一些孔径 a 处计数。

我可以遍历源，但目录有超过 3000 个，这不是很 pythonic 或非常有效(在 8 个波段中插入 3000 个对象需要一段时间)。有没有一种方法可以同时将单个列中的所有阵列插值到相同的孔径？我尝试简单地应用 np.interp，但这引发了 ValueError: object too deep for desired array，以及 np.vectorize(np.interp)，但这引发了 ValueError: object of too small depth for desired array。似乎也应该可以对单个列的内容进行聚合，但我无法理解文档。

有人可以阐明这一点吗？提前致谢!

Answer 1

我不熟悉 astropy 表的格式，但它看起来可以表示为三维 numpy 数组，其中轴为源轴、波段轴和孔径轴。如果是这种情况，您可以使用，例如，scipy.interpolate.interp1d .这是一个简单的例子。

In [51]: from scipy.interpolate import interp1d

制作一些示例数据。 “表格”y 是 3 维的，形状为 (2, 3, 14)。可以将其视为保存 2 个源、3 个波段和 14 个孔径的计数的数组。

In [52]: x = np.array([2, 3, 4, 6, 8, 10, 14, 20, 28, 40, 60, 80, 100, 160])

In [53]: y = np.array([[x, 2*x, 3*x], [x**2, (x+1)**3/400, (x**1.5).astype(int)]])

In [54]: y
Out[54]: 
array([[[    2,     3,     4,     6,     8,    10,    14,    20,    28,
            40,    60,    80,   100,   160],
        [    4,     6,     8,    12,    16,    20,    28,    40,    56,
            80,   120,   160,   200,   320],
        [    6,     9,    12,    18,    24,    30,    42,    60,    84,
           120,   180,   240,   300,   480]],

       [[    4,     9,    16,    36,    64,   100,   196,   400,   784,
          1600,  3600,  6400, 10000, 25600],
        [    0,     0,     0,     0,     1,     3,     8,    23,    60,
           172,   567,  1328,  2575, 10433],
        [    2,     5,     8,    14,    22,    31,    52,    89,   148,
           252,   464,   715,  1000,  2023]]])

创建插值器。默认情况下，这会创建一个线性插值器。 (查看不同插值器的文档字符串。此外，在调用 interp1d 之前，您可能希望以适合线性插值的方式转换数据。)我使用 axis=2 创建孔径轴的插值器。 f 将是一个函数，它接受一个光圈值并返回一个形状为 (2,3) 的数组。

In [55]: f = interp1d(x, y, axis=2)

看看几个 y 切片。这些对应于孔径 2 和 3(即 x[0] 和 x[1])。

In [56]: y[:,:,0]
Out[56]: 
array([[2, 4, 6],
       [4, 0, 2]])

In [57]: y[:,:,1]
Out[57]: 
array([[3, 6, 9],
       [9, 0, 5]])

使用插值器获取孔径 2、2.5 和 3 处的值。正如预期的那样，2 和 3 处的值与 y 中的值匹配。

In [58]: f(2)
Out[58]: 
array([[ 2.,  4.,  6.],
       [ 4.,  0.,  2.]])

In [59]: f(2.5)
Out[59]: 
array([[ 2.5,  5. ,  7.5],
       [ 6.5,  0. ,  3.5]])

In [60]: f(3)
Out[60]: 
array([[ 3.,  6.,  9.],
       [ 9.,  0.,  5.]])