Imaging WCS

image_wcs.py

import numpy as np
from astropy import wcs
from gwcs import wcs as gw
from astropy.modeling.models import Shift, Scale, AffineTransformation2D, RotateNative2Celestial, Pix2Sky_TAN

crpix1 = 512.5
crpix2 = 512.5
cdelt1 = 2.402792
cdelt2 = 2.402792
ctype1 = 'RA---TAN'
ctype2 = 'DEC--TAN'
pc = np.array([[0.99999706, 0.00242302],[-0.00242302, 0.999997056]])
crval1 = 0.00089530541880
crval2 = 0.000384939
w = wcs.WCS()
w.wcs.crpix = [crpix1, crpix2]
w.wcs.cdelt = [cdelt1 / 3600, cdelt2 / 3600]
w.wcs.pc = pc / 3600
w.wcs.ctype = [ctype1, ctype2]
w.wcs.crval = [crval1 / 3600, crval2 /3600]

In [23]: w.printwcs()
WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN'  'DEC--TAN'  
CRVAL : 0.00089530541879999999  0.000384939  
CRPIX : 512.5  512.5  
PC1_1 PC1_2  : 0.99999705999999999  0.0024230200000000001  
PC2_1 PC2_2  : -0.0024230200000000001  0.999997056  
CDELT : 2.4027919999999998  2.4027919999999998  
NAXIS : 0  0

# Looks good - goes to w.wcs.crval
w.all_pix2world(crpix1,crpix2,1)
#Out[24]: array(2.4869594966666667e-07), array(1.0692750151974906e-07)]

# Now try astropy.wcs.with values in arcsec
warcsec = wcs.WCS()
warcsec.wcs.crpix = [crpix1, crpix2]
warcsec.wcs.cdelt = [cdelt1, cdelt2]
warcsec.wcs.pc = pc
warcsec.wcs.ctype = [ctype1, ctype2]
warcsec.wcs.crval = [crval1, crval2]
warcsec.wcs.cunit = ['arcsec', 'arcsec']

# It looks good with astropy.wcs expressed in units of arcsec
warcsec.all_pix2world(crpix1, crpix2, 1)
Out[97]: [array(2.4869594966666667e-07), array(1.0692750151974906e-07)]


# Now do it with modeling without units

shift = Shift(-crpix1) & Shift(-crpix2)
rot = AffineTransformation2D(pc / 3600)
scl = Scale(cdelt1/ 3600) & Scale(cdelt2 / 3600)
tan = Pix2Sky_TAN()
skyrot = RotateNative2Celestial(crval1/3600, crval2/3600, 180)
trans = shift | rot | scl | tan | skyrot

# Looks good wihtout units
trans(crpix1, crpix2)
Out[99]: (2.486959509774181e-07, 1.0692750039155865e-07)


# with units - an error is raised because compound models don't work with units
from astropy import units as u
shiftu = Shift(-crpix1 * u.pix) & Shift(-crpix2 * u.pix)
cdelt = np.array([[cdelt1, 1], [1, cdelt2]])
matrixu = cdelt * pc  
rotu = AffineTransformation2D(matrix=matrixu * u.arcsec, translation=[0,0]*u.arcsec )
rot.input_units_equivalencies={'x':u.pixel_scale(1*u.arcsec/u.pix), 'y': u.pixel_scale(1*u.arcsec/u.pix)}
tanu = Pix2Sky_TAN()
skyrotu = RotateNative2Celestial(crval1*u.arcsec, crval2*u.arcsec, 180*u.deg)
transform_with_units = shiftu | rotu | tanu | skyrotu

transform_with_units(crpix1*u.pix, crpix2*u.pix)
Out[28]: (2.486959509774181e-07, 1.0692750039155865e-07)