From 68fcfa83c9c1061d68c303b87ead83c0de187f03 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Wed, 18 Dec 2024 09:17:45 -0500 Subject: [PATCH] Fix inside/outside calc for inverse with bbox (#536) --- CHANGES.rst | 5 ++++- gwcs/examples.py | 24 +++++++++++++++++++++++- gwcs/tests/conftest.py | 5 +++++ gwcs/tests/test_api.py | 13 +++++++------ gwcs/wcs.py | 16 ++++++++++++++-- 5 files changed, 53 insertions(+), 10 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index d308cb23..79ba331a 100644 --- a/CHANGES.rst +++ b/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) ------------------- diff --git a/gwcs/examples.py b/gwcs/examples.py index 61fc9387..cd5713ae 100644 --- a/gwcs/examples.py +++ b/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) diff --git a/gwcs/tests/conftest.py b/gwcs/tests/conftest.py index 92979383..3bf243e5 100644 --- a/gwcs/tests/conftest.py +++ b/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() diff --git a/gwcs/tests/test_api.py b/gwcs/tests/test_api.py index a676d3fb..cdb49c76 100644 --- a/gwcs/tests/test_api.py +++ b/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): diff --git a/gwcs/wcs.py b/gwcs/wcs.py index a52dc06e..e0cbf191 100644 --- a/gwcs/wcs.py +++ b/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