JP-3824: Fix a bug in wcs_from_sregions when crpix is provided

changes/326.bugfix.rst

@@ -0,0 +1,3 @@
+Fix a bug in ``alignment.util.wcs_from_sregions()`` and
+``alignment.util.wcs_from_footprints()`` functions that was causing incorrect WCS bounding boxes
+when the ``crpix`` argument was provided.

src/stcal/alignment/util.py

@@ -4,7 +4,8 @@
 import functools
 import logging
 import re
-from typing import TYPE_CHECKING
+import typing
+from typing import TYPE_CHECKING, Tuple
 import warnings
 
 if TYPE_CHECKING:
@@ -143,10 +144,16 @@
     return transform
 
 
-def _get_axis_min_and_bounding_box(footprints: list[np.ndarray],
-                                   ref_wcs: gwcs.wcs.WCS) -> tuple:
+def _get_bounding_box_with_offsets(
+        footprints: list[np.ndarray],
+        ref_wcs: gwcs.wcs.WCS,
+        fiducial: Sequence,
+        crpix: Sequence | None,
+        shape: Sequence | None
+        ) -> Tuple[tuple, tuple, astmodels.Model]:
     """
     Calculates axis minimum values and bounding box.
+    Calculates the offsets to the transform.
 
     Parameters
     ----------
@@ -157,28 +164,74 @@
     ref_wcs : ~gwcs.wcs.WCS
         The reference WCS object.
 
+    fiducial : tuple
+        A tuple containing the world coordinates of the fiducial point.
+
+    crpix : list or tuple, optional
+        0-indexed pixel coordinates of the reference pixel.
+
+    shape : tuple, optional
+        Shape (using `numpy.ndarray` convention) of the image array associated
+        with the ``ref_wcs``.
+
+
     Returns
     -------
+
+    tuple
+        A tuple containing the bounding box region in the format
+        ((x0_lower, x0_upper), (x1_lower, x1_upper), ...).
+
     tuple
-        A tuple containing two elements:
-            1 - a :py:class:`np.ndarray` with the minimum value in each axis;
-            2 - a tuple containing the bounding box region in the format
-            ((x0_lower, x0_upper), (x1_lower, x1_upper)).
+        Shape of the image. When ``shape`` argument is `None`, shape is
+        determined from the upper limit of the computed bounding box, otherwise
+        input value of ``shape`` is returned.
+
+    ~astropy.modeling.Model
+        A model with the offsets to be added to the WCS's transform.
+
     """
     domain_bounds = np.hstack([ref_wcs.backward_transform(*f.T) for f in footprints])
-    axis_min_values = np.min(domain_bounds, axis=1)
-    domain_bounds = (domain_bounds.T - axis_min_values).T
-
-    output_bounding_box = []
-    for axis in ref_wcs.output_frame.axes_order:
-        axis_min, axis_max = (
-            domain_bounds[axis].min(),
-            domain_bounds[axis].max(),
+    domain_min = np.min(domain_bounds, axis=1)
+    domain_max = np.max(domain_bounds, axis=1)
+
+    native_crpix = ref_wcs.backward_transform(*fiducial)
+
+    if crpix is None:
+        # shift the coordinates by domain_min so that all input footprints
+        # will project to positive coordinates in the offsetted ref_wcs
+        offsets = tuple(ncrp - dmin for ncrp, dmin in zip(native_crpix, domain_min))
+
+    else:
+        # assume 0-based CRPIX and require that fiducial would map to the user
+        # defined crpix value:
+        offsets = tuple(
+            crp - ncrp for ncrp, crp in zip(native_crpix, crpix)
         )
-        # populate output_bounding_box
-        output_bounding_box.append((axis_min, axis_max))
 
-    return (axis_min_values, output_bounding_box)
+    # Also offset domain limits:
+    domain_min += offsets
+    domain_max += offsets
+
+    if shape is None:
+        shape = tuple(int(dmax + 0.5) for dmax in domain_max[::-1])
+        bounding_box = tuple(
+            (-0.5, s - 0.5) for s in shape[::-1]
+        )
+
+    else:
+        # trim upper bounding box limits
+        shape = tuple(shape)
+        bounding_box = tuple(
+            (max(0, int(dmin + 0.5)) - 0.5, min(int(dmax + 0.5), sz) - 0.5)
+            for dmin, dmax, sz in zip(domain_min, domain_max, shape[::-1])
+        )
+
+    model = astmodels.Shift(-offsets[0], name="crpix1")
+    for k, shift in enumerate(offsets[1:]):
+        model = model.__and__(astmodels.Shift(-shift, name=f"crpix{k + 2:d}"))
+
+    return bounding_box, shape, model
 
 
 def _calculate_fiducial(footprints: list[np.ndarray],
@@ -208,52 +261,6 @@
     return _compute_fiducial_from_footprints(footprints)
 
 
-def _calculate_offsets(fiducial: tuple,
-                       wcs: gwcs.wcs.WCS | None,
-                       axis_min_values: np.ndarray | None,
-                       crpix: Sequence | None) -> astmodels.Model:
-    """
-    Calculates the offsets to the transform.
-
-    Parameters
-    ----------
-    fiducial : tuple
-        A tuple containing the world coordinates of the fiducial point.
-
-    wcs : ~gwcs.wcs.WCS
-        A WCS object. It will be used to determine the
-
-    axis_min_values : np.ndarray
-        A two-elements array containing the minimum pixel value for each axis.
-
-    crpix : list or tuple
-        Pixel coordinates of the reference pixel.
-
-    Returns
-    -------
-    ~astropy.modeling.Model
-        A model with the offsets to be added to the WCS's transform.
-
-    Notes
-    -----
-    If ``crpix=None``, then ``fiducial``, ``wcs``, and ``axis_min_values`` must be
-    provided, in which case, the offsets will be calculated using the WCS object to
-    find the pixel coordinates of the fiducial point and then correct it by the minimum
-    pixel value for each axis.
-    """
-    if crpix is None and fiducial is not None and wcs is not None and axis_min_values is not None:
-        offset1, offset2 = wcs.backward_transform(*fiducial)
-        offset1 -= axis_min_values[0]
-        offset2 -= axis_min_values[1]
-    elif crpix is None:
-        msg = "If crpix is not provided, fiducial, wcs, and axis_min_values must be provided."
-        raise ValueError(msg)
-    else:
-        offset1, offset2 = crpix
-
-    return astmodels.Shift(-offset1, name="crpix1") & astmodels.Shift(-offset2, name="crpix2")
-
-
 def _calculate_new_wcs(wcs: gwcs.wcs.WCS,
                        shape: Sequence | None,
                        footprints: list[np.ndarray],
@@ -283,7 +290,7 @@
         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,20 +309,23 @@
         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(
+
+    bounding_box, shape, offsets = _get_bounding_box_with_offsets(
+        footprints,
+        ref_wcs=wcs_new,
         fiducial=fiducial,
-        wcs=wcs_new,
-        axis_min_values=axis_min_values,
         crpix=crpix,
+        shape=shape
     )
 
-    wcs_new.insert_transform("detector", offsets, after=True)
-    wcs_new.bounding_box = output_bounding_box
-
-    if shape is None:
-        shape = [int(axs[1] - axs[0] + 0.5) for axs in output_bounding_box[::-1]]
+    if any(d < 2 for d in shape):
+        raise ValueError(
+            "Computed shape for the output image using provided "
+            "WCS parameters is too small."
+        )
 
+    wcs_new.insert_transform("detector", offsets, after=True)
+    wcs_new.bounding_box = bounding_box
     wcs_new.pixel_shape = shape[::-1]
     wcs_new.array_shape = shape
     return wcs_new
@@ -826,9 +836,7 @@
     footprint = footprint[:2, :]
 
     # Make sure RA values are all positive
-    negative_ind = footprint[0] < 0
-    if negative_ind.any():
-        footprint[0][negative_ind] = 360 + footprint[0][negative_ind]
+    np.mod(footprint[0], 360, out=footprint[0])
 
     footprint = footprint.T
     return compute_s_region_keyword(footprint)

tests/test_alignment.py

@@ -29,9 +29,8 @@ def _create_wcs_object_without_distortion(
     pscale,
     shape,
 ):
-    # subtract 1 to account for pixel indexing starting at 0
+    # subtract 0 shift to mimic a resampled WCS that does include shift transforms
     shift = models.Shift() & models.Shift()
-
     scale = models.Scale(pscale[0]) & models.Scale(pscale[1])
 
     tan = models.Pix2Sky_TAN()
@@ -171,7 +170,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 +198,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.5, -0.5)))
+    assert all(np.isclose(footprnt[1], wcs(-0.5, 3.5)))
+    assert all(np.isclose(footprnt[2], wcs(3.5, 3.5)))
+    assert all(np.isclose(footprnt[3], wcs(3.5, -0.5)))
 
     # 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)))

tests/test_ramp_fitting_gls_fit.py

@@ -678,12 +678,11 @@ def test_two_groups_unc():
     answer = slopes[2][50, 50]
     check = np.sqrt((deltaDN / gain) / group_time**2 + (rnoise**2 / group_time**2))
     tol = 1.0e-6
-    # print(f"answer = {answer}")
-    # print(f"check = {check}")
+
     np.testing.assert_allclose(answer, check, tol)
 
 
-@pytest.mark.skip(reason="GLS does not comopute VAR_XXX arrays.")
+@pytest.mark.skip(reason="GLS does not compute VAR_XXX arrays.")
 def test_five_groups_unc():
     """ """
     """