Repository navigation
APE 14: pixel_to_world_values ignores CUNIT #16432
Description
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
Activity
astrofrog commentedon May 9, 2024
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 commentedon May 10, 2024
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 commentedon May 10, 2024
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.wcsandgwcs. 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.
[-]pixel_to_world_values ignores CUNIT[/-][+]APE 14: pixel_to_world_values ignores CUNIT[/+]Cadair commentedon May 10, 2024
The only question is do we consider this a breaking API change or a bugfix?
astrofrog commentedon May 10, 2024
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 commentedon May 10, 2024
Are these issues related?
astrofrog commentedon May 10, 2024
Yep!
pllim commentedon May 10, 2024
Are they duplicate issues?
jdavies-st commentedon May 13, 2024
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 commentedon May 13, 2024
One slight complication is that e.g.
to_headeris 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
WCSif a user sets e.g.preserve_units=Truein theWCSconstructor. This has the benefit that everything user-facing is consistent unit-wise.Create a new alternative high-level
FITSWCSclass 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 alsoto_header. This would be a chance to break away from the current API which is a bit complex especially with the public nested.wcsattribute, and make something new and very Pythonic.Cadair commentedon May 13, 2024
I am pro a new class, though it's probably a lot of work.
jdavies-st commentedon May 13, 2024
Yes, a lot of work.
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=Truekwarg that modifies several methods from bothBaseLowLevelWCSandBaseHighLevelWCSsuch aspixel_to_world,pixel_to_world_values,world_to_pixel,world_to_pixel_values,world_axis_units, plus thearray_*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:
Again, that's the APE 14 interface.
One could conceivable extend this to
WCS.to_header()as well, as this already wrapsWCS.wcs.to_header(), which for some reason is also part of the public interface. But anyway, one could add the samepreserve_input_unitskwarg here as well, though I see this as a separate PR perhaps?