在astropy表列中插值数组

Question

我有一个多波段的辐射源目录（如果你想知道，是从SourceExtractor得到的），我把它读入了一个astropy表格，格式如下：

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个元素，这些元素表示在不同距离源中心的情况下，某个源在某个波段的光子计数（也就是用光圈摄影法测得的）。我有一个光圈的数组（2、3、4、6、8、10、14、20、28、40、60、80、100和160像素），我想把这14个样本插值到一个假设的光圈a上。

我可以逐个处理这些源，但目录里有3000多个源，这样做既不够“pythonic”（即不够符合Python的优雅风格），也不够高效（在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

关于“Pythonic”的概念，主要有几个关键点：简单、易读和实用。如果你的情况真的只是偶尔需要这样做（也就是说，你只会进行3000 x 8次插值几次，而不是一百万次），那么最快且最容易理解的解决方案就是用Python的循环来简单地迭代。这里说的最快，是指从你提出问题到代码给出答案的时间。

在人的时间尺度上，循环和调用一个函数24000次的开销是非常小的，肯定比写一个Stack Overflow的帖子要低得多。:-)

Answer 2

我对astropy表格的格式不太熟悉，但看起来它可以用一个三维的numpy数组来表示，分别对应源、波段和孔径。如果是这样的话，你可以使用比如说scipy.interpolate.interp1d。下面是一个简单的例子。

In [51]: from scipy.interpolate import interp1d

先生成一些示例数据。这个“表格”y是三维的，形状是(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.]])