Possible updates to combiner class

[SaOgaz]

Okay so, @eteq recommended I run this stuff by you guys before I started on anything. I've looked through the current stuff in the Combiner class and it looks like there's two things going on.

1) the combinations are being done with numpy masked array arithmetic calls. If these were changed to use the NDData arithmetic it would take care of the uncertainty arrays automatically. Are you guys okay with me changing these calls?

2) in addition to that, the "use the provided uncertainty" arrays when calculating the final uncertainty needs to be hooked into the actual average_combine,sum_combine, median_combine methods one way or another. I'm not sure if this also involves some tweaks to the combine function.

To get what I need out of Combiner I need to do 2 wether or not I do 1. Thoughts?

[SaOgaz]

In reference to #562.

[MSeifert04]

I have a few questions:

• How would you use NDData arithmetic for median_combine?
• Is it useful to use the uncertainty of the inputs? I always thought that in case of combination you're not interested in the propagated uncertainty but the uncertainty of the combined-image (different concept but I can understand why one could want the propagated uncertainty).
• How slow/fast would it be if we replaced the numpy functions with NDData arithmetic operations? combine is quite slow and memory expensive already, if it's a significant slowdown (say ~2-3 times slower or more) then it's probably worth providing it as alternative method/function instead of replacing the original.

[SaOgaz]

Well, median_combine would the obvious exception to using NDData arithmetic.

Using the uncertainties is how a lot of the pipelines for HST data work, and how the IRAF imcombine task is most frequently run for calibration. Uncertainties in this case being the HST Error extensions.

As far as the speed question, I'm not really sure. Any idea about this @eteq?

[MSeifert04]

I thought imcombine also just uses the standard deviation:

Sigma of Combined Pixels The sigma parameter produces optional output sigma image(s). The sigma is the standard deviation, corrected for a finite population, of the input pixel values (excluding rejected pixels) about the output combined pixel values.

[MSeifert04]

I thought about this some more and I think using arithmetic operations with uncertainty propagation isn't correct and the currently implemented way is the "correct one". Because you expect the combination to give a closer representation of the "real value" and then the standard deviation represents how "trustworthy" that combined value is. You don't expect the individual pixels to have a well-defined standard deviation (in most cases it's just an approximation like readout noise, etc.)

However, the more I think about it, we probably need to use ddof=1 when calculating the standard deviation ...

[SaOgaz]

Ahh, you're right, I was conflating imcombine and mscombine (a version of imcombine written specifically for HST data).

[MSeifert04]

@SaraOgaz No problem 😄

I just wondered: Is there some documentation for mscombine regarding the sigma/uncertainty output? I only found http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?mscombine but I couldn't find anything relating to uncertainty propagation there.

[eteq]

Is it useful to use the uncertainty of the inputs? I always thought that in case of combination you're not interested in the propagated uncertainty but the uncertainty of the combined-image (different concept but I can understand why one could want the propagated uncertainty).

I thought about this some more and I think using arithmetic operations with uncertainty propagation isn't correct and the currently implemented way is the "correct one".

@MSeifert04 - I think it very much is useful for at least some combiner use cases. Consider the case of combining 2 or 3 exposures (which is probably ~50% of cases): if you have uncertainties of the individual exposures, using that information is much better information than the stddev of 2 or 3 exposures. Of course sometimes you can't trust the per-image uncertainties, but the combiner machinery should accept the use case where you do.

[eteq]

(Oh and whether or not imcombine or mscombine support this is not critical to our decision here: we want to support IRAF cases, but not necessarily follow all the IRAF conventions...)

[SaOgaz]

@MSeifert04 You can find a little more detail in the gcombine help: http://stsdas.stsci.edu/cgi-bin/gethelp.cgi?gcombine. Mscombine was a re-write/wrapper around gcombine.

[SaOgaz]

@eteq, agreed. @MSeifert04, my main concern here is that the Instrument teams at STScI still use this type of image combination for creating some reference files.

[eteq]

Oh, and to clarify: I see @MSeifert04's point in https://github.com/astropy/ccdproc/issues/569#issuecomment-338396379 that the current implementation is sometimes more correct. So this needs to be an option, as it depends on the situation which is more correct. My inclination is to use the presence/absence of uncertainty as a way to determine which should be used, but it could also be an explicit option to combiners.

[eteq]

(Might also be good to get opinions from @crawfordsm and/or @mwcraig on this as it has some significant implications for some of the science use cases...)

(this = should the combiners support using the uncertainties in individual frames, and how should that be exposed to users)

[MSeifert04]

Ah, thank you for the clarification @eteq. I was thinking too much about bias and flat fields and which of the two approaches is the correct one - because I understood the original issue to propose rewriting the current implementation with the propagate-uncertainty implementation.

Would it make more sense to add a sum and mean/average function unrelated to the Combiner class or are the combine-features like clipping and scaling useful in this context?

[eteq]

The combine-like features are definitely useful in some of these use cases. Also it's less mental burden on users if we don't have two different sets of Combiner-like things. Rather there should be some way to choose the propagation method.

After some conversation with @SaraOgaz I think I lean towards:

• if uncertainty is present and something stddev-like, use the nddata propagation or at least something that uses the uncertainty.
• If its missing, use the approach already here
• if it is present but not stddev-like, raise an error (or maybe a warning, although I favor an error with clear instructions about how to proceed)

[crawfordsm]

I'd agree having the choice of how to propagate the error frame is useful and that it could be right depending on your use case. I'd also agree for the median case that it does not seem to make very much sense so it should only be implemented for the average and sum functionality.

I think my preferred method is to not have the task guess at what it is suppose to do but have a specific keyword passed to it with how it is suppose to be behave. This is far more explicit and would allow people not to have to change their data in case they wanted to use one or the other (ie. like my flats might have error planes but I'd rather have the std measured). This also doesn't break backward compatibility.

The good thing is that average_combine allows you to pass in the function used for averaging, and instead of ~numpy.ma.average, a new function like nddata.average could be passed that handled the combination with propagating the error. For the time being, it doesn't matter where that is implemented.

However, there is some aspects to this that make it non-trivial especially based on the currently construction of Combiner. The uncertainty isn't currently being stored and that would increase the memory limitations. There is also the question of how the code switches between the two. And the 'combine' function would need to be updated. So implementing this would be a non-trivial amount of work.

On the other hand, a case with limited functionality might be easier to implement in combine just using the NDData arithmetic. However, it depends on how important having access to weighting and rejection is to you.

See also my comment on the other thread, but it might time to look at re-factoring all of Combiner to have better performance and handle some of these use cases. It might also also make sense to refactor the xxx_combine tasks to be complete generalized since we are already passing in the function used to do the combination which might be a side step to #578.

[thespacedoctor]

@mwcraig and @crawfordsm - firstly, many thanks for this package and the hours of work you've put into it. Fantastic resource - well done!

Secondly, I'd like to add my +1 to having the option for error propagation of the single-frame error maps into the final-combined map. This was the behaviour I expected and was quite confused about the results I was seeing until I read this issue.