import astropy
from astropy.io import fits
import numpy as np
from astropy import units as uSome notes on how to convert a FITS WCS to a WCS from gwcs (Generalized WCS) to be used within the JWST pipeline
from astropy import wcshdulist = 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}}\)
wWCS Keywords
Number of WCS axes: 2
CTYPE : 'RA---TAN' 'DEC--TAN'
CRVAL : 265.197723389 -28.9921894073
CRPIX : 2048.12 2048.12
NAXIS : 4096 4096Create 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 WCSfrom astropy.modeling import models
from astropy import coordinates as coord
from astropy import units as u
from gwcs import coordinate_frames as cfshift_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)>