Interpolating an array within an astropy table column

Question

I have a multiband catalog of radiation sources (from SourceExtractor, if you care to know), which I have read into an astropy table in the following form:

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

The arrays in FLUX_APER_BAND1, FLUXERR_APER_BAND1, etc. each have 14 elements, which give the number of photon counts for a given source in a given band, within 14 different distances from the center of the source (aperture photometry). I have the array of apertures (2, 3, 4, 6, 8, 10, 14, 20, 28, 40, 60, 80, 100, and 160 pixels), and I want to interpolate the 14 samples into a single (assumed) count at some other aperture a.

I could iterate over the sources, but the catalog has over 3000 of them, and that's not very pythonic or very efficient (interpolating 3000 objects in 8 bands would take a while). Is there a way of interpolating all the arrays in a single column simultaneously, to the same aperture? I tried simply applying np.interp, but that threw ValueError: object too deep for desired array, as well as np.vectorize(np.interp), but that threw ValueError: object of too small depth for desired array. It seems like aggregation should also be possible over the contents of a single column, but I can't make sense of the documentation.

Can someone shed some light on this? Thanks in advance!

Answer 1

I'm not familiar with the format of an astropy table, but it looks like it could be represented as a three-dimensional numpy array, with axes for source, band and aperture. If that is the case, you can use, for example, scipy.interpolate.interp1d. Here's a simple example.

In [51]: from scipy.interpolate import interp1d

Make some sample data. The "table" y is 3-D, with shape (2, 3, 14). Think of it as the array holding the counts for 2 sources, 3 bands and 14 apertures.

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]]])

Create the interpolator. This creates a linear interpolator by default. (Check out the docstring for different interpolators. Also, before calling interp1d, you might want to transform your data in such a way that linear interpolation is appropriate.) I use axis=2 to create an interpolator of the aperture axis. f will be a function that takes an aperture value and returns an array with shape (2,3).

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

Take a look at a couple y slices. These correspond to apertures 2 and 3 (i.e. x[0] and 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]])

Use the interpolator to get the values at apertures 2, 2.5 and 3. As expected, the values at 2 and 3 match the values in 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.]])