Fix inside/outside calc for inverse with bbox (#536)

CHANGES.rst

@@ -14,7 +14,7 @@
 
 - Add ``gwcs.examples`` module, based on the examples located in the testing ``conftest.py``. [#521]
 
-- Force ``bounding_box`` to always be returned as a ``F`` ordered box. [#522] 
+- Force ``bounding_box`` to always be returned as a ``F`` ordered box. [#522]
 
 - Move the bounding box attachment to the forward transform property. [#532]
 
@@ -25,6 +25,9 @@
 - Fixed a bug where evaluating the inverse transform did not
   respect the bounding box. [#498]
 
+- Improved reliability of inside/outside footprint computations when evaluating
+  inverse transform with bounding box. [#536]
+
 
 0.21.0 (2024-03-10)
 -------------------

gwcs/examples.py

@@ -36,7 +36,6 @@ def gwcs_2d_spatial_shift():
     A simple one step spatial WCS, in ICRS with a 1 and 2 px shift.
     """
     pipe = [(DETECTOR_2D_FRAME, MODEL_2D_SHIFT), (ICRC_SKY_FRAME, None)]
-
     return wcs.WCS(pipe)
 
 
@@ -185,6 +184,29 @@ def gwcs_simple_imaging_units():
     return wcs.WCS(pipeline)
 
 
+def gwcs_simple_imaging():
+    shift_by_crpix = models.Shift(-2048) & models.Shift(-1024)
+    matrix = np.array([[1.290551569736E-05, 5.9525007864732E-06],
+                       [5.0226382102765E-06 , -1.2644844123757E-05]])
+    rotation = models.AffineTransformation2D(matrix,
+                                             translation=[0, 0])
+    rotation.inverse = models.AffineTransformation2D(np.linalg.inv(matrix),
+                                                     translation=[0, 0])
+    tan = models.Pix2Sky_TAN()
+    celestial_rotation =  models.RotateNative2Celestial(5.63056810618,
+                                                        -72.05457184279,
+                                                        180)
+    det2sky = shift_by_crpix | rotation | tan | celestial_rotation
+    det2sky.name = "linear_transform"
+
+    detector_frame = cf.Frame2D(name="detector", axes_names=("x", "y"))
+    sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), name='icrs')
+    pipeline = [(detector_frame, det2sky),
+                (sky_frame, None)
+                ]
+    return wcs.WCS(pipeline)
+
+
 def gwcs_stokes_lookup():
     transform = models.Tabular1D([0, 1, 2, 3] * u.pix, [1, 2, 3, 4] * u.one,
                                  method="nearest", fill_value=np.nan, bounds_error=False)

gwcs/tests/conftest.py

@@ -61,6 +61,11 @@ def gwcs_simple_imaging_units():
     return examples.gwcs_simple_imaging_units()
 
 
+@pytest.fixture
+def gwcs_simple_imaging():
+    return examples.gwcs_simple_imaging()
+
+
 @pytest.fixture
 def gwcs_stokes_lookup():
     return examples.gwcs_stokes_lookup()

gwcs/tests/test_api.py

@@ -460,10 +460,11 @@ def test_world_to_pixel(gwcs_2d_spatial_shift, sky_ra_dec):
     assert_allclose(wcsobj.world_to_pixel(sky), wcsobj.invert(ra, dec, with_units=False))
 
 
-def test_world_to_array_index(gwcs_2d_spatial_shift, sky_ra_dec):
-    wcsobj = gwcs_2d_spatial_shift
+def test_world_to_array_index(gwcs_simple_imaging, sky_ra_dec):
+    wcsobj = gwcs_simple_imaging
     sky, ra, dec = sky_ra_dec
-    assert_allclose(wcsobj.world_to_array_index(sky), wcsobj.invert(ra, dec, with_units=False)[::-1])
+
+    assert_allclose(wcsobj.world_to_array_index(sky), wcsobj.invert(ra * u.deg, dec * u.deg, with_units=False)[::-1])
 
 
 def test_world_to_pixel_values(gwcs_2d_spatial_shift, sky_ra_dec):
@@ -473,12 +474,12 @@ def test_world_to_pixel_values(gwcs_2d_spatial_shift, sky_ra_dec):
     assert_allclose(wcsobj.world_to_pixel_values(sky), wcsobj.invert(ra, dec, with_units=False))
 
 
-def test_world_to_array_index_values(gwcs_2d_spatial_shift, sky_ra_dec):
-    wcsobj = gwcs_2d_spatial_shift
+def test_world_to_array_index_values(gwcs_simple_imaging, sky_ra_dec):
+    wcsobj = gwcs_simple_imaging
     sky, ra, dec = sky_ra_dec
 
     assert_allclose(wcsobj.world_to_array_index_values(sky),
-                    wcsobj.invert(ra, dec, with_units=False)[::-1])
+                    wcsobj.invert(ra * u.deg, dec * u.deg, with_units=False)[::-1])
 
 
 def test_ndim_str_frames(gwcs_with_frames_strings):

gwcs/wcs.py

@@ -493,6 +493,7 @@ def outside_footprint(self, world_arrays):
         world_arrays = list(world_arrays)
 
         axes_types = set(self.output_frame.axes_type)
+        axes_phys_types = self.world_axis_physical_types
         footprint = self.footprint()
         not_numerical = False
         if not utils.isnumerical(world_arrays[0]):
@@ -501,14 +502,25 @@ def outside_footprint(self, world_arrays):
         for axtyp in axes_types:
             ind = np.asarray((np.asarray(self.output_frame.axes_type) == axtyp))
 
-            for idim, coord in enumerate(world_arrays):
+            for idim, (coord, phys) in enumerate(zip(world_arrays, axes_phys_types)):
                 coord = _tofloat(coord)
                 if np.asarray(ind).sum() > 1:
                     axis_range = footprint[:, idim]
                 else:
                     axis_range = footprint
                 range = [axis_range.min(), axis_range.max()]
-                outside = (coord < range[0]) | (coord > range[1])
+
+                if (axtyp == 'SPATIAL' and str(phys).endswith((".ra", ".lon"))
+                    and range[1] - range[0] > 180):
+                        # most likely this coordinate is wrapped at 360
+                        d = np.mean(range)
+                        range = [
+                            axis_range[axis_range < d].max(),
+                            axis_range[axis_range > d].min()
+                        ]
+                        outside = (coord >= range[0]) & (coord < range[1])
+                else:
+                    outside = (coord < range[0]) | (coord > range[1])
                 if np.any(outside):
                     if np.isscalar(coord):
                         coord = np.nan