APE 14: pixel_to_world_values ignores CUNIT

[jdavies-st]

Description

When reading in spectral coordinates, the WCS.pixel_to_world_values() method reports the value in different units than are specified in CUNIT for that axis. An example using JWST MIRI MRS data:

In [35]: from astroquery.mast import Observations

In [36]: uri = "mast:JWST/product/jw01864-c1006_t005_miri_ch2-short_s3d.fits"

In [37]: Observations.download_file(uri)
Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:JWST/product/jw01864-c1006_t005_miri_ch2-short_s3d.fits to /Users/jdavies/data/miri_mrs/jw01864-c1006_t005_miri_ch2-short_s3d.fits ...
|============================================================================================================|  28M/ 28M (100.00%)         3s
Out[37]: ('COMPLETE', None, None)

In [39]: from astropy.io import fits

In [40]: from astropy.wcs import WCS

In [41]: with fits.open("jw01864-c1006_t005_miri_ch4-short_s3d.fits", mode="readonly") as hdulist:
    ...:     wcs = WCS(hdulist['SCI'], hdulist)
    ...: 
WARNING: FITSFixedWarning: 'datfix' made the change 'Set DATE-BEG to '2022-12-31T03:46:17.520' from MJD-BEG.
Set DATE-AVG to '2022-12-31T04:16:09.295' from MJD-AVG.
Set DATE-END to '2022-12-31T04:46:01.102' from MJD-END'. [astropy.wcs.wcs]
2024-05-09 23:58:29,193 - stpipe - WARNING - FITSFixedWarning: 'datfix' made the change 'Set DATE-BEG to '2022-12-31T03:46:17.520' from MJD-BEG.
Set DATE-AVG to '2022-12-31T04:16:09.295' from MJD-AVG.
Set DATE-END to '2022-12-31T04:46:01.102' from MJD-END'.
WARNING: FITSFixedWarning: 'obsfix' made the change 'Set OBSGEO-L to   100.587422 from OBSGEO-[XYZ].
Set OBSGEO-B to     7.683995 from OBSGEO-[XYZ].
Set OBSGEO-H to 1691018989.934 from OBSGEO-[XYZ]'. [astropy.wcs.wcs]
2024-05-09 23:58:29,195 - stpipe - WARNING - FITSFixedWarning: 'obsfix' made the change 'Set OBSGEO-L to   100.587422 from OBSGEO-[XYZ].
Set OBSGEO-B to     7.683995 from OBSGEO-[XYZ].
Set OBSGEO-H to 1691018989.934 from OBSGEO-[XYZ]'.

In [42]: wcs.pixel_to_world(12, 12, 6)
Out[42]: 
[<SkyCoord (ICRS): (ra, dec) in deg
     (56.67226241, -52.08471111)>,
 <SpectralCoord 
    (observer: <ICRS Coordinate: (ra, dec, distance) in (deg, deg, m)
                   (100.70828062, 23.42770565, 1.45953725e+11)
                (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
                   (0., 0., 0.)>
     target: <ICRS Coordinate: (ra, dec, distance) in (deg, deg, kpc)
                 (56.67147134, -52.08422501, 1000.)
              (pm_ra_cosdec, pm_dec, radial_velocity) in (mas / yr, mas / yr, km / s)
                 (0., 0., 0.)>
     observer to target (computed from above):
       radial_velocity=0.0 km / s
       redshift=0.0)
   1.77390008e-05 m>]

In [43]: wcs.pixel_to_world_values(12, 12, 6)
Out[43]: (array(56.67226241), array(-52.08471111), array(1.77390008e-05))

In [45]: print(hdulist['SCI'].header['CUNIT*'])
CUNIT1  = 'deg     '           / Axis 1 units                                   
CUNIT2  = 'deg     '           / Axis 2 units                                   
CUNIT3  = 'um      '           / Axis 3 units                                   

As you can see, CUNIT3 is in microns, but the value returned is in meters.

Expected behavior

I would expect that the pixel_to_world_values would report values in units specified by CUNITn. This is what gwcs.WCS does for the same file. First, pip install jwst, then:

In [46]: from jwst import datamodels

In [47]: with datamodels.open("jw01864-c1006_t005_miri_ch2-short_s3d.fits") as s3d:
    ...:     wcs = s3d.meta.wcs
    ...: 

In [48]: wcs.pixel_to_world(12, 12, 6)
Out[48]: 
[<SkyCoord (ICRS): (ra, dec) in deg
     (56.67211301, -52.08454364)>,
 <SpectralCoord 7.51845023 um>]

In [49]: wcs.pixel_to_world_values(12, 12, 6)
Out[49]: (56.67211301091073, -52.08454363698079, 7.51845022890484)

Here it is returned in microns, as expected.

How to Reproduce

See above.

Versions

macOS-10.16-x86_64-i386-64bit
Python 3.11.9 (main, Apr 19 2024, 11:44:45) [Clang 14.0.6 ]
astropy 6.1.0
Numpy 1.26.4
pyerfa 2.0.1.4
Scipy 1.13.0
Matplotlib 3.8.4
gwcs 0.21.0

[astrofrog]

The WCS class converts units to SI (this is actually done by WCSLIB behind the scenes). If you try and print out the WCS you will probably see it is in m. This is also a pet peeve of @Cadair's as WCS converts arcsec to degrees. This isn't technically a bug as we know about it and it is a limitation of WCSLIB.

That's not to say we can't try and solve this though - perhaps we could save the original units in a private attribute in the WCS class and then expose the original units via the APE 14 API and deal with conversions behind the scenes.

[Cadair]

Yes this is super frustrating.

It would also be nice if there was a way to undo this conversion on the result of WCS.to_header().

[jdavies-st]

That's not to say we can't try and solve this though - perhaps we could save the original units in a private attribute in the WCS class and then expose the original units via the APE 14 API and deal with conversions behind the scenes.

Agree. Given that the APE 14 is built on top of WCSLIB, it should be trivial to report the correct units and do the conversion in the wrappers, and more importantly give consistent behavior between astropy.wcs and gwcs. Does APE 14 have an opinion on units? Independent of whatever WCSLIB does?

If my header defines CUNITn, I would expect to get the same CUNITn back with the APE 14 interface. Getting back a unit other than what is in CUNITn has to be a bug, as there's no way to predict it. If my unit is optical depth in a Parma cheese rind, I don't expect that to get that value converted to an SI unit. I expect to get back units of cheese rind depth.

[Cadair]

The only question is do we consider this a breaking API change or a bugfix?

[astrofrog]

Perhaps we should first try and implement it as opt-in via a kwarg and then decide later if we want to make it the default via a deprecation period? Having it opt-in would allow us to iron out any bugs anyway.

[pllim]

Are these issues related?

• WCS updating units should emit warning (CUNITn is not respected) #3658
• WCS.world_to_array_index_values ignores units of inputs #11173

[astrofrog]

Yep!

[pllim]

Are they duplicate issues?

[jdavies-st]

Not duplicate issues, as each issue shows how different methods suffer from the similar issues of internal conversion to specific units and in the case of #3658, how the FITS headers describing the WCS don't round-trip properly because of this internal conversion.

I would argue that APE 14 is a new API, and it specifically uses quantities, so it should respect quantities and the quantity values explicitly. It does not, which is a bug.

Agree that an implementation first would be nice, and then we can decide whether its a bug or a new feature (and for which methods), and whether it just gets implemented for the APE 14 interface or for all the methods.

[astrofrog]

One slight complication is that e.g. to_header is not part of the APE 14 API, so we might need to extend this behavior beyond just the APE 14 methods. It's also going to be a bit confusing if the APE 14 API disagrees unit-wise with e.g. wcs.wcs.cunit.

I have reached out to Mark to get his input on the matter before we proceed with anything on our side (in case it's possible for him to add an option to not convert).

Rather than just change the APE 14 API, there are two alternatives I wanted to mention:

• Create a Wcsprm wrapper which preserves the units and auto-converts, and we can then optionally use that inside WCS if a user sets e.g. preserve_units=True in the WCS constructor. This has the benefit that everything user-facing is consistent unit-wise.

• Create a new alternative high-level FITSWCS class which would provide only the APE 14 API as well as Python properties at the top level of the class to set the various WCS parameters and also to_header. This would be a chance to break away from the current API which is a bit complex especially with the public nested .wcs attribute, and make something new and very Pythonic.

[Cadair]

I am pro a new class, though it's probably a lot of work.

[jdavies-st]

Yes, a lot of work.

Perhaps we should first try and implement it as opt-in via a kwarg and then decide later if we want to make it the default via a deprecation period? Having it opt-in would allow us to iron out any bugs anyway.

Agree this is the way to go for now. And do this for the APE 14 interface, as at least with that interface, people are already thinking about units and want units to round trip. And there's less legacy code using this interface, so deprecations should have lower impact.

The opt-in behavior might be something like a preserve_input_units=True kwarg that modifies several methods from both BaseLowLevelWCS and BaseHighLevelWCS such as pixel_to_world, pixel_to_world_values, world_to_pixel, world_to_pixel_values, world_axis_units, plus the array_* versions of these. The new behavior would utilize the units in CUNITn, converting on the way in and out as needed when passing to the underlying WCSLIB machinery.

FWIW, here's another manifestation of the same problem from my example above:

In [59]: hdulist['SCI'].header['CUNIT3']
Out[59]: 'um'

In [60]: wcs.world_axis_units
Out[60]: ['deg', 'deg', 'm']

Again, that's the APE 14 interface.

One could conceivable extend this to WCS.to_header() as well, as this already wraps WCS.wcs.to_header(), which for some reason is also part of the public interface. But anyway, one could add the same preserve_input_units kwarg here as well, though I see this as a separate PR perhaps?