Create uncertainty mixin and SpectralCube with uncertainty class

[brechmos]

Based on the discussion in https://github.com/radio-astro-tools/spectral-cube/issues/572, there is a desire to have a SpectralCubeWithUncertainty and an uncertainty-attribute-containing mixin.

[brechmos]

@keflavich I am looking into adding the SpectralCubewithUncertainty class (name can be changed). In particular, I am looking at the class method read() and it seems like it would be helpful to check to see if there are 2 (or multiple) fits extensions with data and then provide some "smarts" in order to use one of the extensions as data that holds the uncertainty.

Though, given what I have seen with FITS files, it seems FITS files are a little too, ummm, heterogeneous in how they are formatted in order to provide "smarts" built in, but I suppose I could add a parameter the class method read() that is the extension number of the uncertainty data (but then one would also have to specify if it is variance, standard deviation or inverse variance to understand the type of error is represented in the data).

Do you have an opinion on how to build in "smarts" to the read() class method?

[keflavich]

Right... so, if there are multiple extensions, usually they get names (?). We can search for a second extension, check if it has a name like 'VARIANCE' or similar, or failing that, if it has the same shape and BUNIT (could check for both stddev and variance BUNITs). If it passes those, we declare it an uncertainty and return an uncertainty class. Otherwise, we raise an Exception telling the user they have to specify what the additional extensions are (via kwargs).

[keflavich]

Also, check the dtype; if it's an integer extension, we assume it's some kind of mask and raise an exception. We don't have any concept of bitmap masks at the moment - they kind of have to be defined differently for each case.

[brechmos]

How about "for now" I check to see if one of the kwargs for the classmethod read() is uncertainty, then if it is the value must be a tuple of (HDU number, uncertainty type)?

So, then one could do something like:

cube = SpectralCube.read('adv_with_uncertainty.fits', format='fits', 
                         uncertainty=(1, 'VarianceUncertainty'))

This way, for now, there is far less code to make inferences on things that are going to be difficult to infer. If they pass the "uncertainty" parameter, then the obvious checks can be done (extension exists, data on the HDU extension exists, data is the same size as the flux data.

The second entry of the tuple could be one of VarianceUncertainty, StdDevUncertainty, or InverseVariance to keep it consistent with astropy.nddata.

You good with that "for now"?

[keflavich]

Yeah, that's fine. Should the uncertainty type be an actual astropy class, then?

[brechmos]

There are a couple of options: 1) the uncertainty is a numpy array (or something like that) but then the spectral cube class or the user would have to remember what type it is and/or make assumptions that it is a standard deviation or 2) leverage the uncertainties in astropy.nddata and force the user to state the type of uncertainty it is.

By forcing the user to tell us the type of uncertainty then higher level functions (for example) don't need to make assumptions about the type of uncertainty we have, but would know what it is which would calculations easier.

(I guess it comes down to: at some point we are going to have to know the type of uncertainty, and I would lean to force the issue up front.)

[keflavich]

Right, I meant that instead of a user specifying uncertainty=(1, 'uncertainty_type'), they'd specify uncertainty=(1, UncertaintyClass). However, I don't like that callspec much - it means the user has to construct an appropriately structured tuple with no UI-provided clues. Something like uncertainty_extension=1, uncertainty_type=UncertaintyClass would be better.

[brechmos]

Sorry, I didn't mean they would specify an astropy class/type, just a string. I'll mock it up and create a PR for perusal and further discussion.