From 6b48b174d8a06cbc28f9c3e8bbdf4db09fe977fe Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Sun, 12 May 2024 15:25:19 -0400 Subject: [PATCH 01/19] fix evaluating an inverse with bounding_box --- gwcs/tests/conftest.py | 54 ++++++++++++++++++++++++++++++++++ gwcs/tests/test_api_slicing.py | 4 +-- gwcs/wcs.py | 30 +++++++++++++------ 3 files changed, 77 insertions(+), 11 deletions(-) diff --git a/gwcs/tests/conftest.py b/gwcs/tests/conftest.py index 83fb216b..62a0a2d5 100644 --- a/gwcs/tests/conftest.py +++ b/gwcs/tests/conftest.py @@ -5,6 +5,31 @@ from .. import examples from .. import geometry +import numpy as np + +import astropy.units as u +from astropy.time import Time +from astropy import coordinates as coord +from astropy.modeling import models + +from gwcs import coordinate_frames as cf +from gwcs import spectroscopy as sp +from gwcs import wcs +from gwcs import geometry + +# frames +detector_1d = cf.CoordinateFrame(name='detector', axes_order=(0,), naxes=1, axes_type="detector") +detector_2d = cf.Frame2D(name='detector', axes_order=(0, 1)) +icrs_sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), + axes_order=(0, 1)) + +freq_frame = cf.SpectralFrame(name='freq', unit=u.Hz, axes_order=(0, )) +wave_frame = cf.SpectralFrame(name='wave', unit=u.m, axes_order=(2, ), + axes_names=('lambda', )) + +# transforms +model_2d_shift = models.Shift(1) & models.Shift(2) +model_1d_scale = models.Scale(2) @pytest.fixture @@ -141,3 +166,32 @@ def spher_to_cart(): @pytest.fixture def cart_to_spher(): return geometry.CartesianToSpherical() + + +@pytest.fixture +def gwcs_simple_imaging_no_units(): + 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"), + 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 = wcs.WCS(pipeline) + w.bounding_box = ((2, 100), (5, 500)) + return w diff --git a/gwcs/tests/test_api_slicing.py b/gwcs/tests/test_api_slicing.py index ba0a9b54..626a0e2a 100644 --- a/gwcs/tests/test_api_slicing.py +++ b/gwcs/tests/test_api_slicing.py @@ -262,8 +262,8 @@ def test_celestial_slice(gwcs_3d_galactic_spectral): assert_allclose(wcs.pixel_to_world_values(39, 44), (10.24, 20, 25)) assert_allclose(wcs.array_index_to_world_values(44, 39), (10.24, 20, 25)) - assert_allclose(wcs.world_to_pixel_values(12.4, 20, 25), (39., 44.)) - assert_equal(wcs.world_to_array_index_values(12.4, 20, 25), (44, 39)) + assert_allclose(wcs.world_to_pixel_values(10.24, 20, 25), (39., 44.)) + assert_equal(wcs.world_to_array_index_values(10.24, 20, 25), (44, 39)) assert_equal(wcs.pixel_bounds, [(-2, 45), (5, 50)]) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 2ae714b2..d76e1391 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -396,10 +396,7 @@ def in_image(self, *args, **kwargs): and `False` if input is outside the footprint. """ - kwargs['with_bounding_box'] = True - kwargs['fill_value'] = np.nan - - coords = self.invert(*args, **kwargs) + coords = self.invert(*args, with_bounding_box=False, **kwargs) result = np.isfinite(coords) if self.input_frame.naxes > 1: @@ -480,11 +477,8 @@ def invert(self, *args, **kwargs): except (NotImplementedError, KeyError): args = utils.get_values(self.output_frame.unit, *args) - if 'with_bounding_box' not in kwargs: - kwargs['with_bounding_box'] = True - - if 'fill_value' not in kwargs: - kwargs['fill_value'] = np.nan + with_bounding_box = kwargs.pop('with_bounding_box', True) + fill_value = kwargs.pop('fill_value', np.nan) try: # remove iterative inverse-specific keyword arguments: @@ -493,6 +487,24 @@ def invert(self, *args, **kwargs): except (NotImplementedError, KeyError): result = self.numerical_inverse(*args, **kwargs, with_units=with_units) + if with_bounding_box and self.bounding_box is not None: + bbox = self.bounding_box + if self.input_frame.naxes == 1: + result = [result] + input_shape = self.backward_transform.input_shape(result) + valid_inputs, valid_index, all_out = bbox.prepare_inputs(input_shape, result) + if all_out: + if self.input_frame.naxes == 1: + return bbox._all_out_output(input_shape, fill_value)[0][0] + else: + return bbox._all_out_output(input_shape, fill_value)[0] + else: + if self.input_frame.naxes == 1: + valid_inputs = valid_inputs[0] + result = bbox.prepare_outputs(valid_inputs, valid_index, input_shape, fill_value)[0] + else: + result = tuple(bbox.prepare_outputs(valid_inputs, valid_index, input_shape, fill_value)) + if with_units and self.input_frame: if self.input_frame.naxes == 1: return self.input_frame.coordinates(result) From 55af60bfb8ff0b89645cd453db326ce53f7387e5 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Tue, 14 May 2024 12:35:18 -0400 Subject: [PATCH 02/19] add a change log --- CHANGES.rst | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index cca5803e..d308cb23 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -11,6 +11,7 @@ - Add support for compound bounding boxes and ignored bounding box entries. [#519] + - 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] @@ -19,8 +20,11 @@ - Adjust ``world_to_array_index_values`` to round to integer coordinates as specified by APE 14. [#525] -- Add warning filter to asdf extension to prevent the ``bounding_box`` order warning for gwcs - objects originating from a file. [#526] +- Add warning filter to asdf extension to prevent the ``bounding_box`` order warning for gwcs objects originating from a file. [#526] + +- Fixed a bug where evaluating the inverse transform did not + respect the bounding box. [#498] + 0.21.0 (2024-03-10) ------------------- From 646f0c1ca66cf83efdb58fe8837e86019e4035b4 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Wed, 29 May 2024 07:42:52 -0400 Subject: [PATCH 03/19] test --- gwcs/wcs.py | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index d76e1391..1027dd39 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -472,6 +472,7 @@ def invert(self, *args, **kwargs): if self.output_frame.naxes == 1: args = [args] try: + # uses_quantity constructs the backward_transform if not self.backward_transform.uses_quantity: args = utils.get_values(self.output_frame.unit, *args) except (NotImplementedError, KeyError): @@ -483,27 +484,38 @@ def invert(self, *args, **kwargs): try: # remove iterative inverse-specific keyword arguments: akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS} - result = self.backward_transform(*args, **akwargs) + btrans = self.backward_transform + result = btrans(*args, **akwargs) + #result_shape = btrans.input_shape(result) except (NotImplementedError, KeyError): + btrans = None result = self.numerical_inverse(*args, **kwargs, with_units=with_units) + #result_shape = result[0].shape if with_bounding_box and self.bounding_box is not None: bbox = self.bounding_box if self.input_frame.naxes == 1: result = [result] - input_shape = self.backward_transform.input_shape(result) - valid_inputs, valid_index, all_out = bbox.prepare_inputs(input_shape, result) + if btrans is not None: + result_shape = btrans.input_shape(result) + else: + # numerical_inversse was run - > 2 outputs + if np.isscalar(result[0]): + result_shape = () + else: + result_shape = result[0].shape + valid_inputs, valid_index, all_out = bbox.prepare_inputs(result_shape, result) if all_out: if self.input_frame.naxes == 1: - return bbox._all_out_output(input_shape, fill_value)[0][0] + return bbox._all_out_output(result_shape, fill_value)[0][0] else: - return bbox._all_out_output(input_shape, fill_value)[0] + return bbox._all_out_output(result_shape, fill_value)[0] else: if self.input_frame.naxes == 1: valid_inputs = valid_inputs[0] - result = bbox.prepare_outputs(valid_inputs, valid_index, input_shape, fill_value)[0] + result = bbox.prepare_outputs(valid_inputs, valid_index, result_shape, fill_value)[0] else: - result = tuple(bbox.prepare_outputs(valid_inputs, valid_index, input_shape, fill_value)) + result = tuple(bbox.prepare_outputs(valid_inputs, valid_index, result_shape, fill_value)) if with_units and self.input_frame: if self.input_frame.naxes == 1: From 3f55267993b47a47d79fc686fb6a2ded6e586f0a Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Wed, 29 May 2024 15:41:22 -0400 Subject: [PATCH 04/19] add a fix for running world_to_pixel when the backward transform is not defined --- gwcs/api.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/gwcs/api.py b/gwcs/api.py index 2b1ca22e..981a0ce9 100644 --- a/gwcs/api.py +++ b/gwcs/api.py @@ -132,7 +132,13 @@ def world_to_pixel_values(self, *world_arrays): be returned in the ``(x, y)`` order, where for an image, ``x`` is the horizontal coordinate and ``y`` is the vertical coordinate. """ - world_arrays = self._add_units_input(world_arrays, self.backward_transform, self.output_frame) + try: + backward_transform = self.backward_transform + world_arrays = self._add_units_input(world_arrays, + backward_transform, + self.output_frame) + except NotImplementedError: + pass result = self.invert(*world_arrays, with_units=False) @@ -296,6 +302,10 @@ def _sanitize_pixel_inputs(self, *pixel_arrays): return pixels + def _sanitize_world_inputs(self, *world_arrays): + world_coord = [] + + def pixel_to_world(self, *pixel_arrays): """ Convert pixel values to world coordinates. From e1556e9c400bbd8b21fc6ae9fb0773ae0a94d904 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Wed, 29 May 2024 15:53:11 -0400 Subject: [PATCH 05/19] set detect_divergence to False as a test --- gwcs/wcs.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 1027dd39..068cd757 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -525,8 +525,8 @@ def invert(self, *args, **kwargs): else: return result - def numerical_inverse(self, *args, tolerance=1e-5, maxiter=50, adaptive=True, - detect_divergence=True, quiet=True, with_bounding_box=True, + def numerical_inverse(self, *args, tolerance=1e-4, maxiter=30, adaptive=True, + detect_divergence=False, quiet=True, with_bounding_box=True, fill_value=np.nan, with_units=False, **kwargs): """ Invert coordinates from output frame to input frame using numerical From 631639702e1bba5d9e2baba0e8645280a62aa1dc Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Wed, 29 May 2024 16:14:30 -0400 Subject: [PATCH 06/19] fix bbox with units --- gwcs/wcs.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 068cd757..66279089 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -504,6 +504,13 @@ def invert(self, *args, **kwargs): result_shape = () else: result_shape = result[0].shape + if self.input_frame.naxes > 1: + first_res = result[0] + if not utils.isnumerical(first_res): + result = [i.value for i in result] + else: + if not utils.isnumerical(result): + result = result.value valid_inputs, valid_index, all_out = bbox.prepare_inputs(result_shape, result) if all_out: if self.input_frame.naxes == 1: From 8abeeb7241b676a795f101d22aa4c87d7ec35ba1 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Sun, 2 Jun 2024 16:27:20 -0400 Subject: [PATCH 07/19] revert parameter value changes --- gwcs/tests/test_wcs.py | 4 ++++ gwcs/wcs.py | 4 ++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 5c9af093..1baad46f 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -1199,6 +1199,7 @@ def test_iter_inv(): *w(x, y), adaptive=True, detect_divergence=True, + tolerance=1e-4, maxiter=50, quiet=False ) assert np.allclose((x, y), (xp, yp)) @@ -1218,6 +1219,7 @@ def test_iter_inv(): xp, yp = w.numerical_inverse( *w(x, y), adaptive=True, + tolerance=1e-5, maxiter=50, detect_divergence=False, quiet=False ) @@ -1252,6 +1254,7 @@ def test_iter_inv(): xp, yp = w.numerical_inverse( *w(x, y, with_bounding_box=False), adaptive=False, + tolerance=1e-5, maxiter=50, detect_divergence=True, quiet=False, with_bounding_box=False @@ -1265,6 +1268,7 @@ def test_iter_inv(): xp, yp = w.numerical_inverse( *w(x, y, with_bounding_box=False), adaptive=False, + tolerance=1e-5, maxiter=50, detect_divergence=True, quiet=False, with_bounding_box=False diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 66279089..89444e14 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -532,8 +532,8 @@ def invert(self, *args, **kwargs): else: return result - def numerical_inverse(self, *args, tolerance=1e-4, maxiter=30, adaptive=True, - detect_divergence=False, quiet=True, with_bounding_box=True, + def numerical_inverse(self, *args, tolerance=1e-5, maxiter=30, adaptive=True, + detect_divergence=True, quiet=True, with_bounding_box=True, fill_value=np.nan, with_units=False, **kwargs): """ Invert coordinates from output frame to input frame using numerical From b01a6b45b59f70f5eea07038ab019f2f7d79bf6e Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Sat, 8 Jun 2024 09:00:52 -0400 Subject: [PATCH 08/19] fix --- gwcs/tests/test_api.py | 17 ++---- gwcs/tests/test_coordinate_systems.py | 6 +-- gwcs/tests/test_utils.py | 6 +++ gwcs/utils.py | 11 ++-- gwcs/wcs.py | 76 ++++++++++++--------------- 5 files changed, 51 insertions(+), 65 deletions(-) diff --git a/gwcs/tests/test_api.py b/gwcs/tests/test_api.py index f80c8d8a..a676d3fb 100644 --- a/gwcs/tests/test_api.py +++ b/gwcs/tests/test_api.py @@ -63,7 +63,6 @@ def wcs_ndim_types_units(request): @fixture_all_wcses def test_lowlevel_types(wcsobj): - pytest.importorskip("typeguard") try: # Skip this on older versions of astropy where it dosen't exist. from astropy.wcs.wcsapi.tests.utils import validate_low_level_wcs_types @@ -236,12 +235,12 @@ def test_world_axis_object_classes_4d(gwcs_4d_identity_units): def _compare_frame_output(wc1, wc2): if isinstance(wc1, coord.SkyCoord): assert isinstance(wc1.frame, type(wc2.frame)) - assert u.allclose(wc1.spherical.lon, wc2.spherical.lon) - assert u.allclose(wc1.spherical.lat, wc2.spherical.lat) - assert u.allclose(wc1.spherical.distance, wc2.spherical.distance) + assert u.allclose(wc1.spherical.lon, wc2.spherical.lon, equal_nan=True) + assert u.allclose(wc1.spherical.lat, wc2.spherical.lat, equal_nan=True) + assert u.allclose(wc1.spherical.distance, wc2.spherical.distance, equal_nan=True) elif isinstance(wc1, u.Quantity): - assert u.allclose(wc1, wc2) + assert u.allclose(wc1, wc2, equal_nan=True) elif isinstance(wc1, time.Time): assert u.allclose((wc1 - wc2).to(u.s), 0*u.s) @@ -258,12 +257,6 @@ def _compare_frame_output(wc1, wc2): @fixture_all_wcses def test_high_level_wrapper(wcsobj, request): - if request.node.callspec.params['wcsobj'] in ('gwcs_4d_identity_units', 'gwcs_stokes_lookup'): - pytest.importorskip("astropy", minversion="4.0dev0") - - # Remove the bounding box because the type test is a little broken with the - # bounding box. - del wcsobj._pipeline[0].transform.bounding_box hlvl = HighLevelWCSWrapper(wcsobj) @@ -286,8 +279,6 @@ def test_high_level_wrapper(wcsobj, request): def test_stokes_wrapper(gwcs_stokes_lookup): - pytest.importorskip("astropy", minversion="4.0dev0") - hlvl = HighLevelWCSWrapper(gwcs_stokes_lookup) pixel_input = [0, 1, 2, 3] diff --git a/gwcs/tests/test_coordinate_systems.py b/gwcs/tests/test_coordinate_systems.py index 967657f8..f1e34cb9 100644 --- a/gwcs/tests/test_coordinate_systems.py +++ b/gwcs/tests/test_coordinate_systems.py @@ -190,7 +190,7 @@ def test_temporal_relative(): assert a[1] == Time("2018-01-01T00:00:00") + 20 * u.s -@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher") +#@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher") def test_temporal_absolute(): t = cf.TemporalFrame(reference_frame=Time([], format='isot')) assert t.coordinates("2018-01-01T00:00:00") == Time("2018-01-01T00:00:00") @@ -240,7 +240,7 @@ def test_coordinate_to_quantity_spectral(inp): (Time("2011-01-01T00:00:10"),), (10 * u.s,) ]) -@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher.") +#@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher.") def test_coordinate_to_quantity_temporal(inp): temp = cf.TemporalFrame(reference_frame=Time("2011-01-01T00:00:00"), unit=u.s) @@ -325,7 +325,7 @@ def test_coordinate_to_quantity_frame_2d(): assert_quantity_allclose(output, exp) -@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher.") +#@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher.") def test_coordinate_to_quantity_error(): frame = cf.Frame2D(unit=(u.one, u.arcsec)) with pytest.raises(ValueError): diff --git a/gwcs/tests/test_utils.py b/gwcs/tests/test_utils.py index e69ec536..748b1472 100644 --- a/gwcs/tests/test_utils.py +++ b/gwcs/tests/test_utils.py @@ -6,6 +6,8 @@ from astropy import units as u from astropy import coordinates as coord from astropy.modeling import models +from astropy import table + from astropy.tests.helper import assert_quantity_allclose import pytest from numpy.testing import assert_allclose @@ -104,6 +106,10 @@ def test_isnumerical(): assert gwutils.isnumerical(np.array(0, dtype='>f8')) assert gwutils.isnumerical(np.array(0, dtype='>i4')) + # check a table column + t = table.Table(data=[[1,2,3], [4,5,6]], names=['x', 'y']) + assert not gwutils.isnumerical(t['x']) + def test_get_values(): args = 2 * u.cm diff --git a/gwcs/utils.py b/gwcs/utils.py index 104558cf..dcae0558 100644 --- a/gwcs/utils.py +++ b/gwcs/utils.py @@ -12,6 +12,7 @@ from astropy import coordinates as coords from astropy import units as u from astropy.time import Time, TimeDelta +from astropy import table from astropy.wcs import Celprm @@ -470,14 +471,12 @@ def isnumerical(val): Determine if a value is numerical (number or np.array of numbers). """ isnum = True - if isinstance(val, coords.SkyCoord): - isnum = False - elif isinstance(val, u.Quantity): - isnum = False - elif isinstance(val, (Time, TimeDelta)): + astropy_types=(coords.SkyCoord, u.Quantity, Time, TimeDelta, table.Column, table.Row) + if isinstance(val, astropy_types): isnum = False elif (isinstance(val, np.ndarray) and not np.issubdtype(val.dtype, np.floating) - and not np.issubdtype(val.dtype, np.integer)): + and not np.issubdtype(val.dtype, np.integer) + ): isnum = False return isnum diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 89444e14..ee39c43d 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -467,62 +467,33 @@ def invert(self, *args, **kwargs): """ with_units = kwargs.pop('with_units', False) + try: + btrans = self.backward_transform + except NotImplementedError: + btrans = None + if not utils.isnumerical(args[0]): + # convert astropy objects to numbers and arrays args = self.output_frame.coordinate_to_quantity(*args) if self.output_frame.naxes == 1: args = [args] - try: - # uses_quantity constructs the backward_transform - if not self.backward_transform.uses_quantity: - args = utils.get_values(self.output_frame.unit, *args) - except (NotImplementedError, KeyError): - args = utils.get_values(self.output_frame.unit, *args) + # if the transform does not use units, getthe numerical values + if btrans is not None and not btrans.uses_quantity: + args = utils.get_values(self.output_frame.unit, *args) + with_bounding_box = kwargs.pop('with_bounding_box', True) fill_value = kwargs.pop('fill_value', np.nan) - try: - # remove iterative inverse-specific keyword arguments: + if btrans is not None: akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS} - btrans = self.backward_transform result = btrans(*args, **akwargs) - #result_shape = btrans.input_shape(result) - except (NotImplementedError, KeyError): - btrans = None + else: result = self.numerical_inverse(*args, **kwargs, with_units=with_units) - #result_shape = result[0].shape + # deal with values outside the bounding box if with_bounding_box and self.bounding_box is not None: - bbox = self.bounding_box - if self.input_frame.naxes == 1: - result = [result] - if btrans is not None: - result_shape = btrans.input_shape(result) - else: - # numerical_inversse was run - > 2 outputs - if np.isscalar(result[0]): - result_shape = () - else: - result_shape = result[0].shape - if self.input_frame.naxes > 1: - first_res = result[0] - if not utils.isnumerical(first_res): - result = [i.value for i in result] - else: - if not utils.isnumerical(result): - result = result.value - valid_inputs, valid_index, all_out = bbox.prepare_inputs(result_shape, result) - if all_out: - if self.input_frame.naxes == 1: - return bbox._all_out_output(result_shape, fill_value)[0][0] - else: - return bbox._all_out_output(result_shape, fill_value)[0] - else: - if self.input_frame.naxes == 1: - valid_inputs = valid_inputs[0] - result = bbox.prepare_outputs(valid_inputs, valid_index, result_shape, fill_value)[0] - else: - result = tuple(bbox.prepare_outputs(valid_inputs, valid_index, result_shape, fill_value)) + result = self.out_of_bounds(result) if with_units and self.input_frame: if self.input_frame.naxes == 1: @@ -532,6 +503,25 @@ def invert(self, *args, **kwargs): else: return result + def out_of_bounds(self, pixel_arrays): + if np.isscalar(pixel_arrays) or self.input_frame.naxes == 1: + pixel_arrays = [pixel_arrays] + + pixel_arrays = list(pixel_arrays) + bbox = self.bounding_box + for idim, pix in enumerate(pixel_arrays): + outside = (pix < bbox[idim][0]) | (pix > bbox[idim][1]) + if np.any(outside): + if np.isscalar(pix): + pixel_arrays[idim] = np.nan + else: + pix = pixel_arrays[idim].astype(float, copy=True) + pix[outside] = np.nan + pixel_arrays[idim] = pix + if self.input_frame.naxes == 1: + pixel_arrays = pixel_arrays[0] + return pixel_arrays + def numerical_inverse(self, *args, tolerance=1e-5, maxiter=30, adaptive=True, detect_divergence=True, quiet=True, with_bounding_box=True, fill_value=np.nan, with_units=False, **kwargs): From 3cff9b6fcb83014e8691f0aed752fc9b42fdb4f4 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Wed, 21 Aug 2024 06:55:01 -0400 Subject: [PATCH 09/19] save work --- gwcs/wcs.py | 83 +++++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 71 insertions(+), 12 deletions(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index ee39c43d..b58a400f 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -471,8 +471,9 @@ def invert(self, *args, **kwargs): btrans = self.backward_transform except NotImplementedError: btrans = None - + print(f"args[0], {args[0]}") if not utils.isnumerical(args[0]): + print(f"args1, {args}") # convert astropy objects to numbers and arrays args = self.output_frame.coordinate_to_quantity(*args) if self.output_frame.naxes == 1: @@ -481,19 +482,23 @@ def invert(self, *args, **kwargs): # if the transform does not use units, getthe numerical values if btrans is not None and not btrans.uses_quantity: args = utils.get_values(self.output_frame.unit, *args) - + with_bounding_box = kwargs.pop('with_bounding_box', True) fill_value = kwargs.pop('fill_value', np.nan) + akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS} + print(f"args, {args}") + if with_bounding_box and self.bounding_box is not None: + result = self.outside_footprint(args) if btrans is not None: - akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS} + #akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS} result = btrans(*args, **akwargs) else: result = self.numerical_inverse(*args, **kwargs, with_units=with_units) # deal with values outside the bounding box if with_bounding_box and self.bounding_box is not None: - result = self.out_of_bounds(result) + result = self.out_of_bounds(result, fill_value=fill_value) if with_units and self.input_frame: if self.input_frame.naxes == 1: @@ -502,8 +507,50 @@ def invert(self, *args, **kwargs): return self.input_frame.coordinates(*result) else: return result + + def outside_footprint(self, world_arrays): + # for axis in world_arrays: + # if np.isscalar(axis): + # world_arrays = np.asarray(list(world_arrays), dtype = np.float64) + # world_arrays = [world_arrays] + #print('axis', axis, world_arrays) + #if np.isscalar(axis) or self.output_frame.naxes == 1: + if self.output_frame.naxes == 1: + #print('axis', axis) + #axis = float(axis) + world_arrays = [world_arrays] + #print('axis', axis) + print('world_arrays1', world_arrays) + #world_arrays = np.asarray(list(world_arrays))#, dtype = np.float64) + print('world_arrays2', world_arrays) + + axes_types = set(self.output_frame.axes_type) + for axtyp in axes_types: + footprint = self.footprint(axis_type=axtyp) + + ind = np.asarray((np.asarray(self.output_frame.axes_type) == axtyp)) + #print('ind', ind) + + for idim, coord in enumerate(world_arrays[ind]): + #print('footprint', footprint) + if np.asarray(ind).sum() > 1: + axis_range = footprint[:, idim] + else: + axis_range = footprint + range = [axis_range.min(), axis_range.max()] + #print('idim', idim, coord, range) + outside = (coord < range[0]) | (coord > range[1]) + if np.any(outside): + if np.isscalar(coord): + coord = np.nan + else: + coord[outside] = np.nan + world_arrays[idim] = coord + + return world_arrays - def out_of_bounds(self, pixel_arrays): + + def out_of_bounds(self, pixel_arrays, fill_value=np.nan): if np.isscalar(pixel_arrays) or self.input_frame.naxes == 1: pixel_arrays = [pixel_arrays] @@ -704,7 +751,7 @@ def numerical_inverse(self, *args, tolerance=1e-5, maxiter=30, adaptive=True, >>> import numpy as np >>> filename = get_pkg_data_filename('data/nircamwcs.asdf', package='gwcs.tests') - >>> with asdf.open(filename, memmap=False, lazy_load=False, ignore_missing_extensions=True) as af: + >>> with asdf.open(filename, copy_arrays=True, lazy_load=False, ignore_missing_extensions=True) as af: ... w = af.tree['wcs'] >>> ra, dec = w([1,2,3], [1,1,1]) @@ -1426,23 +1473,33 @@ def footprint(self, bounding_box=None, center=False, axis_type="all"): """ def _order_clockwise(v): + # if self.input_frame.naxes == 1: + # bb = self.bounding_box.bounding_box() + # if isinstance(bb[0], u.Quantity): + # bb = [v.value for v in bb] * bb[0].unit + # return (bb,) return np.asarray([[v[0][0], v[1][0]], [v[0][0], v[1][1]], [v[0][1], v[1][1]], [v[0][1], v[1][0]]]).T if bounding_box is None: if self.bounding_box is None: raise TypeError("Need a valid bounding_box to compute the footprint.") - bb = self.bounding_box + bb = self.bounding_box.bounding_box() else: bb = bounding_box all_spatial = all([t.lower() == "spatial" for t in self.output_frame.axes_type]) - - if all_spatial: + if self.output_frame.naxes == 1: + if isinstance(bb[0], u.Quantity): + bb = np.asarray([b.value for b in bb]) * bb[0].unit + vertices = (bb,) + elif all_spatial: vertices = _order_clockwise(bb) else: vertices = np.array(list(itertools.product(*bb))).T + # workaround an issue with bbox with quantity, interval needs to be a cquantity, not a list of quantities + # strip units if center: vertices = utils._toindex(vertices) @@ -1451,19 +1508,21 @@ def _order_clockwise(v): axis_type = axis_type.lower() if axis_type == 'spatial' and all_spatial: return result.T - + if axis_type != "all": axtyp_ind = np.array([t.lower() for t in self.output_frame.axes_type]) == axis_type if not axtyp_ind.any(): raise ValueError('This WCS does not have axis of type "{}".'.format(axis_type)) - result = np.asarray([(r.min(), r.max()) for r in result[axtyp_ind]]) + if len(axtyp_ind) > 1: + result = np.asarray([(r.min(), r.max()) for r in result[axtyp_ind]]) if axis_type == "spatial": result = _order_clockwise(result) else: result.sort() result = np.squeeze(result) - + if self.output_frame.naxes == 1: + return np.array([result]).T return result.T def fix_inputs(self, fixed): From 49370289e7fe1f2f5ff397230ce38cf658d38f1e Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Wed, 21 Aug 2024 06:55:37 -0400 Subject: [PATCH 10/19] add bbox tests --- gwcs/tests/test_bounding_box.py | 121 ++++++++++++++++++++++++++++++++ 1 file changed, 121 insertions(+) create mode 100644 gwcs/tests/test_bounding_box.py diff --git a/gwcs/tests/test_bounding_box.py b/gwcs/tests/test_bounding_box.py new file mode 100644 index 00000000..5fc13f79 --- /dev/null +++ b/gwcs/tests/test_bounding_box.py @@ -0,0 +1,121 @@ +import numpy as np +from numpy.testing import assert_array_equal, assert_allclose + +from astropy import units as u + +import pytest + + +x = [-1, 2, 4, 13] +y = [np.nan, np.nan, 4, np.nan] +y1 = [np.nan, np.nan, 4, np.nan] + + +@pytest.mark.parametrize((("input", "output")), [((2, 4), (2, 4)), + ((100, 200), (np.nan, np.nan)), + ((x, x), (y, y)) + ]) +def test_2d_spatial(gwcs_2d_spatial_shift, input, output): + w = gwcs_2d_spatial_shift + w.bounding_box = ((-.5, 21), (4, 12)) + + assert_array_equal(w.invert(*w(*input)), output) + assert_array_equal(w.world_to_pixel_values(*w.pixel_to_world_values(*input)), output) + assert_array_equal(w.world_to_pixel(w.pixel_to_world(*input)), output) + + +@pytest.mark.parametrize((("input", "output")), [((2, 4), (2, 4)), + ((100, 200), (np.nan, np.nan)), + ((x, x), (y, y)) + ]) +def test_2d_spatial_coordinate(gwcs_2d_quantity_shift, input, output): + w = gwcs_2d_quantity_shift + w.bounding_box = ((-.5, 21), (4, 12)) + + assert_array_equal(w.invert(*w(*input)), output) + assert_array_equal(w.world_to_pixel_values(*w.pixel_to_world_values(*input)), output) + assert_array_equal(w.world_to_pixel(*w.pixel_to_world(*input)), output) + + +@pytest.mark.parametrize((("input", "output")), [((2, 4), (2, 4)), + ((100, 200), (np.nan, np.nan)), + ((x, x), (y, y)) + ]) +def test_2d_spatial_coordinate_reordered(gwcs_2d_spatial_reordered, input, output): + w = gwcs_2d_spatial_reordered + w.bounding_box = ((-.5, 21), (4, 12)) + + assert_array_equal(w.invert(*w(*input)), output) + assert_array_equal(w.world_to_pixel_values(*w.pixel_to_world_values(*input)), output) + assert_array_equal(w.world_to_pixel(w.pixel_to_world(*input)), output) + + +@pytest.mark.parametrize((("input", "output")), [(2, 2), + ((10, 200), (10, np.nan)), + (x, (np.nan, 2, 4, 13)) + ]) +def test_1d_freq(gwcs_1d_freq, input, output): + w = gwcs_1d_freq + w.bounding_box = (-.5, 21) + print(f"input {input}, {output}") + assert_array_equal(w.invert(w(input)), output) + assert_array_equal(w.world_to_pixel_values(w.pixel_to_world_values(input)), output) + assert_array_equal(w.world_to_pixel(w.pixel_to_world(input)), output) + + +@pytest.mark.parametrize((("input", "output")), [((2, 4, 5), (2, 4, 5)), + ((100, 200, 5), (np.nan, np.nan, np.nan)), + ((x, x, x), (y1, y1, y1)) + ]) +def test_3d_spatial_wave(gwcs_3d_spatial_wave, input, output): + w = gwcs_3d_spatial_wave + w.bounding_box = ((-.5, 21), (4, 12), (3, 21)) + + assert_array_equal(w.invert(*w(*input)), output) + assert_array_equal(w.world_to_pixel_values(*w.pixel_to_world_values(*input)), output) + assert_array_equal(w.world_to_pixel(*w.pixel_to_world(*input)), output) + + +@pytest.mark.parametrize((("input", "output")), [(2, 2), + ((10, 200), (10, np.nan)), + (x, (np.nan, 2, 4, 13)) + ]) +def test_1d_freq_quantity(gwcs_1d_freq_quantity, input, output): + w = gwcs_1d_freq_quantity + #w.bounding_box = (-.5*u.pix, 21*u.pix) + w.bounding_box = (-.5, 21) + + # assert_array_equal(w.invert(w(input)), output) + # assert_array_equal(w.world_to_pixel_values(w.pixel_to_world_values(*input)), output) + # assert_array_equal(w.world_to_pixel(w.pixel_to_world(input)), output) + + +@pytest.mark.parametrize((("input", "output")), [((2, 4), (2, 4)), + ((100, 200), (np.nan, np.nan)), + ((x, x), (y, y)) + ]) +def test_2d_shift_scale_quantity(gwcs_2d_shift_scale_quantity, input, output): + w = gwcs_2d_shift_scale_quantity + w.bounding_box = ((-.5, 21), (4, 12)) + + assert_array_equal(w.invert(*w(*input)), output) + assert_array_equal(w.world_to_pixel_values(*w.pixel_to_world_values(*input)), output) + assert_array_equal(w.world_to_pixel(w.pixel_to_world(*input)), output) + + +@pytest.mark.parametrize((("input", "output")), [((2, 4, 5), (2, 4, 5)), + ((100, 200, 5), (np.nan, np.nan, np.nan)), + ((x, x, x), (y1, y1, y1)) + ]) +def test_3d_identity_units(gwcs_3d_identity_units, input, output): + w = gwcs_3d_identity_units + w.bounding_box = ((-.5, 21), (4, 12), (1, 21)) + + assert_array_equal(w.invert(*w(*input)), output) + assert_array_equal(w.world_to_pixel_values(*w.pixel_to_world_values(*input)), output) + assert_array_equal(w.world_to_pixel(w.pixel_to_world(*input)), output) + + +def test_4d_identity_units(gwcs_4d_identity_units, input, ooutput): + w = gwcs_4d_identity_units + w.bounding_box = ((-.5, 21), (4, 12), (1, 21), (5, 10)) \ No newline at end of file From 189e42ea4ef17c5033936e1f9dc6dab711947024 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Mon, 28 Oct 2024 17:29:50 -0400 Subject: [PATCH 11/19] update --- gwcs/api.py | 8 ++-- gwcs/tests/conftest.py | 82 +-------------------------------- gwcs/tests/test_bounding_box.py | 57 ++++++----------------- gwcs/wcs.py | 46 ++++++------------ 4 files changed, 33 insertions(+), 160 deletions(-) diff --git a/gwcs/api.py b/gwcs/api.py index 981a0ce9..5bf4fff0 100644 --- a/gwcs/api.py +++ b/gwcs/api.py @@ -7,6 +7,7 @@ from astropy.wcs.wcsapi import BaseHighLevelWCS, BaseLowLevelWCS from astropy.modeling import separable +from astropy.wcs.wcsapi.high_level_api import values_to_high_level_objects, high_level_objects_to_values import astropy.units as u from . import utils @@ -302,10 +303,6 @@ def _sanitize_pixel_inputs(self, *pixel_arrays): return pixels - def _sanitize_world_inputs(self, *world_arrays): - world_coord = [] - - def pixel_to_world(self, *pixel_arrays): """ Convert pixel values to world coordinates. @@ -326,8 +323,9 @@ def world_to_pixel(self, *world_objects): """ Convert world coordinates to pixel values. """ + #args = high_level_objects_to_values(*world_objects, low_level_wcs=self) + #result = self.invert(*args) result = self.invert(*world_objects, with_units=True) - if self.input_frame.naxes > 1: first_res = result[0] if not utils.isnumerical(first_res): diff --git a/gwcs/tests/conftest.py b/gwcs/tests/conftest.py index 62a0a2d5..d755e236 100644 --- a/gwcs/tests/conftest.py +++ b/gwcs/tests/conftest.py @@ -113,85 +113,5 @@ def sellmeier_zemax(): @pytest.fixture(scope="function") def gwcs_3d_galactic_spectral(): - return examples.gwcs_3d_galactic_spectral() - - -@pytest.fixture(scope="function") -def gwcs_1d_spectral(): - return examples.gwcs_1d_spectral() - - -@pytest.fixture(scope="function") -def gwcs_spec_cel_time_4d(): - return examples.gwcs_spec_cel_time_4d() - - -@pytest.fixture( - scope="function", - params=[ - (2, 1, 0), - (2, 0, 1), - pytest.param((1, 0, 2), marks=pytest.mark.skip(reason="Fails round-trip for -TAB axis 3")), - ] -) -def gwcs_cube_with_separable_spectral(request): - axes_order = request.param - return examples.gwcs_cube_with_separable_spectral(axes_order) - - -@pytest.fixture( - scope="function", - params=[ - (2, 0, 1), - (2, 1, 0), - pytest.param((0, 2, 1), marks=pytest.mark.skip(reason="Fails round-trip for -TAB axis 2")), - pytest.param((1, 0, 2), marks=pytest.mark.skip(reason="Fails round-trip for -TAB axis 3")), - ] -) -def gwcs_cube_with_separable_time(request): - axes_order = request.param - return examples.gwcs_cube_with_separable_time(axes_order) - - -@pytest.fixture(scope="function") -def gwcs_7d_complex_mapping(): - return examples.gwcs_7d_complex_mapping() - -@pytest.fixture -def spher_to_cart(): - return geometry.SphericalToCartesian() - - -@pytest.fixture -def cart_to_spher(): - return geometry.CartesianToSpherical() - - -@pytest.fixture -def gwcs_simple_imaging_no_units(): - 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"), - 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 = wcs.WCS(pipeline) - w.bounding_box = ((2, 100), (5, 500)) - return w + return examples.gwcs_3d_galactic_spectral() diff --git a/gwcs/tests/test_bounding_box.py b/gwcs/tests/test_bounding_box.py index 5fc13f79..335dcc4c 100644 --- a/gwcs/tests/test_bounding_box.py +++ b/gwcs/tests/test_bounding_box.py @@ -1,8 +1,6 @@ import numpy as np from numpy.testing import assert_array_equal, assert_allclose -from astropy import units as u - import pytest @@ -13,7 +11,7 @@ @pytest.mark.parametrize((("input", "output")), [((2, 4), (2, 4)), ((100, 200), (np.nan, np.nan)), - ((x, x), (y, y)) + ((x, x),(y, y)) ]) def test_2d_spatial(gwcs_2d_spatial_shift, input, output): w = gwcs_2d_spatial_shift @@ -52,7 +50,7 @@ def test_2d_spatial_coordinate_reordered(gwcs_2d_spatial_reordered, input, outpu @pytest.mark.parametrize((("input", "output")), [(2, 2), ((10, 200), (10, np.nan)), - (x, (np.nan, 2, 4, 13)) + (x, (np.nan, 2, 4, 13)) ]) def test_1d_freq(gwcs_1d_freq, input, output): w = gwcs_1d_freq @@ -76,46 +74,19 @@ def test_3d_spatial_wave(gwcs_3d_spatial_wave, input, output): assert_array_equal(w.world_to_pixel(*w.pixel_to_world(*input)), output) -@pytest.mark.parametrize((("input", "output")), [(2, 2), - ((10, 200), (10, np.nan)), - (x, (np.nan, 2, 4, 13)) - ]) -def test_1d_freq_quantity(gwcs_1d_freq_quantity, input, output): - w = gwcs_1d_freq_quantity - #w.bounding_box = (-.5*u.pix, 21*u.pix) - w.bounding_box = (-.5, 21) - - # assert_array_equal(w.invert(w(input)), output) - # assert_array_equal(w.world_to_pixel_values(w.pixel_to_world_values(*input)), output) - # assert_array_equal(w.world_to_pixel(w.pixel_to_world(input)), output) - - -@pytest.mark.parametrize((("input", "output")), [((2, 4), (2, 4)), - ((100, 200), (np.nan, np.nan)), - ((x, x), (y, y)) - ]) -def test_2d_shift_scale_quantity(gwcs_2d_shift_scale_quantity, input, output): - w = gwcs_2d_shift_scale_quantity - w.bounding_box = ((-.5, 21), (4, 12)) - - assert_array_equal(w.invert(*w(*input)), output) - assert_array_equal(w.world_to_pixel_values(*w.pixel_to_world_values(*input)), output) - assert_array_equal(w.world_to_pixel(w.pixel_to_world(*input)), output) - - -@pytest.mark.parametrize((("input", "output")), [((2, 4, 5), (2, 4, 5)), - ((100, 200, 5), (np.nan, np.nan, np.nan)), - ((x, x, x), (y1, y1, y1)) +@pytest.mark.parametrize((("input", "output")), [((1, 2, 3, 4), (1., 2., 3., 4.)), + ((100, 3, 3, 3), (np.nan, 3, 3, 3)), + ((x, x, x, x), [[np.nan, 2., 4., 13.], + [np.nan, 2., 4., 13.], + [np.nan, 2., 4., 13.], + [np.nan, 2., 4., np.nan]]) ]) -def test_3d_identity_units(gwcs_3d_identity_units, input, output): - w = gwcs_3d_identity_units - w.bounding_box = ((-.5, 21), (4, 12), (1, 21)) +def test_gwcs_spec_cel_time_4d(gwcs_spec_cel_time_4d, input, output): + w = gwcs_spec_cel_time_4d - assert_array_equal(w.invert(*w(*input)), output) - assert_array_equal(w.world_to_pixel_values(*w.pixel_to_world_values(*input)), output) - assert_array_equal(w.world_to_pixel(w.pixel_to_world(*input)), output) + assert_allclose(w.invert(*w(*input, with_bounding_box=False)), output, atol=1e-8) -def test_4d_identity_units(gwcs_4d_identity_units, input, ooutput): - w = gwcs_4d_identity_units - w.bounding_box = ((-.5, 21), (4, 12), (1, 21), (5, 10)) \ No newline at end of file +# @pytest.mark.parametrize((("input", "output")), [((2, 4, 5), (2, 4, 5))] +# def test_gwcs_1d_freq_quantity(gwcs_1d_freq_quantity, input, output): +# w = gwcs_1d_freq_quantity \ No newline at end of file diff --git a/gwcs/wcs.py b/gwcs/wcs.py index b58a400f..28102fc6 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -13,7 +13,11 @@ from astropy.modeling.models import (Const1D, Identity, Mapping, Polynomial2D, RotateCelestial2Native, Shift, Sky2Pix_TAN) +from astropy.modeling.parameters import _tofloat from astropy.wcs.utils import celestial_frame_to_wcs, proj_plane_pixel_scales +from astropy.wcs.wcsapi.high_level_api import high_level_objects_to_values + +from astropy import units as u from scipy import linalg, optimize from . import coordinate_frames as cf @@ -471,9 +475,7 @@ def invert(self, *args, **kwargs): btrans = self.backward_transform except NotImplementedError: btrans = None - print(f"args[0], {args[0]}") if not utils.isnumerical(args[0]): - print(f"args1, {args}") # convert astropy objects to numbers and arrays args = self.output_frame.coordinate_to_quantity(*args) if self.output_frame.naxes == 1: @@ -486,12 +488,10 @@ def invert(self, *args, **kwargs): with_bounding_box = kwargs.pop('with_bounding_box', True) fill_value = kwargs.pop('fill_value', np.nan) akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS} - print(f"args, {args}") if with_bounding_box and self.bounding_box is not None: result = self.outside_footprint(args) if btrans is not None: - #akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS} result = btrans(*args, **akwargs) else: result = self.numerical_inverse(*args, **kwargs, with_units=with_units) @@ -507,38 +507,27 @@ def invert(self, *args, **kwargs): return self.input_frame.coordinates(*result) else: return result - + def outside_footprint(self, world_arrays): - # for axis in world_arrays: - # if np.isscalar(axis): - # world_arrays = np.asarray(list(world_arrays), dtype = np.float64) - # world_arrays = [world_arrays] - #print('axis', axis, world_arrays) - #if np.isscalar(axis) or self.output_frame.naxes == 1: - if self.output_frame.naxes == 1: - #print('axis', axis) - #axis = float(axis) - world_arrays = [world_arrays] - #print('axis', axis) - print('world_arrays1', world_arrays) - #world_arrays = np.asarray(list(world_arrays))#, dtype = np.float64) - print('world_arrays2', world_arrays) + world_arrays = list(world_arrays) axes_types = set(self.output_frame.axes_type) + footprint = self.footprint() + world_arrays = [coo.to(unit) for coo, unit in zip(world_arrays, self.output_frame.unit) + if isinstance(coo, u.Quantity)] + world_arrays = [high_level_objects_to_values(coo, low_level_wcs=self) for + coo in world_arrays if not utils.isnumerical(coo)] + for axtyp in axes_types: - footprint = self.footprint(axis_type=axtyp) - ind = np.asarray((np.asarray(self.output_frame.axes_type) == axtyp)) - #print('ind', ind) - for idim, coord in enumerate(world_arrays[ind]): - #print('footprint', footprint) + for idim, coord in enumerate(world_arrays): + coord = _tofloat(coord) if np.asarray(ind).sum() > 1: axis_range = footprint[:, idim] else: axis_range = footprint range = [axis_range.min(), axis_range.max()] - #print('idim', idim, coord, range) outside = (coord < range[0]) | (coord > range[1]) if np.any(outside): if np.isscalar(coord): @@ -1473,11 +1462,6 @@ def footprint(self, bounding_box=None, center=False, axis_type="all"): """ def _order_clockwise(v): - # if self.input_frame.naxes == 1: - # bb = self.bounding_box.bounding_box() - # if isinstance(bb[0], u.Quantity): - # bb = [v.value for v in bb] * bb[0].unit - # return (bb,) return np.asarray([[v[0][0], v[1][0]], [v[0][0], v[1][1]], [v[0][1], v[1][1]], [v[0][1], v[1][0]]]).T @@ -1508,7 +1492,7 @@ def _order_clockwise(v): axis_type = axis_type.lower() if axis_type == 'spatial' and all_spatial: return result.T - + if axis_type != "all": axtyp_ind = np.array([t.lower() for t in self.output_frame.axes_type]) == axis_type if not axtyp_ind.any(): From df75fe2592b1f6f9d11ed420db8bc25705730743 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Mon, 28 Oct 2024 21:54:02 -0400 Subject: [PATCH 12/19] removed unused code --- gwcs/api.py | 1 - 1 file changed, 1 deletion(-) diff --git a/gwcs/api.py b/gwcs/api.py index 5bf4fff0..4ef35cb5 100644 --- a/gwcs/api.py +++ b/gwcs/api.py @@ -7,7 +7,6 @@ from astropy.wcs.wcsapi import BaseHighLevelWCS, BaseLowLevelWCS from astropy.modeling import separable -from astropy.wcs.wcsapi.high_level_api import values_to_high_level_objects, high_level_objects_to_values import astropy.units as u from . import utils From d847767840bac3778b99e579e47135a9ff5af887 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Wed, 30 Oct 2024 16:53:14 -0400 Subject: [PATCH 13/19] enforce 'F' order on bounding_box when calculating the footprint --- gwcs/wcs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 28102fc6..94a10bd5 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -1468,7 +1468,7 @@ def _order_clockwise(v): if bounding_box is None: if self.bounding_box is None: raise TypeError("Need a valid bounding_box to compute the footprint.") - bb = self.bounding_box.bounding_box() + bb = self.bounding_box.bounding_box(order='F') else: bb = bounding_box From 3d3b6c4a550c4087aff2e3d75ec9e89b425a8d25 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Wed, 30 Oct 2024 20:26:56 -0400 Subject: [PATCH 14/19] add newline --- gwcs/tests/test_bounding_box.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/gwcs/tests/test_bounding_box.py b/gwcs/tests/test_bounding_box.py index 335dcc4c..1cc46404 100644 --- a/gwcs/tests/test_bounding_box.py +++ b/gwcs/tests/test_bounding_box.py @@ -85,8 +85,3 @@ def test_gwcs_spec_cel_time_4d(gwcs_spec_cel_time_4d, input, output): w = gwcs_spec_cel_time_4d assert_allclose(w.invert(*w(*input, with_bounding_box=False)), output, atol=1e-8) - - -# @pytest.mark.parametrize((("input", "output")), [((2, 4, 5), (2, 4, 5))] -# def test_gwcs_1d_freq_quantity(gwcs_1d_freq_quantity, input, output): -# w = gwcs_1d_freq_quantity \ No newline at end of file From fb5f70d1035562765f1052d8147f657b47f9c825 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Tue, 10 Dec 2024 12:53:52 -0500 Subject: [PATCH 15/19] adding removed tests --- gwcs/tests/conftest.py | 51 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/gwcs/tests/conftest.py b/gwcs/tests/conftest.py index d755e236..bacbd91e 100644 --- a/gwcs/tests/conftest.py +++ b/gwcs/tests/conftest.py @@ -115,3 +115,54 @@ def sellmeier_zemax(): def gwcs_3d_galactic_spectral(): return examples.gwcs_3d_galactic_spectral() + +@pytest.fixture(scope="function") +def gwcs_1d_spectral(): + return examples.gwcs_1d_spectral() + + +@pytest.fixture(scope="function") +def gwcs_spec_cel_time_4d(): + return examples.gwcs_spec_cel_time_4d() + + +@pytest.fixture( + scope="function", + params=[ + (2, 1, 0), + (2, 0, 1), + pytest.param((1, 0, 2), marks=pytest.mark.skip(reason="Fails round-trip for -TAB axis 3")), + ] +) +def gwcs_cube_with_separable_spectral(request): + axes_order = request.param + return examples.gwcs_cube_with_separable_spectral(axes_order) + + +@pytest.fixture( + scope="function", + params=[ + (2, 0, 1), + (2, 1, 0), + pytest.param((0, 2, 1), marks=pytest.mark.skip(reason="Fails round-trip for -TAB axis 2")), + pytest.param((1, 0, 2), marks=pytest.mark.skip(reason="Fails round-trip for -TAB axis 3")), + ] +) +def gwcs_cube_with_separable_time(request): + axes_order = request.param + return examples.gwcs_cube_with_separable_time(axes_order) + + +@pytest.fixture(scope="function") +def gwcs_7d_complex_mapping(): + return examples.gwcs_7d_complex_mapping() + + +@pytest.fixture +def spher_to_cart(): + return geometry.SphericalToCartesian() + + +@pytest.fixture +def cart_to_spher(): + return geometry.CartesianToSpherical() From c455cab7e9c6e24161a09f5c60cbe8dd08198e63 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Tue, 10 Dec 2024 16:45:17 -0500 Subject: [PATCH 16/19] fix something --- gwcs/tests/test_wcs.py | 6 +++--- gwcs/wcs.py | 30 ++++-------------------------- 2 files changed, 7 insertions(+), 29 deletions(-) diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 1baad46f..5b2a47fb 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -1161,10 +1161,10 @@ def test_in_image(): None)]) w2.bounding_box = [(1, 100), (2, 20)] - assert np.isscalar(w2.in_image(2, 6)) + assert np.isscalar(w2.in_image(2, 6))#[0]) assert not np.isscalar(w2.in_image([2], [6])) - assert w2.in_image(4, 6) - assert not w2.in_image(5, 0) + assert (w2.in_image(4, 6))#.all() + assert not (w2.in_image(5, 0))#, [True, False]) assert np.array_equal( w2.in_image( [[9, 10, 11, 15], [8, 9, 67, 98], [2, 2, np.nan, 102]], diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 94a10bd5..97f2da14 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -400,31 +400,12 @@ def in_image(self, *args, **kwargs): and `False` if input is outside the footprint. """ - coords = self.invert(*args, with_bounding_box=False, **kwargs) + coords = self.invert(*args, **kwargs) result = np.isfinite(coords) if self.input_frame.naxes > 1: result = np.all(result, axis=0) - if self.bounding_box is None or not np.any(result): - return result - - if self.input_frame.naxes == 1: - x1, x2 = self.bounding_box.bounding_box() - - if len(np.shape(args[0])) > 0: - result[result] = (coords[result] >= x1) & (coords[result] <= x2) - elif result: - result = (coords >= x1) and (coords <= x2) - - else: - if len(np.shape(args[0])) > 0: - for c, (x1, x2) in zip(coords, self.bounding_box): - result[result] = (c[result] >= x1) & (c[result] <= x2) - - elif result: - result = all([(c >= x1) and (c <= x2) for c, (x1, x2) in zip(coords, self.bounding_box)]) - return result def invert(self, *args, **kwargs): @@ -489,7 +470,7 @@ def invert(self, *args, **kwargs): fill_value = kwargs.pop('fill_value', np.nan) akwargs = {k: v for k, v in kwargs.items() if k not in _ITER_INV_KWARGS} if with_bounding_box and self.bounding_box is not None: - result = self.outside_footprint(args) + args = self.outside_footprint(args) if btrans is not None: result = btrans(*args, **akwargs) @@ -513,11 +494,8 @@ def outside_footprint(self, world_arrays): axes_types = set(self.output_frame.axes_type) footprint = self.footprint() - world_arrays = [coo.to(unit) for coo, unit in zip(world_arrays, self.output_frame.unit) - if isinstance(coo, u.Quantity)] - world_arrays = [high_level_objects_to_values(coo, low_level_wcs=self) for - coo in world_arrays if not utils.isnumerical(coo)] - + if not utils.isnumerical(world_arrays[0]): + world_arrays = high_level_objects_to_values(*world_arrays, low_level_wcs=self) for axtyp in axes_types: ind = np.asarray((np.asarray(self.output_frame.axes_type) == axtyp)) From ef8f619c96d99e3091228c463265ab979bd6ab91 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Tue, 10 Dec 2024 17:02:29 -0500 Subject: [PATCH 17/19] fix asdf call --- gwcs/wcs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 97f2da14..db0b1a80 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -718,7 +718,7 @@ def numerical_inverse(self, *args, tolerance=1e-5, maxiter=30, adaptive=True, >>> import numpy as np >>> filename = get_pkg_data_filename('data/nircamwcs.asdf', package='gwcs.tests') - >>> with asdf.open(filename, copy_arrays=True, lazy_load=False, ignore_missing_extensions=True) as af: + >>> with asdf.open(filename, lazy_load=False, ignore_missing_extensions=True) as af: ... w = af.tree['wcs'] >>> ra, dec = w([1,2,3], [1,1,1]) From 96de876179f04e4f114d3a1f8593a141e20cd0dc Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Thu, 12 Dec 2024 07:26:37 -0500 Subject: [PATCH 18/19] fix specutils --- gwcs/wcs.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index db0b1a80..a52dc06e 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -15,7 +15,7 @@ Sky2Pix_TAN) from astropy.modeling.parameters import _tofloat from astropy.wcs.utils import celestial_frame_to_wcs, proj_plane_pixel_scales -from astropy.wcs.wcsapi.high_level_api import high_level_objects_to_values +from astropy.wcs.wcsapi.high_level_api import high_level_objects_to_values, values_to_high_level_objects from astropy import units as u from scipy import linalg, optimize @@ -494,7 +494,9 @@ def outside_footprint(self, world_arrays): axes_types = set(self.output_frame.axes_type) footprint = self.footprint() + not_numerical = False if not utils.isnumerical(world_arrays[0]): + not_numerical = True world_arrays = high_level_objects_to_values(*world_arrays, low_level_wcs=self) for axtyp in axes_types: ind = np.asarray((np.asarray(self.output_frame.axes_type) == axtyp)) @@ -513,7 +515,8 @@ def outside_footprint(self, world_arrays): else: coord[outside] = np.nan world_arrays[idim] = coord - + if not_numerical: + world_arrays = values_to_high_level_objects(*world_arrays, low_level_wcs=self) return world_arrays From 08bafd07078c890296437b989169820c389948d6 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Fri, 13 Dec 2024 15:59:24 -0500 Subject: [PATCH 19/19] remove commented out code --- gwcs/api.py | 2 -- gwcs/tests/conftest.py | 25 ------------------------- gwcs/tests/test_coordinate_systems.py | 3 --- gwcs/tests/test_wcs.py | 6 +++--- 4 files changed, 3 insertions(+), 33 deletions(-) diff --git a/gwcs/api.py b/gwcs/api.py index 4ef35cb5..32f6f9c3 100644 --- a/gwcs/api.py +++ b/gwcs/api.py @@ -322,8 +322,6 @@ def world_to_pixel(self, *world_objects): """ Convert world coordinates to pixel values. """ - #args = high_level_objects_to_values(*world_objects, low_level_wcs=self) - #result = self.invert(*args) result = self.invert(*world_objects, with_units=True) if self.input_frame.naxes > 1: first_res = result[0] diff --git a/gwcs/tests/conftest.py b/gwcs/tests/conftest.py index bacbd91e..92979383 100644 --- a/gwcs/tests/conftest.py +++ b/gwcs/tests/conftest.py @@ -5,31 +5,6 @@ from .. import examples from .. import geometry -import numpy as np - -import astropy.units as u -from astropy.time import Time -from astropy import coordinates as coord -from astropy.modeling import models - -from gwcs import coordinate_frames as cf -from gwcs import spectroscopy as sp -from gwcs import wcs -from gwcs import geometry - -# frames -detector_1d = cf.CoordinateFrame(name='detector', axes_order=(0,), naxes=1, axes_type="detector") -detector_2d = cf.Frame2D(name='detector', axes_order=(0, 1)) -icrs_sky_frame = cf.CelestialFrame(reference_frame=coord.ICRS(), - axes_order=(0, 1)) - -freq_frame = cf.SpectralFrame(name='freq', unit=u.Hz, axes_order=(0, )) -wave_frame = cf.SpectralFrame(name='wave', unit=u.m, axes_order=(2, ), - axes_names=('lambda', )) - -# transforms -model_2d_shift = models.Shift(1) & models.Shift(2) -model_1d_scale = models.Scale(2) @pytest.fixture diff --git a/gwcs/tests/test_coordinate_systems.py b/gwcs/tests/test_coordinate_systems.py index f1e34cb9..ac2c7303 100644 --- a/gwcs/tests/test_coordinate_systems.py +++ b/gwcs/tests/test_coordinate_systems.py @@ -190,7 +190,6 @@ def test_temporal_relative(): assert a[1] == Time("2018-01-01T00:00:00") + 20 * u.s -#@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher") def test_temporal_absolute(): t = cf.TemporalFrame(reference_frame=Time([], format='isot')) assert t.coordinates("2018-01-01T00:00:00") == Time("2018-01-01T00:00:00") @@ -240,7 +239,6 @@ def test_coordinate_to_quantity_spectral(inp): (Time("2011-01-01T00:00:10"),), (10 * u.s,) ]) -#@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher.") def test_coordinate_to_quantity_temporal(inp): temp = cf.TemporalFrame(reference_frame=Time("2011-01-01T00:00:00"), unit=u.s) @@ -325,7 +323,6 @@ def test_coordinate_to_quantity_frame_2d(): assert_quantity_allclose(output, exp) -#@pytest.mark.skipif(astropy_version<"4", reason="Requires astropy 4.0 or higher.") def test_coordinate_to_quantity_error(): frame = cf.Frame2D(unit=(u.one, u.arcsec)) with pytest.raises(ValueError): diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 5b2a47fb..23fc82c1 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -1161,10 +1161,10 @@ def test_in_image(): None)]) w2.bounding_box = [(1, 100), (2, 20)] - assert np.isscalar(w2.in_image(2, 6))#[0]) + assert np.isscalar(w2.in_image(2, 6)) assert not np.isscalar(w2.in_image([2], [6])) - assert (w2.in_image(4, 6))#.all() - assert not (w2.in_image(5, 0))#, [True, False]) + assert (w2.in_image(4, 6)) + assert not (w2.in_image(5, 0)) assert np.array_equal( w2.in_image( [[9, 10, 11, 15], [8, 9, 67, 98], [2, 2, np.nan, 102]],