astropy: Presence of SIP keywords leads to ignored PV keywords.

Description

I am working on updating the fits header for one of our telescopes. We wanted to represent the distortions in SIP convention and the projection to be ‘CAR’. While working on this, I noticed if SIP coefficients are present in the header and/or ‘-SIP’ is added to CTYPEia keywords, astropy treats the PV keywords (PV1_0, PV1_1 and PV1_2) as “redundant SCAMP distortions”.

Earlier I could not figure out why the projection weren’t going as I expected, and I suspected a bug in astropy wcs. After talking to Mark Calabretta about the difficulties I’m facing, that suspicion only grew stronger. The header was working as expected with WCSLIB but was giving unexpected behavior in astropy wcs.

The following would be one example header -

header_dict = {
'SIMPLE'  : True, 
'BITPIX'  : -32, 
'NAXIS'   :  2,
'NAXIS1'  : 1024,
'NAXIS2'  : 1024,
'CRPIX1'  : 512.0,
'CRPIX2'  : 512.0,
'CDELT1'  : 0.01,
'CDELT2'  : 0.01,
'CRVAL1'  : 120.0,
'CRVAL2'  : 29.0,
'CTYPE1'  : 'RA---CAR-SIP',
'CTYPE2'  : 'DEC--CAR-SIP',
'PV1_1'   :120.0,
'PV1_2'   :29.0,
'PV1_0'   :1.0,
'A_ORDER' :2,
'A_2_0'   :5.0e-4,
'B_ORDER' :2,
'B_2_0'   :5.0e-4
}
from astropy.io import fits
header = fits.Header(header_dict)

Expected behavior

When you parse the wcs information from this header, the image should be centered at ra=120 and dec=29 with lines of constant ra and dec looking like the following image generated using wcslib -

Actual behavior

If I parse the wcs information using astropy wcs, it throws the following warning - WARNING: FITSFixedWarning: Removed redundant SCAMP distortion parameters because SIP parameters are also present [astropy.wcs.wcs] And the resulting grid is different. Code -

import numpy as np
import matplotlib.pyplot as plt
from astropy.wcs import WCS
w = WCS(header)
ra = np.linspace(116, 126, 25)
dec = np.linspace(25, 34, 25)

for r in ra:
    x, y = w.all_world2pix(np.full_like(dec, r), dec, 0)
    plt.plot(x, y, 'C0')
for d in dec:
    x, y = w.all_world2pix(ra, np.full_like(ra, d), 0)
    plt.plot(x, y, 'C0')

plt.title('Lines of constant equatorial coordinates in pixel space')
plt.xlabel('x')
plt.ylabel('y')

Grid -

The astropy wcs grid/solution does not change whethere we keep or remove the PV keywords. Furthermore, the astropy grid can be recreated in wcslib, by removing the PV keywords.

Steps to Reproduce

1. Initialize the header
2. Parse the header using astropy.wcs.WCS
3. Plot the graticule
4. Remove the PV keywords and run again
5. You will find the same graticule indicating that PV keywords are completely ignored.
6. Additionally, the graticules can be compared with the wcsgrid utility of wcslib.

System Details

Linux-6.0.11-200.fc36.x86_64-x86_64-with-glibc2.35 Python 3.9.12 (main, Apr 5 2022, 06:56:58) [GCC 7.5.0] Numpy 1.21.5 pyerfa 2.0.0 astropy 5.1 Scipy 1.7.3 Matplotlib 3.5.1