Writer needed for cube data

[robelgeda]

How do we save spectral cubes to file formats that are readable by other applications?
What type of api call do we implement for this?

🐱

[nmearl]

https://specutils.readthedocs.io/en/stable/custom_loading.html#creating-a-custom-writer

[robelgeda]

So the user has to create their own writer or are we planning on providing them with some general writers?

[nmearl]

We have a general cube loader, but not a general cube writer. We'll definitely end up providing one, but if a user requires some specific thing to be stored with their output, they'd probably have to create their own.

[dhomeier]

#632 has a writer that supports 2D spectra now, maybe that could help as a template for a cube version – actually I checked that you can throw a 3D spectrum at it as well, that is read in as a cube by SAOImageDS9, but it assumes identical spectra axes for all subspectra along the 2nd axis.

[dhomeier]

Addendum: I've tried to read in the test file created above, which has

SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                  -64 / array data type
NAXIS   =                    3 / number of array dimensions
NAXIS1  =                   10
NAXIS2  =                    3
NAXIS3  =                    2
WCSAXES =                    1 / Number of coordinate axes
BUNIT   = 'W / (m2 nm)'        / [W / (m2 nm)] spectral flux density wav
CRPIX1  =                516.0 / Pixel coordinate of reference point
CDELT1  =         1.5429985523 / [Angstroms] Coordinate increment at reference p
CUNIT1  = 'Angstroms'          / Units of coordinate increment and value
CTYPE1  = 'LINEAR'             / Coordinate type code
CRVAL1  =         4831.4594727 / [Angstroms] Coordinate value at reference point
LATPOLE =                 90.0 / [deg] Native latitude of celestial pole
MJDREFI =                  0.0 / [d] MJD of fiducial time, integer part
MJDREFF =                  0.0 / [d] MJD of fiducial time, fractional part
END

with the generic_cube reader (after uncommenting the @data_loader("Cube", ... line), which resulted in the following exception:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/sw/lib/python3.8/site-packages/astropy/nddata/mixins/ndio.py", line 56, in __call__
    return registry.read(self._cls, *args, **kwargs)
  File "/sw/lib/python3.8/site-packages/astropy/io/registry.py", line 523, in read
    data = reader(*args, **kwargs)
  File "/Users/derek/lib/python3.8/site-packages/specutils/io/default_loaders/generic_cube.py", line 49, in generic_fits
    data = data3[:,iy,ix]
IndexError: index 516 is out of bounds for axis 2 with size 10

so this might not be the desired type of cube, or the comment about generic_cube being "not yet ready" still holds...

[pllim]

I was just looking for this. This is blocking Jdaviz development because when we generate a new spectral cube, we can no longer write it out properly after switching from spectral-cube package to specutils. @keflavich , any advise to push this forward since you are familiar with both packages? 🙏

Problem:

import numpy as np
from astropy import units as u
from astropy.io import fits
from astropy.wcs import WCS
from specutils import Spectrum1D

flux = np.arange(720).reshape((8, 9, 10)).astype(np.float32) * u.nJy
w = WCS({'WCSAXES': 3, 'CRPIX1': 38.0, 'CRPIX2': 38.0, 'CRPIX3': 1.0,
         'CRVAL1': 205.4384, 'CRVAL2': 27.004754, 'CRVAL3': 4.890499866509344,
         'CTYPE1': 'RA---TAN', 'CTYPE2': 'DEC--TAN', 'CTYPE3': 'WAVE',
         'CUNIT1': 'deg', 'CUNIT2': 'deg', 'CUNIT3': 'um',
         'CDELT1': 3.61111097865634E-05, 'CDELT2': 3.61111097865634E-05, 'CDELT3': 0.001000000047497451,  # noqa
         'PC1_1 ': -1.0, 'PC1_2 ': 0.0, 'PC1_3 ': 0,
         'PC2_1 ': 0.0, 'PC2_2 ': 1.0, 'PC2_3 ': 0,
         'PC3_1 ': 0, 'PC3_2 ': 0, 'PC3_3 ': 1,
         'DISPAXIS': 2, 'VELOSYS': -2538.02,
         'SPECSYS': 'BARYCENT', 'RADESYS': 'ICRS', 'EQUINOX': 2000.0,
         'LONPOLE': 180.0, 'LATPOLE': 27.004754,
         'MJDREFI': 0.0, 'MJDREFF': 0.0, 'DATE-OBS': '2014-03-30'})
sp = Spectrum1D(flux=flux, wcs=w)
sp.write("blah3.fits")  # This write out the wrong format

Filename: blah3.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU       4   ()      
  1                1 BinTableHDU     15   8R x 2C   [D, 90E]

What I really want is to write it out like this but using specutils.Spectrum1D.write():

hdu = fits.ImageHDU(flux.value)
hdu.name = "SCI"
hdu.header['BUNIT'] = flux.unit.to_string(format='fits')
hdu.header.update(w.to_header())
hdulist = fits.HDUList([fits.PrimaryHDU(), hdu])
hdulist.writeto("blah4.fits")

Filename: blah4.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU       4   ()      
  1  SCI           1 ImageHDU        35   (10, 9, 8)   float32

Motivation: We want this functionality to be able to roundtrip cube I/O in Jdaviz (Cubeviz).

Workaround: We can stop using specutils.Spectrum1D.write() and make our own writer but given that this is such a common thing that used to be supported by spectral-cube, it would discourage people to use specutils if this isn't provided out of the box.

[dhomeier]

The WCS in that example is 1D? What happens if you write it out with
sp.write("blah3.fits", hdu=0)?

[pllim]

No, it is a 3D WCS. It is a legit spectral cube.

[pllim]

Adding hdu=0 gives an error: ValueError: Wrong number of dimensions in input array. Expected 2.

[pllim]

That example above is a fully reproducible example. I gotta run now but feel free to play with it.

[dhomeier]

No, it is a 3D WCS. It is a legit spectral cube.

I mean 1D-spectral. Didn't think it could become a legit Spectrum1D otherwise.

[keflavich]

I'm listening in here, but unfortunately I don't have the answer; we're not close to having a solution that merges spectral-cube and specutils again. If this is now a priority for STSCI, I'd be interested in talking through how to do that, though.

Spectral-cube's writer is extremely simple, but it is also built on spectral-cubes being closely linked to FITS-WCS.

Writer:
https://github.com/radio-astro-tools/spectral-cube/blob/master/spectral_cube/io/fits.py#L265-L284

HDU builder:
https://github.com/radio-astro-tools/spectral-cube/blob/master/spectral_cube/spectral_cube.py#L2525-L2536

so I imagine most of the effort would be in mangling specutils cubes in to spectral-cube form.

[pllim]

I mean 1D-spectral

Yes, 1D has spectral, while the other 2D is spatial (e.g., RA and Dec).

[dhomeier]

I can get the writer to pass the spectral_axis check by comparing to
wcs.spectral.all_pix2world(np.arange(len(wl)), 0), but that would also only work for FITS-WCS and needs more work not to break other WCSs. And the file it's writing apparently is still no valid 3D WCS:

Header listing for 1. HDU of type IMAGE_HDU:
SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                  -32 / array data type
NAXIS   =                    3 / number of array dimensions
NAXIS1  =                    8
NAXIS2  =                    9
NAXIS3  =                   10
WCSAXES =                    3 / Number of coordinate axes
BUNIT   = 'nJy     '           / [nanoJansky] spectral flux density
CRPIX1  =                  1.0 / Pixel coordinate of reference point
CRPIX2  =                 38.0 / Pixel coordinate of reference point
CRPIX3  =                 38.0 / Pixel coordinate of reference point
PC3_3   =                 -1.0 / Coordinate transformation matrix element
CDELT1  =  1.0000000474975E-09 / [m] Coordinate increment at reference point
CDELT2  =  3.6111109786563E-05 / [deg] Coordinate increment at reference point
CDELT3  =  3.6111109786563E-05 / [deg] Coordinate increment at reference point
CUNIT1  = 'm'                  / Units of coordinate increment and value
CUNIT2  = 'deg'                / Units of coordinate increment and value
CUNIT3  = 'deg'                / Units of coordinate increment and value
CTYPE1  = 'WAVE'               / Vacuum wavelength (linear)
CTYPE2  = 'DEC--TAN'           / Declination, gnomonic projection
CTYPE3  = 'RA---TAN'           / Right ascension, gnomonic projection
CRVAL1  =  4.8904998665093E-06 / [m] Coordinate value at reference point
CRVAL2  =            27.004754 / [deg] Coordinate value at reference point
CRVAL3  =             205.4384 / [deg] Coordinate value at reference point
LONPOLE =                180.0 / [deg] Native longitude of celestial pole
LATPOLE =            27.004754 / [deg] Native latitude of celestial pole
BUNIT   = 'nJy     '           / [nanoJansky] spectral flux density
CRPIX1  =                  1.0 / Pixel coordinate of reference point
CRPIX2  =                 38.0 / Pixel coordinate of reference point
CRPIX3  =                 38.0 / Pixel coordinate of reference point
PC3_3   =                 -1.0 / Coordinate transformation matrix element
CDELT1  =  1.0000000474975E-09 / [m] Coordinate increment at reference point
CDELT2  =  3.6111109786563E-05 / [deg] Coordinate increment at reference point
CDELT3  =  3.6111109786563E-05 / [deg] Coordinate increment at reference point
CUNIT1  = 'm'                  / Units of coordinate increment and value
CUNIT2  = 'deg'                / Units of coordinate increment and value
CUNIT3  = 'deg'                / Units of coordinate increment and value
CTYPE1  = 'WAVE'               / Vacuum wavelength (linear)
CTYPE2  = 'DEC--TAN'           / Declination, gnomonic projection
CTYPE3  = 'RA---TAN'           / Right ascension, gnomonic projection
CRVAL1  =  4.8904998665093E-06 / [m] Coordinate value at reference point
CRVAL2  =            27.004754 / [deg] Coordinate value at reference point
CRVAL3  =             205.4384 / [deg] Coordinate value at reference point
LONPOLE =                180.0 / [deg] Native longitude of celestial pole
LATPOLE =            27.004754 / [deg] Native latitude of celestial pole
MJDREF  =                  0.0 / [d] MJD of fiducial time
DATE-OBS= '2014-03-30'         / ISO-8601 time of observation
MJD-OBS =              56746.0 / [d] MJD of observation
RADESYS = 'ICRS'               / Equatorial coordinate system
SPECSYS = 'BARYCENT'           / Reference frame of spectral coordinates
VELOSYS =             -2538.02 / [m/s] Velocity towards source
END

fails on read with

E astropy.wcs._wcs.InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3105 of file cextern/wcslib/C/wcs.c:
E Unmatched celestial axes.

[pllim]

spectral-cubes being closely linked to FITS-WCS

Maybe this is good enough for now until we have demands for cubes with GWCS. Writer that works for most cases but not all is better than no writer at all?

[dhomeier]

Writer that writes out illegal WCS is not much use either. Although directly reading the header returns a valid WCS, so it's probably something in the loader that needs to be fixed.

[keflavich]

Maybe this is good enough for now until we have demands for cubes with GWCS. Writer that works for most cases but not all is better than no writer at all?

That's certainly my take!

I have also encountered no demand for support for non-FITS-compatible cubes; generally the process of building a cube, at any wavelength, requires putting the data onto a regular grid, and thus FITS can handle it. I'm sure one can come up with use cases that require GWCS, but I haven't seen them in the wild. Solar people probably have counterexamples though.

That said, I'm not clear what's wrong with the output @dhomeier showed. Maybe it's that rogue PC3_3? Otherwise this WCS looks healthy to me?

[pllim]

Re: Illegal WCS -- I am not sure either. I thought the example code I gave to write it back out was fine, as Cubeviz was able to roundtrip with that one.

[dhomeier]

Yes, it was the the loader that "fixed" the WCS when, I think, only 1D (pure spectral) WCS could be passed to Spectrum1D.
That restriction can probably be lifted now. WIP for a fix in #1009.

[pllim]

I have a downstream workaround at spacetelescope/jdaviz#2012 until this issue is resolved. Please have a look to see if I missed anything. Thanks!