astropy: Do not automatically apply SIP correction to alternative WCSes if SIP suffix is missing

Description

If SIP coefficients are present in the header, Astropy applies it to all WCSes. If CTYPE of the WCS in hand doesn’t have a -SIP suffix, Astropy will generate a warning and apply the SIP correction. It may be a valid approach to apply SIP to the primary WCS regardless of -SIP in CTYPE, however, Astropy should not automatically apply SIP correction to alternative WCSes if they don’t have the -SIP ending in their CTYPE and it should not produce any warnings. This warning and default behavior has been a source of confiusion in the past (e.g. issue 9239).

Steps to Reproduce

1. The following script will generate a minimalist FITS file to reproduce the described behavior in the next step:

#! /usr/bin/env python
from astropy.io import fits
import numpy as np
from pathlib import Path

produtsPath = Path().parent / 'products'
produtsPath.mkdir(exist_ok=True)


def get_image_data(npix=99):
    y = np.arange(1, 1 + npix, dtype=np.float32)
    x = y / 100
    x_mesh, y_mesh = np.meshgrid(x, y)
    return x_mesh + y_mesh


def construct_sample_hdul(npix=99):

    data = get_image_data(npix=npix)

    hdul = fits.HDUList([fits.PrimaryHDU()])
    hdul.append(fits.ImageHDU(data=data, name='DATA'))
    hdr = hdul[1].header

    hdr['WCSAXES'] = 2

    hdr['CTYPE1 '] = 'RA---SIN-SIP'
    hdr['CTYPE2 '] = 'DEC--SIN-SIP'

    hdr['CUNIT1 '] = 'deg'
    hdr['CUNIT2 '] = 'deg'

    hdr['CRPIX1 '] = 1  # 1.5 + npix / 2
    hdr['CRPIX2 '] = 1  # 1.5 + npix / 2


    hdr['CRVAL1 '] = 0
    hdr['CRVAL2 '] = 90

    hdr['CDELT1 '] = round(6.2 / 3600, 5)
    hdr['CDELT2 '] = round(6.2 / 3600, 5)

    # i_j: output_input
    hdr['PC1_1  '] = 1
    hdr['PC1_2  '] = 0
    hdr['PC2_1  '] = 0
    hdr['PC2_2  '] = 1

    # ----------------
    # SIP coefficients
    a = 1e-5
    hdr['A_ORDER'] = 2
    hdr['A_0_2  '] = a
    hdr['A_1_1  '] = a
    hdr['A_2_0  '] = a

    hdr['B_ORDER'] = 2
    hdr['B_0_2  '] = a
    hdr['B_1_1  '] = a
    hdr['B_2_0  '] = a

    # Inverse SIP coefficients
    hdr['AP_ORDER'] = 2
    hdr['AP_0_2  '] = a
    hdr['AP_1_1  '] = a
    hdr['AP_2_0  '] = a

    hdr['BP_ORDER'] = 2
    hdr['BP_0_2  '] = a
    hdr['BP_1_1  '] = a
    hdr['BP_2_0  '] = a

    # ----------------
    hdr['WCSAXESW'] = 2

    hdr['CTYPE1W '] = 'WAVE'
    hdr['CTYPE2W '] = 'WAVE'

    hdr['CUNIT1W '] = 'm'
    hdr['CUNIT2W '] = 'm'

    hdr['CRPIX1W '] = 0
    hdr['CRPIX2W '] = 0

    hdr['CRVAL1W '] = 0
    hdr['CRVAL2W '] = 0

    hdr['CDELT1W '] = 1
    hdr['CDELT2W '] = 1

    # i_j: output_input
    hdr['PC1_1W  '] = 1
    hdr['PC1_2W  '] = 0
    hdr['PC2_1W  '] = 0
    hdr['PC2_2W  '] = 1
    return hdul


if __name__ == '__main__':
    hdul = construct_sample_hdul()
    hdul.writeto(produtsPath / 'minimal.fits', overwrite=True)

1. Initialize and read file header:

# Standar Dependencies
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
import warnings
from pathlib import Path

# Notebook
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

# Read File
detector_id = 2
produtsPath = Path().parent / 'products'
# filePath = produtsPath / f'realistic_lvf_wcs_2p2d_array-{detector_id}.fits'
filePath = produtsPath / f'minimal.fits'
hdul = fits.open(filePath)

hdr = hdul[1].header

1. Reconstruct the alternative WCS without suppressing SIP. To suppress, uncomment wcs.sip = None.

# Create a wcs object
with warnings.catch_warnings():
    warnings.filterwarnings(
        'ignore', message='.*The WCS transformation has more axes.*')
    wcs = WCS(hdr, hdul, 'W')
    # !!! toggle next line
    # wcs.sip = None

1. Depending on whether SIP was manually suppressed or not, the following pixel to world transformations would be different.

# Query wavelength and bandwidth for a set of pixel coordinates
if True:
    b = 99
    x1 = [0, b, b, 0]
    y1 = [0, 0, b, b]
    corner_wl_out, corner_bandwidth_out = wcs.pixel_to_world(x1, y1)

    # Generate a report
    print(f'\nFile was read: {filePath.stem}')
    print('Input x, y Pixel Coord (counting from 0, 0):')
    print(x1, y1)
    print('Output Wavelength:')
    print(corner_wl_out)
    print('Output Bandwidth:')
    print(corner_bandwidth_out)
    print('\n------------------------------------\n')

Output with SIP suppression:

File was read: minimal
Input x, y Pixel Coord (counting from 0, 0):
[0, 99, 99, 0] [0, 0, 99, 99]
Output Wavelength:
[  1. 100. 100.   1.] m
Output Bandwidth:
[  1.   1. 100. 100.] m

Output without SIP suppression:

File was read: minimal
Input x, y Pixel Coord (counting from 0, 0):
[0, 99, 99, 0] [0, 0, 99, 99]
Output Wavelength:
[  1.00003 100.10101 100.3       1.10101] m
Output Bandwidth:
[  1.00003   1.10101 100.3     100.10101] m

This warning is generated regardless:

INFO: 
                Inconsistent SIP distortion information is present in the FITS header and the WCS object:
                SIP coefficients were detected, but CTYPE is missing a "-SIP" suffix.
                astropy.wcs is using the SIP distortion coefficients,
                therefore the coordinates calculated here might be incorrect.

                If you do not want to apply the SIP distortion coefficients,
                please remove the SIP coefficients from the FITS header or the
                WCS object.  As an example, if the image is already distortion-corrected
                (e.g., drizzled) then distortion components should not apply and the SIP
                coefficients should be removed.

                While the SIP distortion coefficients are being applied here, if that was indeed the intent,
                for consistency please append "-SIP" to the CTYPE in the FITS header or the WCS object.

                 [astropy.wcs.wcs]

Desired Behavior

No warning, no SIP correction as default behavior for alternative WCSes if CTYPE is missing a -SIP suffix.

System Details

Linux-5.16.15-76051615-generic-x86_64-with-glibc2.34
Python 3.9.7 | packaged by conda-forge | (default, Sep  2 2021, 17:58:34) 
[GCC 9.4.0]
Numpy 1.22.3
pyerfa 2.0.0.1
astropy 5.0.4
Scipy 1.8.0
Matplotlib 3.5.1