Disable bounding box in calc_pixmap (#164)

Allow disabling the bounding box on one or both input WCS

Author: mcara

.github/workflows/ci.yml

@@ -22,15 +22,13 @@ jobs:
       envs: |
         - linux: check-style
         - linux: check-security
-
         - linux: py310-xdist
         - linux: py311-xdist
         - linux: py313-xdist
         - macos: py312-xdist
         - windows: py312-xdist
         - linux: py312-xdist-cov
           coverage: codecov
-
         - linux: py312-xdist-devdeps
     secrets:
       CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}

CHANGES.rst

@@ -4,6 +4,15 @@
 Release Notes
 =============
 
+2.0.1 (Unreleased)
+==================
+
+- Update ``utils.calc_pixmap`` code to be ready for upcoming changes in GWCS
+  due to which inverse WCS transformations will respect bounding box by
+  allowing the caller of ``utils.calc_pixmap`` to disable the bounding box(es)
+  on both or one of the input WCS when computing the pixel map. [#164]
+
+
 2.0.0 (2024-10-23)
 ==================
 

drizzle/tests/test_utils.py

@@ -102,6 +102,63 @@ def test_translated_map(wcs_type):
     assert_almost_equal(pixmap[2:, 2:], ok_pixmap[2:, 2:], decimal=5)
 
 
+def test_disable_gwcs_bbox():
+    """
+    Map a pixel array to a translated version ofitself.
+    """
+    first_wcs = wcs_from_file(
+        "j8bt06nyq_sip_flt.fits",
+        ext=1,
+        wcs_type="gwcs"
+    )
+    second_wcs = wcs_from_file(
+        "j8bt06nyq_sip_flt.fits",
+        ext=1,
+        crpix_shift=(-2, -2),  # shift loaded WCS by adding this to CRPIX
+        wcs_type="gwcs"
+    )
+
+    ok_pixmap = np.indices(first_wcs.array_shape, dtype='float64') - 2.0
+    ok_pixmap = ok_pixmap.transpose()
+
+    # Mapping an array to a translated array
+
+    # disable both bounding boxes:
+    pixmap = calc_pixmap(first_wcs, second_wcs, disable_bbox="both")
+    assert_almost_equal(pixmap[2:, 2:], ok_pixmap[2:, 2:], decimal=5)
+    assert np.all(np.isfinite(pixmap[:2, :2]))
+    assert np.all(np.isfinite(pixmap[-2:, -2:]))
+    # check bbox was restored
+    assert first_wcs.bounding_box is not None
+    assert second_wcs.bounding_box is not None
+
+    # disable "from" bounding box:
+    pixmap = calc_pixmap(second_wcs, first_wcs, disable_bbox="from")
+    assert_almost_equal(pixmap[:-2, :-2], ok_pixmap[:-2, :-2] + 4.0, decimal=5)
+    assert np.all(np.logical_not(np.isfinite(pixmap[-2:, -2:])))
+    # check bbox was restored
+    assert first_wcs.bounding_box is not None
+    assert second_wcs.bounding_box is not None
+
+    # disable "to" bounding boxes:
+    pixmap = calc_pixmap(first_wcs, second_wcs, disable_bbox="to")
+    assert_almost_equal(pixmap[2:, 2:], ok_pixmap[2:, 2:], decimal=5)
+    assert np.all(np.isfinite(pixmap[:2, :2]))
+    assert np.all(pixmap[:2, :2] < 0.0)
+    assert np.all(np.isfinite(pixmap[-2:, -2:]))
+    # check bbox was restored
+    assert first_wcs.bounding_box is not None
+    assert second_wcs.bounding_box is not None
+
+    # enable all bounding boxes:
+    pixmap = calc_pixmap(first_wcs, second_wcs, disable_bbox="none")
+    assert_almost_equal(pixmap[2:, 2:], ok_pixmap[2:, 2:], decimal=5)
+    assert np.all(np.logical_not(np.isfinite(pixmap[:2, :2])))
+    # check bbox was restored
+    assert first_wcs.bounding_box is not None
+    assert second_wcs.bounding_box is not None
+
+
 def test_estimate_pixel_scale_ratio():
     w = wcs_from_file("j8bt06nyq_flt.fits", ext=1)
     pscale = estimate_pixel_scale_ratio(w, w, w.wcs.crpix, (0, 0))

drizzle/utils.py

@@ -7,7 +7,7 @@
 _DEG2RAD = math.pi / 180.0
 
 
-def calc_pixmap(wcs_from, wcs_to, shape=None):
+def calc_pixmap(wcs_from, wcs_to, shape=None, disable_bbox="to"):
     """
     Calculate a discretized on a grid mapping between the pixels of two images
     using provided WCS of the original ("from") image and the destination ("to")
@@ -35,6 +35,14 @@ def calc_pixmap(wcs_from, wcs_to, shape=None):
         ``numpy.ndarray`` order. When provided, it takes precedence over the
         ``wcs_from.array_shape`` property.
 
+    disable_bbox : {"to", "from", "both", "none"}, optional
+        Indicates whether to use or not to use the bounding box of either
+        (both) ``wcs_from`` or (and) ``wcs_to`` when computing pixel map. When
+        ``disable_bbox`` is "none", pixel coordinates outside of the bounding
+        box are set to `NaN` only if ``wcs_from`` or (and) ``wcs_to`` sets
+        world coordinates to NaN when input pixel coordinates are outside of
+        the bounding box.
+
     Returns
     -------
     pixmap : numpy.ndarray
@@ -57,35 +65,62 @@ def calc_pixmap(wcs_from, wcs_to, shape=None):
     If ``bounding_box`` is not available, a `ValueError` will be raised.
 
     """
+    if (bbox_from := getattr(wcs_from, "bounding_box", None)) is not None:
+        try:
+            # to avoid dependency on astropy just to check whether
+            # the bounding box is an instance of
+            # modeling.bounding_box.ModelBoundingBox, we try to
+            # directly use and bounding_box(order='F') and if it fails,
+            # fall back to converting the bounding box to a tuple
+            # (of intervals):
+            bbox_from = bbox_from.bounding_box(order='F')
+        except AttributeError:
+            bbox_from = tuple(bbox_from)
+
+    if (bbox_to := getattr(wcs_to, "bounding_box", None)) is not None:
+        try:
+            # to avoid dependency on astropy just to check whether
+            # the bounding box is an instance of
+            # modeling.bounding_box.ModelBoundingBox, we try to
+            # directly use and bounding_box(order='F') and if it fails,
+            # fall back to converting the bounding box to a tuple
+            # (of intervals):
+            bbox_to = bbox_to.bounding_box(order='F')
+        except AttributeError:
+            bbox_to = tuple(bbox_to)
+
     if shape is None:
         shape = wcs_from.array_shape
-        if shape is None:
-            if (bbox := getattr(wcs_from, "bounding_box", None)) is not None:
-                try:
-                    # to avoid dependency on astropy just to check whether
-                    # the bounding box is an instance of
-                    # modeling.bounding_box.ModelBoundingBox, we try to
-                    # directly use and bounding_box(order='F') and if it fails,
-                    # fall back to converting the bounding box to a tuple
-                    # (of intervals):
-                    bbox = bbox.bounding_box(order='F')
-                except AttributeError:
-                    bbox = tuple(bbox)
-
-                if (nd := np.ndim(bbox)) == 1:
-                    bbox = (bbox, )
-                if nd > 1:
-                    shape = tuple(
-                        int(math.ceil(lim[1] + 0.5)) for lim in bbox[::-1]
-                    )
+        if shape is None and bbox_from is not None:
+            if (nd := np.ndim(bbox_from)) == 1:
+                bbox_from = (bbox_from, )
+            if nd > 1:
+                shape = tuple(
+                    int(math.ceil(lim[1] + 0.5)) for lim in bbox_from[::-1]
+                )
 
     if shape is None:
         raise ValueError(
             'The "from" WCS must have pixel_shape property set.'
         )
 
     y, x = np.indices(shape, dtype=np.float64)
-    x, y = wcs_to.world_to_pixel_values(*wcs_from.pixel_to_world_values(x, y))
+
+    # temporarily disable the bounding box for the "from" WCS:
+    if disable_bbox in ["from", "both"] and bbox_from is not None:
+        wcs_from.bounding_box = None
+    if disable_bbox in ["to", "both"] and bbox_to is not None:
+        wcs_to.bounding_box = None
+    try:
+        x, y = wcs_to.world_to_pixel_values(
+            *wcs_from.pixel_to_world_values(x, y)
+        )
+    finally:
+        if bbox_from is not None:
+            wcs_from.bounding_box = bbox_from
+        if bbox_to is not None:
+            wcs_to.bounding_box = bbox_to
+
     pixmap = np.dstack([x, y])
     return pixmap