[BUG] Wrong WCS mapping when reading from Spectrum1D

[WenkeRen]

Jdaviz component

Cubeviz

Description

Following the instruction of Importing Data into Cubeviz, I found there are some inconsistency between the wcs in the Spectrum1D module and that passed to Cubeviz. To me, it seems that Cubeviz reversed the RA and DEC diagonally. Personally, I suspect it could due to the Spectrum1D will somehow rearrange the data sequence from [WAVE, DEC, RA] to [RA, DEC, WAVE], but I do not have ability to locate this issue further.

I gave a simple example for your quick look below.

How to Reproduce

# Import packages I used
import numpy as np
from astropy import wcs
from astropy import units as u
from astropy.coordinates import SpectralCoord, SkyCoord
from specutils import Spectrum1D
import jdaviz
from jdaviz import Cubeviz, Specviz

# Build a random datacube following the instruction of Spectrum1D documents: https://specutils.readthedocs.io/en/stable/spectrum1d.html#slicing
w = wcs.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,
         '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})
spec = Spectrum1D(flux=np.random.default_rng(12345).random((20, 5, 10)) * u.Jy, wcs=w)

# According to this setting, we would expect this data have 20 pixels in wavelength axis, 10 pixels along RA and 5 pixels along DEC.
# We can also do a quick check showing that inside the Spectrum1D, the wcs is not messed:

# According to the wcs settings, the following cut only cut a small slice along RA axis. The data shape would reduced in RA axis
lower = [SpectralCoord(4.89, unit=u.um), SkyCoord(ra=205.439900, dec=26, unit=u.deg)]
upper = [SpectralCoord(4.91, unit=u.um), SkyCoord(ra=205.439697, dec=27.5, unit=u.deg)]
print(f'New datacube shape is: {spec.crop(lower, upper).data.shape}. Only reduced in RA axis')

# For this datacube, we should expect 10 pixel along RA direction and 5 pixel along DEC direction. Put this data into cubeviz we got:
cubeviz = Cubeviz()
cubeviz.load_data(spec)
cubeviz.show()

Expected behavior

Put these code into a jupyter notebook, you will get a interactive plot from cubeviz. From the flux-viewer, we can see the image have 5 pixels along RA direction and 10 along DEC direction. But the input data cube is in reverse.

Browser

Chrome 126.0.6478.127 (arm64)

Jupyter

IPython : 8.26.0
ipykernel : 6.29.5
ipywidgets : 8.1.3
jupyter_client : 8.6.2
jupyter_core : 5.7.2
jupyter_server : 2.14.2
jupyterlab : 4.2.3
nbclient : 0.7.4
nbconvert : 7.16.4
nbformat : 5.10.4
notebook : 7.2.1
qtconsole : 5.5.2
traitlets : 5.14.3

Software versions

macOS-14.4.1-arm64-arm-64bit
Python 3.11.9 | packaged by conda-forge | (main, Apr 19 2024, 18:34:54) [Clang 16.0.6 ]
Numpy 1.26.4
astropy 6.1.1
matplotlib 3.8.4
scipy 1.14.0
scikit-image 0.24.0
asdf 3.3.0
stdatamodels 2.0.0
gwcs 0.21.0
regions 0.9
specutils 1.15.0
specreduce 1.3.0
photutils 1.13.0
astroquery 0.4.7
pyyaml 6.0.1
asteval 1.0.1
idna 3.7
traitlets 5.14.3
bqplot 0.12.43
bqplot-image-gl 1.4.11
glue-core 1.21.1
glue-jupyter 0.22.0
glue-astronomy 0.10.0
echo 0.9.0
ipyvue 1.11.1
ipyvuetify 1.9.4
ipysplitpanes 0.2.0
ipygoldenlayout 0.4.0
ipypopout 1.2.1
Jinja2 3.1.4
voila 0.4.5
vispy 0.14.3
sidecar 0.7.0
Jdaviz 3.10.2

[pllim]

I think this is because Spectrum1D likes to move the spectral axis to a certain place. As a result, the spatial indexing gets messed up in a way that it is (nx, ny) when sliced, while viewer naively expected (ny, nx) in a Pythonic way. Perhaps this won't be such a problem anymore with specutils 2.0, @rosteen ?

Do the test cases here help explain?

jdaviz/jdaviz/conftest.py

Line 102 in 3de1b08

def image_cube_hdu_obj_microns():

jdaviz/jdaviz/conftest.py

Line 206 in 3de1b08

def _create_spectrum1d_cube_with_fluxunit(fluxunit=u.Jy, shape=(2, 2, 4), with_uncerts=False):

[rosteen]

@WenkeRen Thanks for the detailed and reproducible bug report, I'm looking into this. As @pllim said, currently specutils swaps the spectral axis to last, and there's a small possibility we haven't been handling this correctly for the display - most of our test cases are pretty square in the spatial directions and I want to make sure this didn't slip past us.

Either way, yes @pllim specutils 2.0 should make this simpler to handle.

[pllim]

Pretty sure we are transposing internally. Example:

jdaviz/jdaviz/configs/imviz/plugins/coords_info/coords_info.py

Lines 228 to 234 in 3de1b08

	def _image_shape_inds(self, image):
	# return the indices in image.shape for the x and y dimension, respectively
	if image.ndim == 3:
	# cubeviz case
	return (0, 1) # (ix_shape, iy_shape)
	elif image.ndim == 2:
	return (1, 0) # (ix_shape, iy_shape)

[rosteen]

I think this is a real bug, but specifically with the _parse_spectrum1d_3d parser method (the parsers loading FITS files and hdulists and such are fine). I think that the code posted by @WenkeRen shows that this line is incorrectly swapping the RA and Declination in the flux cube:

flux = np.moveaxis(flux, 1, 0)

This is a very easy fix, and in fact I remove that line in the specutils 2.0 compatibility PR (see here), but it seems like it never should have been there in the first place.

I'm going to do some more thinking to make sure it's really not needed.

[pllim]

Oh, thanks for investigating. If we remove it, I think we also have to remove the workarounds we been putting in to compensate for it. 😅

Though, I think it was put in because user was complaining image looked transposed compared to a different cube viewer. We should definitely also double check that too.