Fixes for bounding box units

CHANGES.rst

@@ -16,6 +16,9 @@
 - Refactor ``WCS`` to use a ``Pipeline`` base class which adds basic checks to ensure that the pipeline is valid. These
   include checking for duplicate frame names and that the last transform is ``None``. [#545]
 
+- Bugfix for ``WCS.invert`` and ``WCS.to_fits`` that prevented evaluation when the attached bounding box happened to have
+  units on its values. [#554]
+
 
 0.22.0 (2024-12-19)
 -------------------

gwcs/tests/test_wcs.py

@@ -1759,3 +1759,48 @@ def test_error_with_not_none_last():
         ValueError, match="The last step in the pipeline must have a None transform."
     ):
         wcs.WCS(pipeline)
+
+
+def test_bounding_box_with_units():
+    """
+    Test that the invert method works when a bounding box has units.
+    """
+
+    # GWCS that is adapted from its Getting Started.
+    shift_by_crpix = models.Shift(-(5 - 1) * u.pix) & models.Shift(-(5 - 1) * u.pix)
+    matrix = np.array(
+        [
+            [1.290551569736e-05, 5.9525007864732e-06],
+            [5.0226382102765e-06, -1.2644844123757e-05],
+        ]
+    )
+    rotation = models.AffineTransformation2D(matrix * u.deg, translation=[0, 0] * u.deg)
+    rotation.input_units_equivalencies = {
+        "x": u.pixel_scale(1 * (u.deg / u.pix)),
+        "y": u.pixel_scale(1 * (u.deg / u.pix)),
+    }
+    rotation.inverse = models.AffineTransformation2D(
+        np.linalg.inv(matrix) * u.pix, translation=[0, 0] * u.pix
+    )
+    rotation.inverse.input_units_equivalencies = {
+        "x": u.pixel_scale(1 * (u.pix / u.deg)),
+        "y": u.pixel_scale(1 * (u.pix / u.deg)),
+    }
+    tan = models.Pix2Sky_TAN()
+    celestial_rotation = models.RotateNative2Celestial(
+        3.581704851882 * u.deg, -30.39197867265 * u.deg, 180 * u.deg
+    )
+    det2sky = shift_by_crpix | rotation | tan | celestial_rotation
+    det2sky.name = "linear_transform"
+    detector_frame = cf.Frame2D(
+        name="detector", axes_names=("x", "y"), unit=(u.pix, u.pix)
+    )
+    sky_frame = cf.CelestialFrame(
+        reference_frame=coord.ICRS(), name="icrs", unit=(u.deg, u.deg)
+    )
+    pipeline = [(detector_frame, det2sky), (sky_frame, None)]
+    w_gwcs = wcs.WCS(pipeline)
+    w_gwcs.bounding_box = ((0, 8), (0, 10)) * u.pix  # x, y
+
+    w_gwcs.invert(4 * u.deg, 5 * u.deg)
+    w_gwcs.to_fits(bounding_box=([0, 100] * u.pix, [0, 100] * u.pix))

gwcs/wcs/_wcs.py

@@ -403,7 +403,8 @@ def outside_footprint(self, world_arrays):
                     max_ax = axis_range[~m].min()
                     outside = (coord > min_ax) & (coord < max_ax)
                 else:
-                    outside = (coord < min_ax) | (coord > max_ax)
+                    coord_ = self._remove_units_input([coord], self.output_frame)[0]
+                    outside = (coord_ < min_ax) | (coord_ > max_ax)
                 if np.any(outside):
                     if np.isscalar(coord):
                         coord = np.nan
@@ -1177,6 +1178,7 @@ def footprint(self, bounding_box=None, center=False, axis_type="all"):
         """
 
         def _order_clockwise(v):
+            v = [self._remove_units_input(vv, self.input_frame) for vv in v]
             return np.asarray(
                 [
                     [v[0][0], v[1][0]],
@@ -1556,6 +1558,11 @@ def _to_fits_sip(
         if bounding_box is None:
             bounding_box = self.bounding_box
 
+        first_bound = bounding_box[0][0]
+        if isinstance(first_bound, u.Quantity):
+            bounding_box = [
+                self._remove_units_input(bb, self.input_frame) for bb in bounding_box
+            ]
         bb_center = np.mean(bounding_box, axis=1)
 
         fixi_dict = {
@@ -1581,54 +1588,90 @@ def _to_fits_sip(
             crpix1 = crpix[0] - 1
             crpix2 = crpix[1] - 1
 
+        if isinstance(first_bound, u.Quantity):
+            crpix1 = u.Quantity(crpix1, first_bound.unit)
+            crpix2 = u.Quantity(crpix2, first_bound.unit)
+
         # check that the bounding box has some reasonable size:
         if (xmax - xmin) < 1 or (ymax - ymin) < 1:
             msg = "Bounding box is too small for fitting a SIP polynomial"
             raise ValueError(msg)
 
         lon, lat = transform(crpix1, crpix2)
+        pole = 180
+        if isinstance(lon, u.Quantity):
+            pole = u.Quantity(pole, lon.unit)
 
         # Now rotate to native system and deproject. Recall that transform
         # expects pixels in the original coordinate system, but the SIP
         # transform is relative to crpix coordinates, thus the initial shift.
         ntransform = (
             (Shift(crpix1) & Shift(crpix2))
             | transform
-            | RotateCelestial2Native(lon, lat, 180)
+            | RotateCelestial2Native(lon, lat, pole)
             | sky2pix_proj
         )
 
         # standard sampling:
-        u, v = make_sampling_grid(
-            npoints, tuple(bounding_box[k] for k in input_axes), crpix=[crpix1, crpix2]
+        crpix_ = [crpix1, crpix2]
+        if isinstance(crpix1, u.Quantity):
+            crpix_ = self._remove_units_input(crpix_, self.input_frame)
+        u_grid, v_grid = make_sampling_grid(
+            npoints, tuple(bounding_box[k] for k in input_axes), crpix=crpix_
         )
-        undist_x, undist_y = ntransform(u, v)
+        if isinstance(crpix1, u.Quantity):
+            u_grid = u.Quantity(u_grid, crpix1.unit)
+            v_grid = u.Quantity(v_grid, crpix2.unit)
+
+        undist_x, undist_y = ntransform(u_grid, v_grid)
 
         # Double sampling to check if sampling is sufficient.
         ud, vd = make_sampling_grid(
             2 * npoints,
             tuple(bounding_box[k] for k in input_axes),
-            crpix=[crpix1, crpix2],
+            crpix=crpix_,
         )
+        if isinstance(crpix1, u.Quantity):
+            ud = u.Quantity(ud, crpix1.unit)
+            vd = u.Quantity(vd, crpix2.unit)
+
         undist_xd, undist_yd = ntransform(ud, vd)
 
+        input_0 = 0
+        input_1 = 1
+        if isinstance(crpix1, u.Quantity):
+            input_0 = u.Quantity(0, crpix1.unit)
+            input_1 = u.Quantity(1, crpix1.unit)
+
         # Determine approximate pixel scale in order to compute error threshold
         # from the specified pixel error. Computed at the center of the array.
-        x0, y0 = ntransform(0, 0)
-        xx, xy = ntransform(1, 0)
-        yx, yy = ntransform(0, 1)
+        x0, y0 = ntransform(input_0, input_0)
+        xx, xy = ntransform(input_1, input_0)
+        yx, yy = ntransform(input_0, input_1)
         pixarea = np.abs((xx - x0) * (yy - y0) - (xy - y0) * (yx - x0))
         plate_scale = np.sqrt(pixarea)
 
+        plate_scale = (
+            plate_scale.value if isinstance(plate_scale, u.Quantity) else plate_scale
+        )
+        u_grid = u_grid.value if isinstance(u_grid, u.Quantity) else u_grid
+        v_grid = v_grid.value if isinstance(v_grid, u.Quantity) else v_grid
+        undist_x = undist_x.value if isinstance(undist_x, u.Quantity) else undist_x
+        undist_y = undist_y.value if isinstance(undist_y, u.Quantity) else undist_y
+        ud = ud.value if isinstance(ud, u.Quantity) else ud
+        vd = vd.value if isinstance(vd, u.Quantity) else vd
+        undist_xd = undist_xd.value if isinstance(undist_xd, u.Quantity) else undist_xd
+        undist_yd = undist_yd.value if isinstance(undist_yd, u.Quantity) else undist_yd
+
         # The fitting section.
         if verbose:
             sys.stdout.write("\nFitting forward SIP ...")
         fit_poly_x, fit_poly_y, max_resid = fit_2D_poly(
             degree,
             max_pix_error,
             plate_scale,
-            u,
-            v,
+            u_grid,
+            v_grid,
             undist_x,
             undist_y,
             ud,
@@ -1658,8 +1701,8 @@ def _to_fits_sip(
                 1,
                 U,
                 V,
-                u - U,
-                v - V,
+                u_grid - U,
+                v_grid - V,
                 Ud,
                 Vd,
                 ud - Ud,
@@ -1669,8 +1712,14 @@ def _to_fits_sip(
 
         # create header with WCS info:
         w = celestial_frame_to_wcs(frame.reference_frame, projection=projection)
-        w.wcs.crval = [lon, lat]
-        w.wcs.crpix = [crpix1 + 1, crpix2 + 1]
+        w.wcs.crval = [
+            lon.value if isinstance(lon, u.Quantity) else lon,
+            lat.value if isinstance(lat, u.Quantity) else lat,
+        ]
+        w.wcs.crpix = [
+            crpix1.value if isinstance(crpix1, u.Quantity) else crpix1 + 1,
+            crpix2.value if isinstance(crpix2, u.Quantity) else crpix2 + 1,
+        ]
         w.wcs.pc = cdmat if nlon < nlat else cdmat[::-1]
         w.wcs.set()
         hdr = w.to_header(True)