Representing a HEALPix map with data

[astrofrog]

At the moment, the HEALPix class doesn't include data internally, and instead data should be passed when interpolating. However, it would be nice to be able to represent a HEALPix map including potentially multiple fields, so that we can easily read/write from/to e.g. HEALPix FITS files.

Here is an idea of how to do this (but looking for feedback): we change __init__ on HEALPix so as to optionally also accept a data argument. This can be either a single data array or a dict-like object of arrays (which can include Table). Then when interpolating, if there is a single array, then that array gets interpolated, otherwise a dictionary of interpolated values will be returned (and there could be an option to specify the field to interpolate to save computing time).

We can then use the Astropy I/O registry to register FITS readers/writers (which could also take HDU objects).

@woodmd @cdeil - any thoughts on this?

[woodmd]

Sounds reasonable to me although I would go for a tuple rather than a dict argument.

With respect to FITS I/O it would be nice if this would support the HEALPix format conventions that we created for gamma-ray analysis:
http://gamma-astro-data-formats.readthedocs.io/en/latest/skymaps/healpix/index.html

Here we also have the concept of a partial-sky HEALPix map where we define the geometry of the map with a HPX_REG string that we write to the FITS header. Right now we've only defined disk and single pixel geometries but one could also define more complicated geometries (square or arbitrary list of healpix pixels). Do you think it would be feasible to add support for partial-sky maps in this package?

[astrofrog]

Sounds reasonable to me although I would go for a tuple rather than a dict argument.

Do you mean a named tuple or really just a tuple? (the issue with tuple would be that we'd lose the field names from the FITS table).

Regarding partial maps, I think it would be fine to include here (depending on the amount of code overhead)

[woodmd]

Ah I see you intended that each element would map to a column. I would still prefer something that allows tuple-like access (i.e. by index) so I guess named tuple might be a good compromise.

For the partial maps the overhead should be fairly small. What we should try to decide now is what convention to use for specifying the geometry of a partial-sky map. In the gamma-astro-data-formats spec we use a string so one option would be to add an optional string argument to all methods that can accept a partial sky map.

[astrofrog]

Thinking about this more, a named tuple won't work because some column names aren't necessarily valid Python names.

Could you provide an example file that you'd like to use and show how you would expect the index access to behave? Just to clarify, in some cases HEALPIX maps have a single column, e.g. TEMPERATURE, which is the a vector column - I'm suggesting this would be one field in a dictionary and that you would then indeed use item access to access rows/columns of that field. In my mind it doesn't make sense to access e.g. TEMPERATURE using an index since it would be far more reliable to simply use TEMPERATURE in case the number/order of columns changes in future.

Just to clarify, if I have a file with one column which is then a vector column,

[woodmd]

In gamma-ray applications we often have columns which correspond to image planes of a higher dimensional data structure -- for instance a cube where the third axis is energy. In that case we usually name the columns according to the image plane number (e.g. CHANNEL1, CHANNEL2, etc.). In this scenario being able to access the columns by index rather than name would be convenient although not absolutely critical. Here and here are files that illustrate what this looks like in the case of all- and partial-sky maps.

[astrofrog]

Ah that makes sense - I think we could easily have some kind of object that can be accessed by name or by index.