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
- Initialize the header
- Parse the header using astropy.wcs.WCS
- Plot the graticule
- Remove the PV keywords and run again
- You will find the same graticule indicating that PV keywords are completely ignored.
- 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
About this issue
- Original URL
- State: closed
- Created a year ago
- Comments: 19 (10 by maintainers)
Just commenting out https://github.com/astropy/astropy/blob/966be9fedbf55c23ba685d9d8a5d49f06fa1223c/astropy/wcs/wcs.py#L793 solves the issue for me. But, I don’t know if that would be desirable as we might be back to square one with the old PTF images.
Once the appropriate approach for fixing this is decided, I can try to make a small PR.
PV
keywords are not optional keywords in CAR projection to relate the native spherical coordinates with celestial coordinates (RA, Dec). By default they have values equal to zero, but in my case I need to define these parameters.