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)
((-world[0][1]+ world[1][1]) * u.deg).to(u.arcsec)
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
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)
((-__[1][0]+ __[1][1])).to(u.arcsec)
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)>