astropy: Incompatible results from WCS.pixel_to_world between astropy 3.x and 4.x (distortion in WCSLIB 6.3 + Cutout2D)

EDIT: The cause identified in https://github.com/astropy/astropy/issues/11216#issuecomment-755754911

Description

The use of WCS.pixel_to_world gives different results when using astropy 3.x or 4.x (same environment).

Expected behavior

The two astropy versions should give the same results, or expected incompatibilities should be quoted explicitly in the documentation.

Actual behavior

The two astropy versions give different results, and I could not find a documentation on it.

Steps to Reproduce

The common.fits can be found here (37 KB).

from astropy.io import fits
from astropy.wcs import WCS

fn = "common.fits"

header = fits.open(fn)[0].header
w = WCS(header)

print(w.pixel_to_world([[0, 0]], 0))
# astropy 3.2.2: 
# <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg
#   [[(250.51824533, -75.47419369), (250.51824533, -75.47419369)]]>
# astropy 4.0.1.post1:
# <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg
#   [[(251.99786293, -75.83653121), (251.99786293, -75.83653121)]]>

I am not sure what’s happening here. I have the impression the information from the header is used differently in the two astropy versions, but I cannot find what has changed explicitly in the documentation.

System Details

Darwin-20.2.0-x86_64-i386-64bit
Python 3.7.1 (default, Oct 23 2018, 14:07:42)
[Clang 4.0.1 (tags/RELEASE_401/final)]
Numpy 1.16.2
astropy 3.2.2
Scipy 1.4.1
Matplotlib 3.0.1

Darwin-20.2.0-x86_64-i386-64bit
Python 3.7.1 (default, Oct 23 2018, 14:07:42)
[Clang 4.0.1 (tags/RELEASE_401/final)]
Numpy 1.16.2
astropy 4.0.1.post1
Scipy 1.4.1
Matplotlib 3.0.1

Checked also using the official docker images (https://hub.docker.com/u/continuumio/): 2019.10 (astropy 3.2.2), 2020.02 (astropy 4.0), 2020.07 (astropy 4.0.1.post1), 2020.11 (astropy 4.0.2).

About this issue

  • Original URL
  • State: closed
  • Created 3 years ago
  • Comments: 54 (38 by maintainers)

Most upvoted comments

I note that the distortion-free parts of the WCS (using the header without the CQDIS/DQ cards) return identical results: 3.2.2

>>> w = wcs.WCS(fits.Header.fromstring(hdr[:1802], sep='\n'))
>>> wd = wcs.WCS(fits.Header.fromstring(hdr, sep='\n'))
>>> w.pixel_to_world(0, 0)
<SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg
    (250.51820607, -75.47417884)>
>>> wd.pixel_to_world(0, 0).separation(w.pixel_to_world(0, 0))
<Angle 1.78130314e-05 deg>

4.2

>>> w = wcs.WCS(fits.Header.fromstring(hdr[:1802], sep='\n'))
>>> wd = wcs.WCS(fits.Header.fromstring(hdr, sep='\n'))
>>> w.pixel_to_world(0, 0)
<SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg
    (250.51820607, -75.47417884)>
>>> wd.pixel_to_world(0, 0).separation(w.pixel_to_world(0, 0))
<Angle 0.51542227 deg>
>>> w.pixel_to_world(736.32251317135, 806.0426298396).separation(w.pixel_to_world(0, 0))
<Angle 0.51548823 deg>
>>> wd.pixel_to_world(736.32251317135, 806.0426298396).separation(wd.pixel_to_world(0, 0))
<Angle 1.0309105 deg>

IOW in the 4.x version the distortion seems to offset the point by almost exactly the full distance to CRPIX1, CRPIX2 again (instead of ~64 mas, blowing up the field by a factor of 2)! Unless this is an almost pathologically distorted WCS (does not look like it, with diagonal coefficients close to 1.0 AFAICS), this looks like a serious bug.

+2
dhomeier on Jan 6, 2021

Hi @pllim – thanks for the follow-up. I did not get more into the details since last year, and did not encounter this problem again in my work.

+1
JulienPeloton on Jun 3, 2022

Interestingly enough, I cannot find an open issue here to indicate any plans of adding distortion support in Cutout2D…

+1
pllim on Jan 7, 2021

From what I understand,

  1. WCSLIB 6.3 does not bring inconsistencies, but it changed how distortion is handled.
  2. Cutout2D does not handle cases where the input WCS object contains distortion lookup tables
  3. So cutouts created from files with distortion gets affected by the change of WCSLIB versions.

It is not a genuine bug as it is mentioned in the doc (I would call it a limitation perhaps), and this should be resolved the day Cutout2D gets updated to include the distortion part (if this is the plan one day).

+1
JulienPeloton on Jan 7, 2021

Also,

In [56]: w.wcs.p2s([[2, 7]], 1)['world']
Out[56]: array([[251.99254424, -75.8309369 ]])

In [57]: w.wcs_pix2world([[2, 7]], 1)
Out[57]: array([[251.99254424, -75.8309369 ]])

In [58]: w.all_pix2world([[2, 7]], 1)
Out[58]: array([[251.99254424, -75.8309369 ]])

There is no difference between all_pix2world and wcs_pix2world which means that all_* does nothing special to handle distortions: everything is handled internally to WCSLIB.

+1
mcara on Jan 7, 2021

I added a “Close?” label but I’ll keep it open for a bit more in case others have opinions on how to improve the user experience regarding WCSLIB upgrades. 😸

Thanks for your patience and understanding!

+1
pllim on Jan 6, 2021

Thanks! I think one can close the issue then.

+1
JulienPeloton on Jan 6, 2021

Actually, that info is available but it is kinda hidden. 😬

>>> from astropy.wcs import _wcs
>>> _wcs.__version__
'7.3.0'
+1
pllim on Jan 6, 2021

Thanks @dhomeier, I think you found where the issue is coming from (distortion).

By bisecting, I found the PR that introduced behavior changes: https://github.com/astropy/astropy/pull/9125: update wcslib to 6.4

Actually, this PR made an update from WCSLIB 6.2 to 6.4 directly. And looking at WCSLIB changelog, one can read for 6.3:

https://github.com/astropy/astropy/blob/4b242e6f5465410aa575f844e1a0b81db4534c94/cextern/wcslib/CHANGES#L160-L169

So the way distortions are treated changed (for good or bad - that I do not know). It is not a bug in astropy then, but this should probably be documented somewhere else than just in the changelog of WCSLIB.

+1
JulienPeloton on Jan 6, 2021

~At the very least, w.to_header() might give a clue.~ (Thanks for the info!)

Also, what about the dev version of astropy? 💭

+1
pllim on Jan 6, 2021

Thanks - which information in particular do you have in mind? here is the header:

SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                   16 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                  100
NAXIS2  =                  150
WCSAXES =                    2 / Number of coordinate axes
CRPIX1  =      737.32251317135 / Pixel coordinate of reference point
CRPIX2  =       807.0426298396 / Pixel coordinate of reference point
PC1_1   =    0.025278924159706 / Coordinate transformation matrix element
PC1_2   =  0.00037445145190975 / Coordinate transformation matrix element
PC2_1   = -0.00037626636323055 / Coordinate transformation matrix element
PC2_2   =    0.025284403492745 / Coordinate transformation matrix element
CDELT1  =    -0.01867555394435 / [deg] Coordinate increment at reference point
CDELT2  =     0.01867555394435 / [deg] Coordinate increment at reference point
CUNIT1  = 'deg'                / Units of coordinate increment and value
CUNIT2  = 'deg'                / Units of coordinate increment and value
CTYPE1  = 'RA---TAN'           / Right ascension, gnomonic projection
CTYPE2  = 'DEC--TAN'           / Declination, gnomonic projection
CRVAL1  =      249.10970770833 / [deg] Coordinate value at reference point
CRVAL2  =     -75.102948144444 / [deg] Coordinate value at reference point
LONPOLE =                180.0 / [deg] Native longitude of celestial pole
LATPOLE =     -75.102948144444 / [deg] Native latitude of celestial pole
WCSNAME = 'DSS PLATEID 0194'   / Coordinate system title
RADESYS = 'FK5'                / Equatorial coordinate system
EQUINOX =               2000.0 / [yr] Equinox of equatorial coordinates
MJD-OBS =              42867.0 / [d] MJD of observation matching DATE-OBS
DATE-OBS= '1976-03-30'         / ISO-8601 observation date matching MJD-OBS
CQDIS1  = 'TPD'                / Q = sequent, Template Polynomial Distortion
DQ1     = 'NAXES:  2'          / Two independent variables
DQ1     = 'AXIS.1: 1'          / 1st independent variable: axis 1 (= u)
DQ1     = 'AXIS.2: 2'          / 2nd independent variable: axis 2 (= v)
DQ1     = 'AUX.1.COEFF.0:   -0.030369554325293' / TPD: x = c0 + c1*u + c2*v
DQ1     = 'AUX.1.COEFF.1:    -0.99999816063804' / TPD: x = c0 + c1*u + c2*v
DQ1     = 'AUX.1.COEFF.2:    0.014809554960054' / TPD: x = c0 + c1*u + c2*v
DQ1     = 'AUX.2.COEFF.0:   0.0013505128079274' / TPD: y = d0 + d1*u + d2*v
DQ1     = 'AUX.2.COEFF.1:    0.014881334702971' / TPD: y = d0 + d1*u + d2*v
DQ1     = 'AUX.2.COEFF.2:     0.99978145301581' / TPD: y = d0 + d1*u + d2*v
DQ1     = 'TPD.FWD.0:    -0.030382917644437' / TPD coefficient: 1
DQ1     = 'TPD.FWD.1:     -0.99978145301581' / TPD coefficient: x
DQ1     = 'TPD.FWD.2:     0.014809554960054' / TPD coefficient: y
DQ1     = 'TPD.FWD.4:   -3.321209956147E-08' / TPD coefficient: xx
DQ1     = 'TPD.FWD.5:   2.1215435903186E-07' / TPD coefficient: xy
DQ1     = 'TPD.FWD.6:   2.2014764207626E-07' / TPD coefficient: yy
DQ1     = 'TPD.FWD.7:  -3.3806674482692E-08' / TPD coefficient: xxx
DQ1     = 'TPD.FWD.8:   9.4293173357228E-10' / TPD coefficient: xxy
DQ1     = 'TPD.FWD.9:    -3.54701586968E-08' / TPD coefficient: xyy
DQ1     = 'TPD.FWD.10:  1.2163819157713E-09' / TPD coefficient: yyy
CQDIS2  = 'TPD'                / Q = sequent, Template Polynomial Distortion
DQ2     = 'NAXES:  2'          / Two independent variables
DQ2     = 'AXIS.1: 1'          / 1st independent variable: axis 1 (= u)
DQ2     = 'AXIS.2: 2'          / 2nd independent variable: axis 2 (= v)
DQ2     = 'AUX.1.COEFF.0:   0.0013505128079274' / TPD: x = c0 + c1*u + c2*v
DQ2     = 'AUX.1.COEFF.1:    0.014881334702971' / TPD: x = c0 + c1*u + c2*v
DQ2     = 'AUX.1.COEFF.2:     0.99978145301581' / TPD: x = c0 + c1*u + c2*v
DQ2     = 'AUX.2.COEFF.0:   -0.030369554325293' / TPD: y = d0 + d1*u + d2*v
DQ2     = 'AUX.2.COEFF.1:    -0.99999816063804' / TPD: y = d0 + d1*u + d2*v
DQ2     = 'AUX.2.COEFF.2:    0.014809554960054' / TPD: y = d0 + d1*u + d2*v
DQ2     = 'TPD.FWD.0:  -0.00089857082115081' / TPD coefficient: 1
DQ2     = 'TPD.FWD.1:      0.99999816063804' / TPD coefficient: x
DQ2     = 'TPD.FWD.2:     0.014881334702971' / TPD coefficient: y
DQ2     = 'TPD.FWD.4:  -3.2511261465127E-07' / TPD coefficient: xx
DQ2     = 'TPD.FWD.5:   1.5287025105602E-07' / TPD coefficient: xy
DQ2     = 'TPD.FWD.6:  -2.2662497174755E-07' / TPD coefficient: yy
DQ2     = 'TPD.FWD.7:   3.5215497254629E-08' / TPD coefficient: xxx
DQ2     = 'TPD.FWD.8:   1.8773137889193E-09' / TPD coefficient: xxy
DQ2     = 'TPD.FWD.9:     3.35202384015E-08' / TPD coefficient: xyy
DQ2     = 'TPD.FWD.10: -2.5608500457886E-09' / TPD coefficient: yyy
+1
JulienPeloton on Jan 6, 2021
Read more comments on GitHub
← astropy: Importing astropy.coordinates fails with python -OO flag
astropy: Incorrect WCS conversion with SIP distortions (error ~1 arcmin!) →