Support GALPROP spectral cubes?

[cdeil]

In gamma-ray astronomy spectral cubes are used where the spectral axis is log(energy), but a BINTABLE extension is used to define the spectral bins ... I think the spectral WCS info from the cube FITS extension is never used.

E.g. the Fermi LAT and GALPROP are using this format.

Here's an example:

$ wget http://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/gll_iem_v02.fit
$ ipython
...
In [1]: from spectral_cube import SpectralCube
WARNING: ConfigurationDefaultMissingWarning: Could not determine version of package spectral_cube Cannot install default profile. If you are importing from source, this is expected. [spectral_cube._astropy_init]

In [2]: SpectralCube.read("gll_iem_v02.fit")
/Users/deil/Library/Python/2.7/lib/python/site-packages/spectral_cube-0.2.dev389-py2.7.egg/spectral_cube/cube_utils.py:80: UserWarning: No spectral axis found; header may be non-compliant.
  warnings.warn("No spectral axis found; header may be non-compliant.")
Out[2]: 
SpectralCube with shape=(30, 360, 720):
 n_x: 720  type_x: GLON-CAR  unit_x: deg
 n_y: 360  type_y: GLAT-CAR  unit_y: deg
 n_s: 30  type_s: photon energy  unit_s: MeV

In [4]: from astropy.table import Table

In [5]: Table.read('gll_iem_v02.fit')
Out[5]: 
<Table rows=30 names=('Energy')>
array([(49.99999999999999,), (64.98283059038403,), (84.45536543077102,),
       (109.76297408473539,), (142.6541750009015,), (185.40144174189157,),
       (240.9582095985259,), (313.16293027406516,), (407.00427290375654,),
       (528.9657943133443,), (687.4738919994396,), (893.4799891822312,),
       (1161.217175458542,), (1509.1835798293337,), (1961.4204179567785,),
       (2549.173014732112,), (3313.049563238307,), (4305.826770109215,),
       (5596.096231070952,), (7273.003467023406,), (9452.407043617155,),
       (12284.883311734517,), (15966.129821381603,), (20750.486187338385,),
       (26968.506571598286,), (35049.7978763566,), (45552.70155252967,),
       (59202.86975844718,), (76943.40111955488,), (100000.0,)], 
      dtype=[('Energy', '>f8')])

Is it possible to support this format in SpectralCube? Or should I subclass SpectralCube in gammapy?

cc @ellisowen @adonath

[astrofrog]

@cdeil - the I/O framework in spectralcube will transition to using the astropy I/O framework eventually, which will make it easy for you to specify your own readers/writers. I guess at the moment I'm not sure how we can store the table of energies in the WCS object (in theory wcslib supports tabulated WCS, but I'm not sure how it's used).

How do you envisage a reader for this format working? Can you currently write a reader function that returns a SpectralCube, or is the way the data is stored in SpectralCube not suited to this file format?

[cdeil]

I just started looking at spectral-cube an hour ago ... the answer to all your questions above is currently "I have no idea".

Since I need to work with these gamma-ray spectral cubes now, what I'll do as a temporary solution is sub-class or wrap SpectralCube in a GammaSpectralCube class in gammapy and implement what I need there, adding the energy table as a data member.

[cdeil]

I've implemented a first version of the GammaSpectralCube and LogEnergyAxis with some of the functionality we need here: gammapy/gammapy#118

I did not manage to re-use code from SpectralCube so far, probably mainly because I did not manage to put my energy axis into the WCS object but have it as a separate data member, but also because the core functionality I needed right now (interpolating all axes and integrating the energy axis in a specific way) wasn't available in SpectralCube yet.

If you have any comments in the (probably pretty bad) GammaSpectralCube and LogEnergyAxis implementation or hints how I might re-use SpectralCube, that would be very welcome.

Related feature request: astropy/astropy#2362