Some notes on how to convert a FITS WCS to a WCS from gwcs (Generalized WCS) to be used within the JWST pipeline

import astropy
from astropy.io import fits
import numpy as np
from astropy import units as u
from astropy import wcs
hdulist = fits.open("example_field/iris_sim_gc_filterKN3.fits")

The conventional FITS WCS is defined by keywords in the FITS file and is automatically parsed by astropy.wcs.WCS

hdulist[0].header[:22]
SIMPLE  =                    T / Written by IDL:  Mon Dec 16 16:40:46 2019      
BITPIX  =                  -64 /  IEEE double precision floating point          
NAXIS   =                    2 /                                                
NAXIS1  =                 4096 /                                                
NAXIS2  =                 4096 /                                                
INSTR   = 'IRIS    '           /                                                
SCALE   =           0.00400000 / pixel scale (arcsec)                           
UNITS   = 'electrons'          /                                                
COORD_SY= 'C       '           /                                                
RADECSYS= 'FK5     '           /                                                
CTYPE1  = 'RA--LINEAR'         /                                                
CTYPE2  = 'DEC-LINEAR'         /                                                
CUNIT1  = 'deg     '           /                                                
CUNIT2  = 'deg     '           /                                                
CRPIX1  =              2048.12 /                                                
CRPIX2  =              2048.12 /                                                
CRVAL1  =        265.197723389 /                                                
CRVAL2  =       -28.9921894073 /                                                
CDELT1  =    1.11111013731E-06 /0.00400000 pixel scale                          
CDELT2  =    1.11111013731E-06 /0.00400000 pixel scale                          
CROTA1  =                    0 /                                                
CROTA2  =                    0 /                                                

Cannot make LINEAR to work, so let's instead replace it with Gnomonic

hdulist[0].header["CTYPE1"] = "RA---TAN"
hdulist[0].header["CTYPE2"] = "DEC--TAN"
w = wcs.WCS(hdulist[0])
WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / 
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
WARNING: FITSFixedWarning: EPOCH = '2019-12-17T00:40:46.48107737302794Z' / 
a floating-point value was expected. [astropy.wcs.wcs]
w.to_header()
WCSAXES =                    2 / Number of coordinate axes                      
CRPIX1  =              2048.12 / Pixel coordinate of reference point            
CRPIX2  =              2048.12 / Pixel coordinate of reference point            
CDELT1  =    1.11111013731E-06 / [deg] Coordinate increment at reference point  
CDELT2  =    1.11111013731E-06 / [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  =        265.197723389 / [deg] Coordinate value at reference point      
CRVAL2  =       -28.9921894073 / [deg] Coordinate value at reference point      
LONPOLE =                180.0 / [deg] Native longitude of celestial pole       
LATPOLE =       -28.9921894073 / [deg] Native latitude of celestial pole        
MJDREFI =                  0.0 / [d] MJD of fiducial time, integer part         
MJDREFF =                  0.0 / [d] MJD of fiducial time, fractional part      
RADESYS = 'FK5'                / Equatorial coordinate system                   
EQUINOX =               2000.0 / [yr] Equinox of equatorial coordinates         

We can then convert between the pixel indices and the coordinates in the sky

pixcrd = np.array([[0, 0], [0, 4096],[4096, 4096], [4096,0]], dtype=np.float64)
world = w.wcs_pix2world(pixcrd, 0)
print(world)
[[265.19512288 -28.99446396]
 [265.195123   -28.98991285]
 [265.20032602 -28.98991285]
 [265.20032613 -28.99446396]]

We can calculate the size of the instrument field of view

((-world[0][0]+ world[-1][0]) * u.deg).to(u.arcsec)
$18.731693 \; \mathrm{{}^{\prime\prime}}$
((-world[0][1]+ world[1][1]) * u.deg).to(u.arcsec)
$16.383986 \; \mathrm{{}^{\prime\prime}}$
w
WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN'  'DEC--TAN'  
CRVAL : 265.197723389  -28.9921894073  
CRPIX : 2048.12  2048.12  
NAXIS : 4096  4096

Create a gwcs WCS object

We want now to use astropy.modeling to build a transformation that is equivalent to the FITS WCS transformation defined above

from gwcs import WCS
from astropy.modeling import models
from astropy import coordinates as coord
from astropy import units as u
from gwcs import coordinate_frames as cf
shift_by_crpix = models.Shift(-(hdulist[0].header["CRPIX1"] - 1)*u.pix) & models.Shift(-(hdulist[0].header["CRPIX2"] - 1)*u.pix)
tan = models.Pix2Sky_TAN()
celestial_rotation =  models.RotateNative2Celestial(
    hdulist[0].header["CRVAL1"]*u.deg, hdulist[0].header["CRVAL2"]*u.deg, 180*u.deg)
tan.input_units_equivalencies = {"x": u.pixel_scale(hdulist[0].header["CDELT1"] *u.deg/u.pix),
                                      "y": u.pixel_scale(hdulist[0].header["CDELT2"] *u.deg/u.pix)}
det2sky = shift_by_crpix | tan | celestial_rotation
det2sky.name = "linear_transform"
detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"),
                            unit=(u.pix, u.pix))
sky_frame = cf.CelestialFrame(reference_frame=coord.FK5(), name='fk5',
                              unit=(u.deg, u.deg))
pipeline = [(detector_frame, det2sky),
            (sky_frame, None)
           ]
wcsobj = WCS(pipeline)
print(wcsobj)
  From      Transform    
-------- ----------------
detector linear_transform
     fk5             None
 wcsobj(pixcrd[:,0]*u.pix, pixcrd[:,1]*u.pix, with_units=False)
(<Quantity [265.19512288, 265.195123  , 265.20032602, 265.20032613] deg>,
 <Quantity [-28.99446396, -28.98991285, -28.98991285, -28.99446396] deg>)
((-_[0][0]+ _[0][-1])).to(u.arcsec)
$18.731693 \; \mathrm{{}^{\prime\prime}}$
((-__[1][0]+ __[1][1])).to(u.arcsec)
$16.383986 \; \mathrm{{}^{\prime\prime}}$
wcsobj
<WCS(output_frame=fk5, input_frame=detector, forward_transform=Model: CompoundModel
Name: linear_transform
Inputs: ('x0', 'x1')
Outputs: ('alpha_C', 'delta_C')
Model set size: 1
Expression: [0] & [1] | [2] | [3]
Components: 
    [0]: <Shift(offset=-2047.12 pix)>

    [1]: <Shift(offset=-2047.12 pix)>

    [2]: <Pix2Sky_Gnomonic()>

    [3]: <RotateNative2Celestial(lon=265.19772339 deg, lat=-28.99218941 deg, lon_pole=180. deg)>
Parameters:
    offset_0 offset_1     lon_3            lat_3        lon_pole_3
      pix      pix         deg              deg            deg    
    -------- -------- ------------- ------------------- ----------
    -2047.12 -2047.12 265.197723389 -28.992189407300003      180.0)>