Rebin Behavior

[joshwalawender]

The rebin function isn't doing what I expect from a binning process on a CCD. Here's an example:

In [31]: arr1 = CCDData([[1,2,3,4], [11,12,13,14], [21,22,23,24], [31,32,33,34]], unit=u.adu)

In [32]: arr1
Out[32]:
CCDData([[ 1,  2,  3,  4],
         [11, 12, 13, 14],
         [21, 22, 23, 24],
         [31, 32, 33, 34]])

In [33]: rebin(arr1, (2,2))
Out[33]:
CCDData([[ 1,  3],
         [21, 23]])

It looks like it is taking only the upper left element of the original pixel grid and using it as the value for the output pixel. I would expect this to do an average or some other combination on the pixel values to get the value for the output pixel. i.e. I would expect the output of the above example to be:

CCDData([[ 6.5,  8.5],
         [26.5, 28.5]])

Am I misunderstanding the role of this function or is this the expected behavior?

[MSeifert04]

I think rebin is probably the wrong name for the function. From the code it seems like it's more a replicate-function because it doesn't apply any function to subsections but only takes the calculated elements. For example it works well if the output-shape is bigger than the initial shape:

import numpy as np
from ccdproc import rebin as ccdproc_rebin

a = np.array([[1,2,3,4], [11,12,13,14], [21,22,23,24], [31,32,33,34]], dtype=np.float32)
ccdproc_rebin(a, (8,8))
array([[  1.,   1.,   2.,   2.,   3.,   3.,   4.,   4.],
       [  1.,   1.,   2.,   2.,   3.,   3.,   4.,   4.],
       [ 11.,  11.,  12.,  12.,  13.,  13.,  14.,  14.],
       [ 11.,  11.,  12.,  12.,  13.,  13.,  14.,  14.],
       [ 21.,  21.,  22.,  22.,  23.,  23.,  24.,  24.],
       [ 21.,  21.,  22.,  22.,  23.,  23.,  24.,  24.],
       [ 31.,  31.,  32.,  32.,  33.,  33.,  34.,  34.],
       [ 31.,  31.,  32.,  32.,  33.,  33.,  34.,  34.]], dtype=float32)

Just for comparison I think the functionality you want is avaiable for "numpy-arrays"-only in astropy/skikit-image and astroscrappy:

from skimage.measure import block_reduce as skimage_rebin  # this is also used by astropy.nddata.utils.block_reduce
from astroscrappy import rebin as astroscrappy_rebin  # assumes 2x2 block-size

astroscrappy_rebin(a)  # can only handle float32 and assumes 2x2 blocks
array([[  6.5,   8.5],
       [ 26.5,  28.5]], dtype=float32)

skimage_rebin(a, (2,2), np.mean)  # (2,2) is block size not final shape! instead of mean also sum could be used.
array([[  6.5,   8.5],
       [ 26.5,  28.5]], dtype=float32)

[crawfordsm]

Yes the rebinning in ccdproc uses slicing. There is a private _blkavg function that would do more of what you want but I think we didn't include it as a public function at the time as we weren't sure exactly how to handle the uncertainty and masking (and hence the rebinning uses slicing . It is probably something that would be good to revisit now, but at the very least, the doc string to rebin should be updated to indicate that it is slicing and not an average.

[joshwalawender]

Thanks for the info. FYI, it looks like the astroscrappy_rebin is faster than the skimage_rebin:

In [90]: %timeit skimage_rebin(arr1, (2,2), np.mean)
10000 loops, best of 3: 123 µs per loop

In [91]: %timeit astroscrappy_rebin(arr1.data.astype(np.float32))
The slowest run took 9.63 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 2.99 µs per loop

Interestingly, skimage_rebin took a CCDData object as input.

I'll play with uncertainty and masking and if I get something that works for me, I'll probably do a PR because I think this would be a common task that people want to do.

[mwcraig]

@joshwalawender -- any progress on this? We are currently thinking that we will deprecate the current rebin, replacing with block_reduce and block_replicate from astropy.nddata.utils.

Those completely ignore masking and uncertainty, though...and we are having trouble figuring a generic way to do that and are thinking instead about documenting how you would apply reduce/replicate to uncertainty and/or mask.

[joshwalawender]

@mwcraig -- No progress. I put in a work around in my code for now, but haven't made any progress on a general tool.