JP-1585: Wrap first angle at 360 in forward V2V3 to sky and at 180 fo…

…r inverse. (#5206)

CHANGES.rst

@@ -259,6 +259,14 @@ stpipe
 
 - Enable prefetch of pars reference files for associations. [#5249]
 
+transforms
+----------
+
+- Wrap first spherical angle ("RA") at 360 degrees in the forward ``V23ToSky``
+  transformation and to 180 degrees for the inverse transformation ("V2").
+  This is now done using models defined in ``astropy`` and ``gwcs`` packages
+  replacing ``V23ToSky`` model in JWST's WCS pipeline. [#5206]
+
 wavecorr
 --------
 

jwst/assign_wcs/pointing.py

@@ -2,29 +2,36 @@
 from astropy import units as u
 from astropy import coordinates as coords
 from astropy.modeling import models as astmodels
-from ..datamodels import DataModel
+from astropy.modeling.models import Scale, RotationSequence3D
+from astropy.coordinates.matrix_utilities import rotation_matrix, matrix_product
 from gwcs import utils as gwutils
+from gwcs.geometry import SphericalToCartesian, CartesianToSpherical
 from gwcs import coordinate_frames as cf
 from gwcs import wcs
-from ..transforms.models import V23ToSky
+
+from ..datamodels import DataModel
 
 
 __all__ = ["compute_roll_ref", "frame_from_model", "fitswcs_transform_from_model"]
 
 
-def v23tosky(input_model):
+def v23tosky(input_model, wrap_v2_at=180, wrap_lon_at=360):
     v2_ref = input_model.meta.wcsinfo.v2_ref / 3600
     v3_ref = input_model.meta.wcsinfo.v3_ref / 3600
     roll_ref = input_model.meta.wcsinfo.roll_ref
     ra_ref = input_model.meta.wcsinfo.ra_ref
     dec_ref = input_model.meta.wcsinfo.dec_ref
 
-    angles = [-v2_ref, v3_ref, -roll_ref, -dec_ref, ra_ref]
+    angles = np.array([v2_ref, -v3_ref, roll_ref, dec_ref, -ra_ref])
     axes = "zyxyz"
-    sky_rotation = V23ToSky(angles, axes_order=axes, name="v23tosky")
+    rot = RotationSequence3D(angles, axes_order=axes)
+
     # The sky rotation expects values in deg.
     # This should be removed when models work with quantities.
-    return astmodels.Scale(1/3600) & astmodels.Scale(1/3600) | sky_rotation
+    m = ((Scale(1/3600) & Scale(1/3600)) | SphericalToCartesian(wrap_lon_at=wrap_v2_at)
+         | rot | CartesianToSpherical(wrap_lon_at=wrap_lon_at))
+    m.name = 'v23tosky'
+    return m
 
 
 def compute_roll_ref(v2_ref, v3_ref, roll_ref, ra_ref, dec_ref, new_v2_ref, new_v3_ref):
@@ -57,13 +64,15 @@ def compute_roll_ref(v2_ref, v3_ref, roll_ref, ra_ref, dec_ref, new_v2_ref, new_
     v3_ref = v3_ref / 3600
 
     if np.isclose(v2_ref, 0, atol=1e-13) and np.isclose(v3_ref, 0, atol=1e-13):
-        angles = [-roll_ref, -dec_ref, - ra_ref]
+        angles = [roll_ref, dec_ref, ra_ref]
         axes = "xyz"
     else:
-        angles = [-v2_ref, v3_ref, -roll_ref, -dec_ref, ra_ref]
+        angles = [v2_ref, -v3_ref, roll_ref, dec_ref, -ra_ref]
         axes = "zyxyz"
-    M = V23ToSky._compute_matrix(np.deg2rad(angles), axes)
-    return _roll_angle_from_matrix(M, v2, v3)
+
+    matrices = [rotation_matrix(a, ax) for a, ax in zip(angles, axes)]
+    m = matrix_product(*matrices[::-1])
+    return _roll_angle_from_matrix(m, v2, v3)
 
 
 def _roll_angle_from_matrix(matrix, v2, v3):

jwst/lib/set_telescope_pointing.py

@@ -10,6 +10,7 @@
 import numpy as np
 
 from ..assign_wcs.util import update_s_region_keyword, calc_rotation_matrix
+from ..assign_wcs.pointing import v23tosky
 from ..datamodels import Level1bModel
 from ..lib.engdb_tools import ENGDB_Service
 from .exposure_types import IMAGING_TYPES, FGS_GUIDE_EXP_TYPES
@@ -530,12 +531,7 @@ def update_s_region(model, siaf):
         "Vertices for aperture {0}: {1}".format(model.meta.aperture.name, vertices)
     )
     # Execute IdealToV2V3, followed by V23ToSky
-    from ..transforms.models import IdealToV2V3, V23ToSky
-    v2_ref_deg = model.meta.wcsinfo.v2_ref / 3600  # in deg
-    v3_ref_deg = model.meta.wcsinfo.v3_ref / 3600  # in deg
-    roll_ref = model.meta.wcsinfo.roll_ref
-    ra_ref = model.meta.wcsinfo.ra_ref
-    dec_ref = model.meta.wcsinfo.dec_ref
+    from ..transforms.models import IdealToV2V3
     vparity = model.meta.wcsinfo.vparity
     v3yangle = model.meta.wcsinfo.v3yangle
 
@@ -547,15 +543,10 @@ def update_s_region(model, siaf):
     )
     v2, v3 = idltov23(xvert, yvert)  # in arcsec
 
-    # Convert to deg
-    v2 = v2 / 3600  # in deg
-    v3 = v3 / 3600  # in deg
-    angles = [-v2_ref_deg, v3_ref_deg, -roll_ref, -dec_ref, ra_ref]
-    axes = "zyxyz"
-    v23tosky = V23ToSky(angles, axes_order=axes)
-    ra_vert, dec_vert = v23tosky(v2, v3)
-    negative_ind = ra_vert < 0
-    ra_vert[negative_ind] = ra_vert[negative_ind] + 360
+    # hardcode wrapping angles for V2 and RA here. Could be made more
+    # flexible if needed.
+    v23tosky_tr = v23tosky(model, wrap_v2_at=180, wrap_lon_at=360)
+    ra_vert, dec_vert = v23tosky_tr(v2, v3)
     # Do not do any sorting, use the vertices in the SIAF order.
     footprint = np.array([ra_vert, dec_vert]).T
     update_s_region_keyword(model, footprint)