3

我有一个辐射源的多波段目录(来自 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 个,这不是非常 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。似乎也应该可以对单个列的内容进行聚合,但我无法理解文档。

有人可以对此有所了解吗?提前致谢!

4

2 回答 2

2

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

In [51]: from scipy.interpolate import interp1d

制作一些样本数据。“表”y是 3D 的,形状为 (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.]])
于 2014-07-12T00:46:42.380 回答
1

关于 Pythonic,其关键方面是简单性、可读性和实用性。如果您的案例确实是一次性的(即您将进行 3000 x 8 插值几次而不是一百万次),那么最快和最容易理解的解决方案将是简单的一个,即使用 Python 循环进行迭代. 我所说的最快是指从你知道你的问题到你从代码中得到答案的那一刻。

在人类/天文学家的时间尺度上,循环和调用函数 24000 次的开销非常小,而且绝对比编写堆栈溢出帖子要低得多。:-)

于 2014-07-12T14:23:05.600 回答