Add block_reduce_hdu function

[cdeil]

@keflavich @astrofrog I'd like to propose to move the block_reduce_hdu function from gammapy.image to this imageutils repo.

I think it's a good example to have what we want in this package (see the description I wrote here and feel free to change it if it's not what you want).

A basic decision for the imageutils API is if we want to take ImageHDU objects as input and output or whether we want to use ndarray and wcs objects separately?

[cdeil]

@astrofrog Unrelated to this issue ... I think travis-ci is scheduling builds for each Github user or org ... for gammapy/gammapy the builds always start immediately, but for this astropy/imageutils they don't, I think because there's a queue from astropy/astropy and the other repos under the astropy org. Not too bad, but could be a reason to not put many more affiliated packages under the astropy org (which or course has other advantages).

[astrofrog]

@cdeil - yeah, I noticed that too... Maybe you could move imageutils under your user for now and move it to the astropy org once it's more mature? (what I'm planning on doing for reproject)

[keflavich]

@cdeil - I think we need to support both APIs - ImageHDU and ndarray + wcs. They should probably live in different sub-namespaces. e.g.:

• imageutils.hdu_utils for wrappers that accept HDUs (block_reduce_hdu)
• imageutils.array_utils for the low-level array manipulation routines (e.g., block_reduce)
• imageutils.wcs_utils for wrappers that accept wcs objects and arrays.

I think we actually need all 3, and maybe even a 4th level to work on filename strings that would be useful for writing command-line scripts. This is the approach I've taken in FITS_tools and I've found it useful so far.

[cdeil]

@astrofrog Let's leave imageutils here for now ... it's a bit annoying, but not really a limiting factor for development. :-)

[cdeil]

@keflavich I'm :+1: on a sub-package for array_utils for functions that only operate on arrays.

I'm not so sure how to best handle the image_hdu vs. (ndarray, wcs) input option. In principle "there should be one way to do it" and having a smaller API, i.e. fewer functions is better.

I guess going back and forth between image_hdu and (ndarray, wcs) is not without it's subtleties?

It looks like @astrofrog prefers to have one function that accepts either (see reproject) and does some checking on input ... any problem with doing that? (it does mean having much fewer functions and I think is OK in terms of complexity for users, especially if done consistently for all our imageutils functions)

[keflavich]

The subtleties of converting between HDU and ndarray+wcs are actually quite awful: lots of metadata is lost in the Header <-> WCS interchange, but most importantly, units are not preserved.

I'm OK with having a top-level wrapper that does the input checking, but I strongly prefer to have a wrapper function that does the checking and then farms out the actions to the appropriate lower-level function. It is important that the low-level functions still be accessible. i.e.:

def magical_wrapper(x):
    if isinstance(x,str):
        return string_function(x)
    elif isinstance(x, Header):
         return hdu_function(x)

[cdeil]

@keflavich Do you have a (preferably simple) example in your codes that is implemented as you prefer (I'm looking for something small we can review and merge first to get us all on the same page API-wise)

[keflavich]

Sort of, yeah. These are the three levels of functions I've implemented for cube regridding; there's no magic top-level wrapper yet though

https://github.com/keflavich/FITS_tools/blob/master/FITS_tools/cube_regrid.py#L26 https://github.com/keflavich/FITS_tools/blob/master/FITS_tools/cube_regrid.py#L63 https://github.com/keflavich/FITS_tools/blob/master/FITS_tools/cube_regrid.py#L100

[cdeil]

@keflavich Hmmm ... I'm not sure having these three regrid_cube functions is useful. It seems that regrid_fits_cube is a few-line wrapper for regrid_cube_hdu which just uses filenames instead of HDUList objects. Then regrid_cube_hdu is a wrapper for regrid_cube which uses HDU objects for input / output instead of having data and header separately? Plus it optionally does smoothing which the other two don't.

I'm not sure having extra functions is worth the API bloat just to save the user from typing a few lines to load their data from file. But the fundamental question is whether most of the image utility functions we write should use (besides the array data) a Header or WCS object or always support both.

Let's see what @astrofrog thinks ... if we are still unsure in a few days we could always ask this one point on the astropy-dev mailing list.

[astrofrog]

It's always a tough question as to how convenient to make it. For example you could have a function that can accept multiple input types and have a convenience function that always sanitizes the input and returns the same thing, so that the overhead of accepting various types of objects is not high. But to some extent that can be added later, so maybe best to focus on writing the core function now and worrying about conveniences later?

[cdeil]

@astrofrog Thanks for your comment.

@keflavich Maybe we should drop the aim to get the API as good as possible for everything we merge here now. Instead we accrete our image utils functions (and do some API discussion / cleanup, but e.g. both ImageHDU and Header and WCS objects are acceptable) and then have a hangout for this API discussion a month or two from now when the core functionality is there already? (could be part of the review when proposing this to be moved into the Astropy core.)

[keflavich]

:+1: As long as we have "low-level" functions that work on arrays and a higher-level version that works on any of [HDU, Header, WCS], I think the approach you suggest is fine. There is some ambiguity in the language I've used here - short answer is, yes, I agree.

[cdeil]

As long as we have "low-level" functions that work on arrays

Here my preference would be to directly contribute those to scikit-image if it's something substantial (of course we'll still have our own utility functions) instead of building our own little scikit-image. Actually this is what I wrote here ... let me know if you disagree with this ... we can certainly still change this. As an example, the scikit-image block-reduce function was added by them after I made a feature request a few months ago and now I only want to have a wrapper here that also block-reduces the Header correctly.

[cdeil]

Should we do a hangout next week for 1 hour max to discuss this? I'd be available Mo to We.

[keflavich]

We can try for that; I'm waiting to hear about the schedule for another telecon first.

Does scikit-image always work nicely with floats?

[cdeil]

We can try for that; I'm waiting to hear about the schedule for another telecon first.

OK ... I'll be offline soon until Sunday evening ... I'll check back with you on Monday if you have time.

Does scikit-image always work nicely with floats?

From what I've seen, yes. They sometimes scale inputs or outputs to the 0 ... 1 range for floats, but usually they have an option to not do this. photutils already uses scikit-image, so I don't think there will be strong objections to having it as a dependency for astropy.image.

[keflavich]

For dealing with a WCS object instead of a header, you need to account for the possibility of wcs having a pc array rather than a cd. @astrofrog - did we have a better idiom for this?

if hasattr(w.wcs, 'cd'):
    w.wcs.cdelt = w.wcs.cdelt * block_size
else:
    assert np.all(w.wcs.cdelt == 1.0) # I think this has to be true for the PC case
    w.wcs.pc = w.wcs.pc * block_size
w.wcs.crpix =  ((w.wcs.crpix-0.5)/block_size)+0.5

I think that approach is right... but it's worth verifying!