>>> from astropy.modeling.models import (Shift, Scale, Rotation2D, Pix2Sky_TAN, RotateNative2Celestial, Mapping, Polynomial2D, AffineTransformation2D)
>>> import gwcs
>>> from gwcs import coordinate_frames as cf
>>> import numpy as np
>>> from astropy import coordinates as coord
>>>
>>> CRPIX1 = 2048.0
>>> CRPIX2 = 1024.0
>>> CRVAL1 = 5.63056810618
>>> CRVAL2 = -72.0545718428
>>> LONPOLE = 180
>>> PC1_1 = 1.29058668e-05
>>> PC1_2 = 5.95320246e-06
>>> PC2_1 = 5.02215196e-06
>>> PC2_2 = -1.26450104e-05
>>> CTYPE1 = 'RA---TAN'
>>> CTYPE2 = 'DEC---TAN'
>>> dist_x = Polynomial2D(2, c0_0=0.0013, c1_0=0.5, c2_0=0, c0_1=0.823, c0_2=1.4, c1_1=1.7, name='x_distortion')
>>> dist_y = Polynomial2D(2, c0_0=0.03, c1_0=0.25, c2_0=1.2, c0_1=0.3, c0_2=0.4, c1_1=0.7, name='y_distortion')
>>> distortion = dist_x & dist_y
>>> distortion_mapping = Mapping((0, 1, 0, 1), name='distortion_mapping')
>>> distortion = distortion_mapping | dist_x & dist_y
>>> shift = Shift(CRPIX1, name="crpix1") & Shift(CRPIX2, name="crpix2")
>>> plane_rotation = AffineTransformation2D(matrix=np.array([[PC1_1, PC1_2], [PC2_1, PC2_2]]))
>>> tangent = Pix2Sky_TAN()
>>> sky_rotation = RotateNative2Celestial(CRVAL1, CRVAL2, LONPOLE)
>>> focal2sky = shift | plane_rotation | tangent | sky_rotation
>>>
>>> inframe = cf.Frame2D(name='detector')
>>> focal = cf.Frame2D(name='focal')
>>> outframe = cf.CelestialFrame(reference_frame=coord.ICRS(), name='CelestialFrame')
>>> pipeline = [(inframe, distortion), (focal, focal2sky), (outframe, None)]
>>> memw = gwcs.wcs.WCS(pipeline)
>>>
>>> memw(1,1)
(5.736394270266226, -72.0572205208842)
>>> tr = memw.get_transform('focal', 'CelestialFrame')
>>> tr
<CompoundModel6(offset_0=2048.0, offset_1=1024.0, matrix_2=[[ 0.00001291, 0.00000595], [ 0.00000502,-0.00001265]], translation_2=[ 0., 0.], lon_4=5.63056810618, lat_4=-72.0545718428, lon_pole_4=180.0)>
>>> tr(1,1)
(5.7362144489763285, -72.05721404287078)
>>> tr['crpix1'].offset=0
>>> tr
<CompoundModel6(offset_0=0.0, offset_1=1024.0, matrix_2=[[ 0.00001291, 0.00000595], [ 0.00000502,-0.00001265]], translation_2=[ 0., 0.], lon_4=5.63056810618, lat_4=-72.0545718428, lon_pole_4=180.0)>
>>> tr(1,1)
(5.650428514863572, -72.0675269477636)
>>>
>>> # memw before replacing focal->celestial transformation:
>>> memw
<WCS(output_frame=CelestialFrame, input_frame=detector, forward_transform=Model: CompoundModel11
Inputs: ('x0', 'x1')
Outputs: (u'alpha_C', u'delta_C')
Model set size: 1
Expression: [0] | [1] & [2] | [3] & [4] | [5] | [6] | [7]
Components:
[0]: <Mapping((0, 1, 0, 1), name=distortion_mapping)>
[1]: <Polynomial2D(2, c0_0=0.0013, c1_0=0.5, c2_0=0.0, c0_1=0.823, c0_2=1.4, c1_1=1.7, name='x_distortion')>
[2]: <Polynomial2D(2, c0_0=0.03, c1_0=0.25, c2_0=1.2, c0_1=0.3, c0_2=0.4, c1_1=0.7, name='y_distortion')>
[3]: <Shift(offset=2048.0, name='crpix1')>
[4]: <Shift(offset=1024.0, name='crpix2')>
[5]: <AffineTransformation2D(matrix=[[ 0.00001291, 0.00000595], [ 0.00000502,-0.00001265]], translation=[ 0., 0.])>
[6]: <Pix2Sky_Gnomonic()>
[7]: <RotateNative2Celestial(lon=5.63056810618, lat=-72.0545718428, lon_pole=180.0)>
Parameters:
c0_0_1 c1_0_1 c2_0_1 c0_1_1 c0_2_1 c1_1_1 c0_0_2 c1_0_2 c2_0_2 c0_1_2 c0_2_2 c1_1_2 offset_3 offset_4 matrix_5 [2,2] translation_5 [2] lon_7 lat_7 lon_pole_7
------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ -------- -------- --------------------------------- ----------------- ------------- -------------- ----------
0.0013 0.5 0.0 0.823 1.4 1.7 0.03 0.25 1.2 0.3 0.4 0.7 2048.0 1024.0 1.29058668e-05 .. -1.26450104e-05 0.0 .. 0.0 5.63056810618 -72.0545718428 180.0)>
>>> memw(1,1)
(5.736394270266226, -72.0572205208842)
>>>
>>> # Now replace focal->celestial using the modified 'tr' which we have checked above that it works:
>>> memw.set_transform('focal', 'CelestialFrame', tr)
>>> memw(1,1)
(5.736394270266226, -72.0572205208842)
>>> # Above result is identical to the result obtained before WCS was modified.
>>>
>>> # Let's check that focal->celestial transformation in the WCS is the one that we just set:
>>> memw.get_transform('focal', 'CelestialFrame')
<CompoundModel6(offset_0=0.0, offset_1=1024.0, matrix_2=[[ 0.00001291, 0.00000595], [ 0.00000502,-0.00001265]], translation_2=[ 0., 0.], lon_4=5.63056810618, lat_4=-72.0545718428, lon_pole_4=180.0)>
>>> # Indeed, this is identical to 'tr':
>>> tr
<CompoundModel6(offset_0=0.0, offset_1=1024.0, matrix_2=[[ 0.00001291, 0.00000595], [ 0.00000502,-0.00001265]], translation_2=[ 0., 0.], lon_4=5.63056810618, lat_4=-72.0545718428, lon_pole_4=180.0)>
>>>
>>> # However, wcs object itself (full forward transformation) does not show modified 'CRPIX' value:
>>> memw
<WCS(output_frame=CelestialFrame, input_frame=detector, forward_transform=Model: CompoundModel14
Inputs: ('x0', 'x1')
Outputs: (u'alpha_C', u'delta_C')
Model set size: 1
Expression: [0] | [1] & [2] | [3] & [4] | [5] | [6] | [7]
Components:
[0]: <Mapping((0, 1, 0, 1), name=distortion_mapping)>
[1]: <Polynomial2D(2, c0_0=0.0013, c1_0=0.5, c2_0=0.0, c0_1=0.823, c0_2=1.4, c1_1=1.7, name='x_distortion')>
[2]: <Polynomial2D(2, c0_0=0.03, c1_0=0.25, c2_0=1.2, c0_1=0.3, c0_2=0.4, c1_1=0.7, name='y_distortion')>
[3]: <Shift(offset=2048.0, name='crpix1')>
[4]: <Shift(offset=1024.0, name='crpix2')>
[5]: <AffineTransformation2D(matrix=[[ 0.00001291, 0.00000595], [ 0.00000502,-0.00001265]], translation=[ 0., 0.])>
[6]: <Pix2Sky_Gnomonic()>
[7]: <RotateNative2Celestial(lon=5.63056810618, lat=-72.0545718428, lon_pole=180.0)>
Parameters:
c0_0_1 c1_0_1 c2_0_1 c0_1_1 c0_2_1 c1_1_1 c0_0_2 c1_0_2 c2_0_2 c0_1_2 c0_2_2 c1_1_2 offset_3 offset_4 matrix_5 [2,2] translation_5 [2] lon_7 lat_7 lon_pole_7
------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ ------ -------- -------- --------------------------------- ----------------- ------------- -------------- ----------
0.0013 0.5 0.0 0.823 1.4 1.7 0.03 0.25 1.2 0.3 0.4 0.7 2048.0 1024.0 1.29058668e-05 .. -1.26450104e-05 0.0 .. 0.0 5.63056810618 -72.0545718428 180.0)>