Fix a bug in wcs_from_sregions when crpix is provided

changes/326.bugfix.rst

@@ -0,0 +1,4 @@
+Fix a bug in ``alignment.util.wcs_from_sregions()`` and
+``alignment.util.wcs_from_footprints()`` functions due to which, when the
+``crpix`` argument was provided, the bounding box of the created WCS was
+not adjusted to account for this ``crpix`` value.

src/stcal/alignment/util.py

@@ -227,13 +227,18 @@ def _calculate_offsets(fiducial: tuple,
         A two-elements array containing the minimum pixel value for each axis.
 
     crpix : list or tuple
-        Pixel coordinates of the reference pixel.
+        0-indexed pixel coordinates of the reference pixel.
 
     Returns
     -------
     ~astropy.modeling.Model
         A model with the offsets to be added to the WCS's transform.
 
+    tuple
+        A tuple of offset *values* that are to be *subtracted* from input
+        coordinates (negative values of offsets in the returned transform).
+        Note: these are equivalent to an effective 0-indexed "crpix".
+
     Notes
     -----
     If ``crpix=None``, then ``fiducial``, ``wcs``, and ``axis_min_values`` must be
@@ -249,9 +254,10 @@ def _calculate_offsets(fiducial: tuple,
         msg = "If crpix is not provided, fiducial, wcs, and axis_min_values must be provided."
         raise ValueError(msg)
     else:
+        # assume 0-based CRPIX
         offset1, offset2 = crpix
 
-    return astmodels.Shift(-offset1, name="crpix1") & astmodels.Shift(-offset2, name="crpix2")
+    return astmodels.Shift(-offset1, name="crpix1") & astmodels.Shift(-offset2, name="crpix2"), (offset1, offset2)
 
 
 def _calculate_new_wcs(wcs: gwcs.wcs.WCS,
@@ -283,7 +289,7 @@ def _calculate_new_wcs(wcs: gwcs.wcs.WCS,
         coordinate system.
 
     crpix : tuple, optional
-        The coordinates of the reference pixel.
+        0-indexed coordinates of the reference pixel.
 
     transform : ~astropy.modeling.Model
         An optional transform to be prepended to the transform constructed by the
@@ -302,14 +308,23 @@ def _calculate_new_wcs(wcs: gwcs.wcs.WCS,
         transform=transform,
         input_frame=wcs.input_frame,
     )
-    axis_min_values, output_bounding_box = _get_axis_min_and_bounding_box(footprints, wcs_new)
-    offsets = _calculate_offsets(
+    axis_min_values, bbox = _get_axis_min_and_bounding_box(footprints, wcs_new)
+    offsets, shifts = _calculate_offsets(
         fiducial=fiducial,
         wcs=wcs_new,
         axis_min_values=axis_min_values,
         crpix=crpix,
     )
 
+    if crpix is None:
+        output_bounding_box = bbox
+    else:
+        output_bounding_box = []
+        for axis_range, shift in zip(bbox, shifts):
+            output_bounding_box.append(
+                (axis_range[0] + shift, axis_range[1] + shift)
+            )
+
     wcs_new.insert_transform("detector", offsets, after=True)
     wcs_new.bounding_box = output_bounding_box
 

tests/test_alignment.py

@@ -30,7 +30,7 @@ def _create_wcs_object_without_distortion(
     shape,
 ):
     # subtract 1 to account for pixel indexing starting at 0
-    shift = models.Shift() & models.Shift()
+    shift = models.Shift(0) & models.Shift(0)
 
     scale = models.Scale(pscale[0]) & models.Scale(pscale[1])
 
@@ -171,7 +171,7 @@ def test_sregion_to_footprint():
 
     footprint = _sregion_to_footprint(s_region)
 
-    assert footprint.shape == (4,2)
+    assert footprint.shape == (4, 2)
     assert np.allclose(footprint, expected_footprint)
 
 
@@ -199,18 +199,27 @@ def test_wcs_from_footprints(s_regions):
     )
     dm_2 = _create_wcs_and_datamodel(fiducial_world, shape, pscale)
     wcs_2 = dm_2.meta.wcs
+
     if s_regions:
         footprints = [wcs_1.footprint(), wcs_2.footprint()]
         wcs = wcs_from_sregions(footprints, wcs_1, dm_1.meta.wcsinfo.instance)
     else:
         wcs_list = [wcs_1, wcs_2]
-        wcs = wcs_from_footprints(wcs_list, wcs_1, dm_1.meta.wcsinfo.instance)
-
-    # check that all elements of footprint match the *vertices* of the new combined WCS
-    assert all(np.isclose(wcs.footprint()[0], wcs(0, 0)))
-    assert all(np.isclose(wcs.footprint()[1], wcs(0, 4)))
-    assert all(np.isclose(wcs.footprint()[2], wcs(4, 4)))
-    assert all(np.isclose(wcs.footprint()[3], wcs(4, 0)))
+        msg = "wcs_from_footprints is deprecated and will be removed"
+        with pytest.warns(DeprecationWarning, match=msg):
+            wcs = wcs_from_footprints(
+                wcs_list,
+                wcs_1,
+                dm_1.meta.wcsinfo.instance
+            )
+
+    # check that all elements of footprint match the *vertices* of the new
+    # combined WCS
+    footprnt = wcs.footprint()
+    assert all(np.isclose(footprnt[0], wcs(0, 0)))
+    assert all(np.isclose(footprnt[1], wcs(0, 4)))
+    assert all(np.isclose(footprnt[2], wcs(4, 4)))
+    assert all(np.isclose(footprnt[3], wcs(4, 0)))
 
     # check that fiducials match their expected coords in the new combined WCS
     assert all(np.isclose(wcs_1(0, 0), wcs(2.5, 1.5)))