BUG: Cubeviz parser gives wrong sky coordinates for non-flux cube under certain condition #1991

pllim · 2023-02-06T22:40:31Z

Should probably wait till after this PR is merged before attempting to fix:

refactor mouseover coords info display #1976

Previous failed attempt to address this problem:

Fix Cubeviz coordinates panel with inaccurate sky coordinates sometimes #2002

Related to this headache:

How to use PaddedSpectrumWCS on 3D cube? glue-viz/glue-astronomy#75

Might cause more headache later:

Allow spectral axis to be anywhere, instead of forcing it to be last astropy/specutils#999

Consider this test case:

jdaviz/jdaviz/configs/cubeviz/plugins/tests/test_parsers.py

Line 26 in 3fc9a92

    
           def test_fits_image_hdu_with_microns(image_cube_hdu_obj_microns, cubeviz_helper):

jdaviz/jdaviz/conftest.py

Line 95 in 3fc9a92

def image_cube_hdu_obj_microns():

import numpy as np
from astropy.io import fits

from jdaviz import Cubeviz


def image_cube_hdu_obj_microns():
    flux_hdu = fits.ImageHDU(np.ones((8, 9, 10)).astype(np.float32))
    flux_hdu.name = 'FLUX'

    uncert_hdu = fits.ImageHDU(np.zeros((8, 9, 10)).astype(np.float32))
    uncert_hdu.name = 'ERR'

    mask_hdu = fits.ImageHDU(np.ones((8, 9, 10)).astype(np.uint16))
    mask_hdu.name = 'MASK'

    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'}

    flux_hdu.header.update(wcs)
    flux_hdu.header['BUNIT'] = '1E-17 erg*s^-1*cm^-2*Angstrom^-1'

    uncert_hdu.header['BUNIT'] = '1E-17 erg*s^-1*cm^-2*Angstrom^-1'

    return fits.HDUList([fits.PrimaryHDU(), flux_hdu, uncert_hdu, mask_hdu])


cubeviz_helper = Cubeviz()
cubeviz_helper.load_data(image_cube_hdu_obj_microns())

Once loaded, you will see that "flux" cube has PaddedSpectrumWCS from glue-astronomy but not the "uncert" cube.

>>> cubeviz_helper.app.data_collection[0].coords.__class__.__name__
'PaddedSpectrumWCS'
>>> cubeviz_helper.app.data_collection[1].coords.__class__.__name__
'WCS'

I forgot exactly what are the series of decisions but I think this was originally added in #1480 .

jdaviz/jdaviz/configs/cubeviz/plugins/parsers.py

Line 153 in 3fc9a92

metadata['_orig_wcs'] = wcs

In Cubeviz coordinates info display, _orig_wcs is favored (or moment map and also maybe because the PaddedSpectrumWCS does not have a functioning pixel_to_world method (as required by APE 14).

>>> cubeviz_helper.app.data_collection[0].coords.pixel_to_world(0, 0, 0)
[<SpectralCoord 4.89049987 um>, <Quantity 0. pix>, <Quantity 1. pix>]
>>> cubeviz_helper.app.data_collection[1].coords.pixel_to_world(0, 0, 0)
[<SpectralCoord 
    (target: <ICRS Coordinate: (ra, dec, distance) in (deg, deg, kpc)
                 (205.4384, 27.004754, 1000.)
              (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
                 (0., 0., 0.)>)
   4.89049987e-06 m>,
 <SkyCoord (ICRS): (ra, dec) in deg
     (205.4398996, 27.00341788)>]

Unfortunately, somewhere in the parsing (possibly by specutils), the WCS axes are swapped. So now we have this:

>>> cubeviz_helper.app.data_collection[0].meta['_orig_wcs']
WCS Keywords

Number of WCS axes: 3
CTYPE : 'RA---TAN'  'DEC--TAN'  'WAVE'  
CRVAL : 205.4384  27.004754  4.890499866509344e-06  
CRPIX : 38.0  38.0  1.0  
PC1_1 PC1_2 PC1_3  : -1.0  0.0  0.0  
PC2_1 PC2_2 PC2_3  : 0.0  1.0  0.0  
PC3_1 PC3_2 PC3_3  : 0.0  0.0  1.0  
CDELT : 3.61111097865634e-05  3.61111097865634e-05  1.000000047497451e-09  
NAXIS : 10  9  8

>>> cubeviz_helper.app.data_collection[1].coords
WCS Keywords

Number of WCS axes: 3
CTYPE : 'WAVE'  'DEC--TAN'  'RA---TAN'  
CRVAL : 4.890499866509344e-06  27.004754  205.4384  
CRPIX : 1.0  38.0  38.0  
PC1_1 PC1_2 PC1_3  : 1.0  0.0  0.0  
PC2_1 PC2_2 PC2_3  : 0.0  1.0  0.0  
PC3_1 PC3_2 PC3_3  : 0.0  0.0  -1.0  
CDELT : 1.000000047497451e-09  3.61111097865634e-05  3.61111097865634e-05  
NAXIS : 8  9  10

As a result, assuming (x, y, z) input is now wrong for "uncert" (and by similar logic, "mask") cube:

>>> # x=1 y=0 z=0
>>> cubeviz_helper.app.data_collection[0].meta['_orig_wcs'].pixel_to_world(1, 0, 0)
[<SkyCoord (ICRS): (ra, dec) in deg
     (205.43985907, 27.00341788)>,
 <SpectralCoord 
    (target: <ICRS Coordinate: (ra, dec, distance) in (deg, deg, kpc)
                 (205.4384, 27.004754, 1000.)
              (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
                 (0., 0., 0.)>)
   4.89049987e-06 m>]
>>> # WRONG: z=1 y=0 x=0
>>> cubeviz_helper.app.data_collection[1].coords.pixel_to_world(1, 0, 0)
[<SpectralCoord 
    (target: <ICRS Coordinate: (ra, dec, distance) in (deg, deg, kpc)
                 (205.4384, 27.004754, 1000.)
              (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
                 (0., 0., 0.)>)
   4.89149987e-06 m>,
 <SkyCoord (ICRS): (ra, dec) in deg
     (205.4398996, 27.00341788)>]
>>> # SHOULD HAVE BEEN: z=0 y=0 x=1
>>> cubeviz_helper.app.data_collection[1].coords.pixel_to_world(0, 0, 1)
[<SpectralCoord 
    (target: <ICRS Coordinate: (ra, dec, distance) in (deg, deg, kpc)
                 (205.4384, 27.004754, 1000.)
              (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
                 (0., 0., 0.)>)
   4.89049987e-06 m>,
 <SkyCoord (ICRS): (ra, dec) in deg
     (205.43985907, 27.00341788)>]

Proposed fix: Use coords.celestial component to do pixel_to_world translation though one needs to figure out if it is always (x, y) by definition. Also does GWCS have this component?

If proposed fix does not work, will astropy/specutils#999 naturally solve this upstream?

🐱

The text was updated successfully, but these errors were encountered:

pllim · 2023-02-07T23:00:56Z

For GWCS, Mihai suggested I look at https://gwcs.readthedocs.io/en/latest/api/gwcs.wcs.WCS.html#gwcs.wcs.WCS.available_frames though I am not optimistic those lower level attributes follow APE 14.

pllim · 2023-02-09T16:44:47Z

.celestial isn't as magical as I hoped. It has the same problem with swapped axes.

>>> wcs0 = cubeviz_helper.app.data_collection[0].meta['_orig_wcs']
>>> wcs0.celestial.pixel_to_world(1, 0)
<SkyCoord (ICRS): (ra, dec) in deg
    (205.43985907, 27.00341788)>
>>> wcs1 = cubeviz_helper.app.data_collection[1].coords
>>> wcs1.celestial.pixel_to_world(1, 0)
<SkyCoord (ICRS): (ra, dec) in deg
    (205.4398996, 27.00345399)>
>>> wcs1.celestial.pixel_to_world(0, 1)
<SkyCoord (ICRS): (ra, dec) in deg
    (205.43985907, 27.00341788)>

Maybe the proper fix is to propagate .meta['_orig_wcs'] consistently also to uncert and mask data in the parser and not worry about .celestial at all.

pllim · 2023-02-10T23:00:42Z

When I change the image_cube_hdu_obj_microns above to write out np.arange(720).reshape((8, 9, 10)), save the HDUList out to a FITS file, read the flux cube back into Ginga, then go to slice=1 (second slice), this is what I see in Ginga for the bottom left pixel that is x=0 y=0 in Cubeviz but x=1 y=1 in Ginga:

Sky coordinates agree more or less with flux cube coordinates display, but is a bit off for uncert.

pllim · 2023-02-10T23:34:37Z

This is what I don't understand... If you read that file made in #1991 (comment) back in:

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

filename = "image_cube_hdu_obj_microns.fits"
pf = fits.open(filename)
w = WCS(pf[1].header, pf)
flux = pf[1].data * u.Unit(pf[1].header['BUNIT'])

# UserWarning: Input WCS indicates that the spectral axis is not last.
# Reshaping arrays to put spectral axis last.
sp = Spectrum1D(flux=flux, wcs=w)

This is the original WCS (w) where flux.shape is (8, 9, 10):

WCS Keywords

Number of WCS axes: 3
CTYPE : 'RA---TAN'  'DEC--TAN'  'WAVE'  
CRVAL : 205.4384  27.004754  4.890499866509344e-06  
CRPIX : 38.0  38.0  1.0  
PC1_1 PC1_2 PC1_3  : -1.0  0.0  0.0  
PC2_1 PC2_2 PC2_3  : 0.0  1.0  0.0  
PC3_1 PC3_2 PC3_3  : 0.0  0.0  1.0  
CDELT : 3.61111097865634e-05  3.61111097865634e-05  1.000000047497451e-09  
NAXIS : 10  9  8

This is the WCS from Spectrum1D (sp.wcs) that supposedly goes with sp (sp.shape is (10, 9, 8)):

WCS Keywords

Number of WCS axes: 3
CTYPE : 'WAVE'  'DEC--TAN'  'RA---TAN'  
CRVAL : 4.890499866509344e-06  27.004754  205.4384  
CRPIX : 1.0  38.0  38.0  
PC1_1 PC1_2 PC1_3  : 1.0  0.0  0.0  
PC2_1 PC2_2 PC2_3  : 0.0  1.0  0.0  
PC3_1 PC3_2 PC3_3  : 0.0  0.0  -1.0  
CDELT : 1.000000047497451e-09  3.61111097865634e-05  3.61111097865634e-05  
NAXIS : 8  9  10

So when specutils say "reshaping arrays to put spectral axis last," I would expect the order is now (x, y, z) where z is the spectral axis. So let's say we try to sample the WCS at the lower left pixel of the first slice as per #1991 (comment) above, why is the original WCS (the one that specutils didn't touch) agree with Ginga but not the one that is supposed to go with the actual object? Did I misunderstand what "put spectral axis last" mean?

>>> w.pixel_to_world(0, 0, 1)[0].to_string('hmsdms')
'13h41m45.57590354s +27d00m12.3043716s'
>>> sp.wcs.pixel_to_world(0, 0, 1)[-1].to_string('hmsdms')
'13h41m45.56617642s +27d00m12.30437312s'

pllim · 2023-02-11T02:34:02Z

I think I understand now, the .meta["_orig_wcs"] was correct because it didn't get swap and I was assuming the opposite axes order, so it was a case of two wrongs make it right. See #2009 for hopefully a better fix.

stscijgbot-jwql · 2024-07-09T18:43:08Z

This issue is tracked on JIRA as JDAT-3070.

stscijgbot-jwql · 2024-07-09T18:43:09Z

Comment by Pey-Lian Lim on JIRA:

I have a proposed solution at #2009

pllim added bug Something isn't working cubeviz labels Feb 6, 2023

pllim mentioned this issue Feb 7, 2023

refactor mouseover coords info display #1976

Merged

9 tasks

pllim mentioned this issue Feb 9, 2023

Fix Cubeviz coordinates panel with inaccurate sky coordinates sometimes #2002

Closed

13 tasks

This was referenced Feb 10, 2023

local variable 'intercept' referenced before assignment ejeschke/ginga#1041

Closed

Fix Cubeviz mouseover bugs, add WCS to moment map #2009

Merged

pllim closed this as completed in #2009 Feb 27, 2023

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

BUG: Cubeviz parser gives wrong sky coordinates for non-flux cube under certain condition #1991

BUG: Cubeviz parser gives wrong sky coordinates for non-flux cube under certain condition #1991

pllim commented Feb 6, 2023 •

edited

Loading

pllim commented Feb 7, 2023

pllim commented Feb 9, 2023

pllim commented Feb 10, 2023 •

edited

Loading

pllim commented Feb 10, 2023 •

edited

Loading

pllim commented Feb 11, 2023

stscijgbot-jwql commented Jul 9, 2024

stscijgbot-jwql commented Jul 9, 2024

BUG: Cubeviz parser gives wrong sky coordinates for non-flux cube under certain condition #1991

BUG: Cubeviz parser gives wrong sky coordinates for non-flux cube under certain condition #1991

Comments

pllim commented Feb 6, 2023 • edited Loading

pllim commented Feb 7, 2023

pllim commented Feb 9, 2023

pllim commented Feb 10, 2023 • edited Loading

pllim commented Feb 10, 2023 • edited Loading

pllim commented Feb 11, 2023

stscijgbot-jwql commented Jul 9, 2024

stscijgbot-jwql commented Jul 9, 2024

pllim commented Feb 6, 2023 •

edited

Loading

pllim commented Feb 10, 2023 •

edited

Loading

pllim commented Feb 10, 2023 •

edited

Loading