Repository navigation
WCS updating units should emit warning (CUNITn is not respected) #3658
Description
In the following example, a warning should really be emitted indicating that the units are being changed:
In [14]: import numpy as np
In [15]: from astropy import wcs
In [16]: w2 = wcs.WCS(naxis=2)
In [17]: w2.wcs.crpix = np.array([ 2060.389893, 2037.25 ])
In [18]: w2.wcs.cdelt = np.array([ 0.600109, 0.600109])
In [19]: w2.wcs.crval = np.array([ 0., 0.])
In [20]: w2.wcs.ctype = ['HPLN-TAN', 'HPLT-TAN']
In [21]: w2.wcs.cunit = ['arcsec', 'arcsec']
In [22]: w2.wcs.cunit
Out[22]: ['arcsec', 'arcsec']
In [23]: w2.wcs_pix2world([[2000,2000]], 0)
Out[23]: array([[ 3.59990100e+02, -6.04276412e-03]])
In [24]: w2.wcs.cunit
Out[24]: ['deg', 'deg']
This is based on an example that @Cadair showed me.
cc @mdboom
Activity
Cadair commentedon Apr 2, 2015
Interesting that it only updates post the
wcs_pix2worldcall.mdboom commentedon Apr 2, 2015
Everything gets normalized before computation -- so it's converting the arcsec to deg internally (and adding the appropriate scale). I'm not sure if it's feasible to warn every time wcslib normalizes the data structure (it pretty much does that as a matter of course for the majority of WCSes). This is essentially a FAQ at this point --
astropy.wcsshould be used to create or parse WCS data structures, but it is not suitable for editing FITS headers directly where you need to specifically control they key/value pairs (there are similar issues for CD vs PC normalization and a number of other things). For that, one should useastropy.io.fits. Probably need a big red note in the docs about that.Part of the problem, perhaps, here is that unlike the regular C API which encourages
setto be called before computation,astropy.wcscalls it implicitly on the user's behalf.nden commentedon Apr 2, 2015
Is the issue here that the input header is in units of 'arcsec' but since wcslib works internally with degrees, it returns the result in degrees? I've ran into this before. And I agree that it is unexpected (I would expect to get back a result in arcsec.)
I consider this one more case of "FITS knows better than you what you want".
mdboom commentedon Apr 2, 2015
Yeah -- this is just was wcslib does: It normalizes all units before computation. It's fairly reasonable if you want to load a WCS and use it. It's not so helpful if you're trying to control what ends up in the FITS header.
astrofrog commentedon Apr 2, 2015
I still feel like the call to set should warn about the unit change. I'm familiar with the normalization behavior, and we're not going to change it now given all the code that relies on this, but I think that in the specific case where one creates a WCS programmatically and explicitly sets the units to a certain pair of values, the fact that the call to wcs_pix2world results in
setbeing called and changes the values should definitely emit a warning when the units are changed.astropy.io.fitscurrently emits warnings for non-standard things if I'm not mistaken, so I guess I don't see the downside of emitting warnings for things that are changed during the WCS normalization? E.g.This would help avoid bad surprises by users.
astrofrog commentedon Apr 2, 2015
As a side note, I was thinking that it would be cool to have a
WCSclass that understands and also returnsQuantityobjects. Call itQWCS(likeQTable). This is different from previous suggestions to make it understandSkyCoord. Is there any reason we could do theQWCSsubclass?mdboom commentedon Apr 2, 2015
But how do we detect that specific case from loading from a header? I'm not sure there's a robust way...
QWCSsounds like a good idea. Does it need to be a separate class? If you pass in pixel quantities, you get degree quantities out... pass in raw arrays (as now) get raw arrays out...pllim commentedon Mar 31, 2022
Ugh... I think I just ran into this. I saw a
CUNIT3 = 'Angstrom'but it keeps giving me backwcs.wcs.cunit = ['deg', 'deg', 'm']. I thought I was going crazy until @keflavich enlightened me.https://github.com/keflavich/Orion_ALMA_2016.1.00165.S/blob/2ffef30d3b0eb7e7bb8a1c7a36170386a2974fbd/analysis/kinematic_analysis_pv.py#L202-L203
Can we turn off unit conversion altogether and respect the
CUNITunit?keflavich commentedon Mar 31, 2022
I'll also point to this:
https://github.com/radio-astro-tools/spectral-cube/blob/7076f7376d687c6334f2485242227a28c4856f55/spectral_cube/spectral_axis.py#L200
as the machinery we had to build in spectral-cube to get spectral axes in the units we wanted.
[-]WCS updating units should emit warning[/-][+]WCS updating units should emit warning (CUNITn is not respected)[/+]