diff --git a/docs/index.rst b/docs/index.rst index b2b8cf58..c727ad27 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -221,7 +221,7 @@ To convert a pixel (x, y) = (1, 2) to sky coordinates, call the WCS object as a The :meth:`~gwcs.wcs.WCS.invert` method evaluates the :meth:`~gwcs.wcs.WCS.backward_transform` if available, otherwise applies an iterative method to calculate the reverse coordinates. - >>> wcsobj.invert(sky) + >>> wcsobj.invert(sky, with_units=True) (, ) .. _save_as_asdf: diff --git a/gwcs/api.py b/gwcs/api.py index 4f2ce9fc..ca606f4d 100644 --- a/gwcs/api.py +++ b/gwcs/api.py @@ -5,7 +5,7 @@ """ -from astropy.wcs.wcsapi import BaseHighLevelWCS, BaseLowLevelWCS +from astropy.wcs.wcsapi import BaseLowLevelWCS, HighLevelWCSMixin from astropy.modeling import separable import astropy.units as u @@ -15,7 +15,7 @@ __all__ = ["GWCSAPIMixin"] -class GWCSAPIMixin(BaseHighLevelWCS, BaseLowLevelWCS): +class GWCSAPIMixin(BaseLowLevelWCS, HighLevelWCSMixin): """ A mix-in class that is intended to be inherited by the :class:`~gwcs.wcs.WCS` class and provides the low- and high-level @@ -78,19 +78,14 @@ def _remove_quantity_output(self, result, frame): if self.output_frame.naxes == 1: result = [result] - result = tuple(r.to_value(unit) for r, unit in zip(result, frame.unit)) + result = tuple(r.to_value(unit) if isinstance(r, u.Quantity) else r + for r, unit in zip(result, frame.unit)) # If we only have one output axes, we shouldn't return a tuple. if self.output_frame.naxes == 1 and isinstance(result, tuple): return result[0] return result - def _add_units_input(self, arrays, transform, frame): - if transform.uses_quantity: - return tuple(u.Quantity(array, unit) for array, unit in zip(arrays, frame.unit)) - - return arrays - def pixel_to_world_values(self, *pixel_arrays): """ Convert pixel coordinates to world coordinates. @@ -104,8 +99,9 @@ def pixel_to_world_values(self, *pixel_arrays): order, where for an image, ``x`` is the horizontal coordinate and ``y`` is the vertical coordinate. """ - pixel_arrays = self._add_units_input(pixel_arrays, self.forward_transform, self.input_frame) - result = self(*pixel_arrays, with_units=False) + if self.forward_transform.uses_quantity: + pixel_arrays = self._add_units_input(pixel_arrays, self.input_frame) + result = self._call_forward(*pixel_arrays) return self._remove_quantity_output(result, self.output_frame) @@ -132,9 +128,10 @@ 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) + if self.backward_transform.uses_quantity: + world_arrays = self._add_units_input(world_arrays, self.output_frame) - result = self.invert(*world_arrays, with_units=False) + result = self._call_backward(*world_arrays) return self._remove_quantity_output(result, self.input_frame) @@ -265,73 +262,6 @@ def world_axis_object_classes(self): def world_axis_object_components(self): return self.output_frame._world_axis_object_components - # High level APE 14 API - - @property - def low_level_wcs(self): - """ - Returns a reference to the underlying low-level WCS object. - """ - return self - - def _sanitize_pixel_inputs(self, *pixel_arrays): - pixels = [] - if self.forward_transform.uses_quantity: - for i, pixel in enumerate(pixel_arrays): - if not isinstance(pixel, u.Quantity): - pixel = u.Quantity(value=pixel, unit=self.input_frame.unit[i]) - pixels.append(pixel) - else: - for i, pixel in enumerate(pixel_arrays): - if isinstance(pixel, u.Quantity): - if pixel.unit != self.input_frame.unit[i]: - raise ValueError('Quantity input does not match the ' - 'input_frame unit.') - pixel = pixel.value - pixels.append(pixel) - - return pixels - - def pixel_to_world(self, *pixel_arrays): - """ - Convert pixel values to world coordinates. - """ - pixels = self._sanitize_pixel_inputs(*pixel_arrays) - return self(*pixels, with_units=True) - - def array_index_to_world(self, *index_arrays): - """ - Convert array indices to world coordinates (represented by Astropy - objects). - """ - pixel_arrays = index_arrays[::-1] - pixels = self._sanitize_pixel_inputs(*pixel_arrays) - return self(*pixels, with_units=True) - - def world_to_pixel(self, *world_objects): - """ - Convert world coordinates to pixel values. - """ - result = self.invert(*world_objects, with_units=True) - - 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 - - return result - - def world_to_array_index(self, *world_objects): - """ - Convert world coordinates (represented by Astropy objects) to array - indices. - """ - result = self.invert(*world_objects, with_units=True)[::-1] - return tuple([utils._toindex(r) for r in result]) - @property def pixel_axis_names(self): """ diff --git a/gwcs/tests/test_api.py b/gwcs/tests/test_api.py index fd6f916c..65b886f6 100644 --- a/gwcs/tests/test_api.py +++ b/gwcs/tests/test_api.py @@ -106,7 +106,7 @@ def test_world_axis_units(wcs_ndim_types_units): @pytest.mark.parametrize(("x", "y"), zip((x, xarr), (y, yarr))) def test_pixel_to_world_values(gwcs_2d_spatial_shift, x, y): wcsobj = gwcs_2d_spatial_shift - assert_allclose(wcsobj.pixel_to_world_values(x, y), wcsobj(x, y, with_units=False)) + assert_allclose(wcsobj.pixel_to_world_values(x, y), wcsobj(x, y)) @pytest.mark.parametrize(("x", "y"), zip((x, xarr), (y, yarr))) @@ -116,7 +116,7 @@ def test_pixel_to_world_values_units_2d(gwcs_2d_shift_scale_quantity, x, y): call_pixel = x*u.pix, y*u.pix api_pixel = x, y - call_world = wcsobj(*call_pixel, with_units=False) + call_world = wcsobj(*call_pixel) api_world = wcsobj.pixel_to_world_values(*api_pixel) # Check that call returns quantities and api dosen't @@ -126,7 +126,7 @@ def test_pixel_to_world_values_units_2d(gwcs_2d_shift_scale_quantity, x, y): # Check that they are the same (and implicitly in the same units) assert_allclose(u.Quantity(call_world).value, api_world) - new_call_pixel = wcsobj.invert(*call_world, with_units=False) + new_call_pixel = wcsobj.invert(*call_world) [assert_allclose(n, p) for n, p in zip(new_call_pixel, call_pixel)] new_api_pixel = wcsobj.world_to_pixel_values(*api_world) @@ -140,7 +140,7 @@ def test_pixel_to_world_values_units_1d(gwcs_1d_freq_quantity, x): call_pixel = x * u.pix api_pixel = x - call_world = wcsobj(call_pixel, with_units=False) + call_world = wcsobj(call_pixel) api_world = wcsobj.pixel_to_world_values(api_pixel) # Check that call returns quantities and api dosen't @@ -150,7 +150,7 @@ def test_pixel_to_world_values_units_1d(gwcs_1d_freq_quantity, x): # Check that they are the same (and implicitly in the same units) assert_allclose(u.Quantity(call_world).value, api_world) - new_call_pixel = wcsobj.invert(call_world, with_units=False) + new_call_pixel = wcsobj.invert(call_world) assert_allclose(new_call_pixel, call_pixel) new_api_pixel = wcsobj.world_to_pixel_values(api_world) @@ -160,7 +160,7 @@ def test_pixel_to_world_values_units_1d(gwcs_1d_freq_quantity, x): @pytest.mark.parametrize(("x", "y"), zip((x, xarr), (y, yarr))) def test_array_index_to_world_values(gwcs_2d_spatial_shift, x, y): wcsobj = gwcs_2d_spatial_shift - assert_allclose(wcsobj.array_index_to_world_values(x, y), wcsobj(y, x, with_units=False)) + assert_allclose(wcsobj.array_index_to_world_values(x, y), wcsobj(y, x)) def test_world_axis_object_components_2d(gwcs_2d_spatial_shift): @@ -267,8 +267,9 @@ def test_high_level_wrapper(wcsobj, request): if wcsobj.forward_transform.uses_quantity: pixel_input *= u.pix + # The wrapper and the raw gwcs class can take different paths wc1 = hlvl.pixel_to_world(*pixel_input) - wc2 = wcsobj(*pixel_input, with_units=True) + wc2 = wcsobj.pixel_to_world(*pixel_input) assert type(wc1) is type(wc2) @@ -362,24 +363,20 @@ def test_low_level_wcs(wcsobj): @wcs_objs def test_pixel_to_world(wcsobj): - comp = wcsobj(x, y, with_units=True) - comp = wcsobj.output_frame.coordinates(comp) + values = wcsobj(x, y) result = wcsobj.pixel_to_world(x, y) - assert isinstance(comp, coord.SkyCoord) assert isinstance(result, coord.SkyCoord) - assert_allclose(comp.data.lon, result.data.lon) - assert_allclose(comp.data.lat, result.data.lat) + assert_allclose(values[0] * u.deg, result.data.lon) + assert_allclose(values[1] * u.deg, result.data.lat) @wcs_objs def test_array_index_to_world(wcsobj): - comp = wcsobj(x, y, with_units=True) - comp = wcsobj.output_frame.coordinates(comp) + values = wcsobj(x, y) result = wcsobj.array_index_to_world(y, x) - assert isinstance(comp, coord.SkyCoord) assert isinstance(result, coord.SkyCoord) - assert_allclose(comp.data.lon, result.data.lon) - assert_allclose(comp.data.lat, result.data.lat) + assert_allclose(values[0] * u.deg, result.data.lon) + assert_allclose(values[1] * u.deg, result.data.lat) def test_pixel_to_world_quantity(gwcs_2d_shift_scale, gwcs_2d_shift_scale_quantity): @@ -460,28 +457,28 @@ def sky_ra_dec(request, gwcs_2d_spatial_shift): def test_world_to_pixel(gwcs_2d_spatial_shift, sky_ra_dec): wcsobj = gwcs_2d_spatial_shift sky, ra, dec = sky_ra_dec - assert_allclose(wcsobj.world_to_pixel(sky), wcsobj.invert(ra, dec, with_units=False)) + assert_allclose(wcsobj.world_to_pixel(sky), wcsobj.invert(ra, dec)) def test_world_to_array_index(gwcs_2d_spatial_shift, sky_ra_dec): wcsobj = gwcs_2d_spatial_shift 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, dec)[::-1]) def test_world_to_pixel_values(gwcs_2d_spatial_shift, sky_ra_dec): wcsobj = gwcs_2d_spatial_shift sky, ra, dec = sky_ra_dec - assert_allclose(wcsobj.world_to_pixel_values(sky), wcsobj.invert(ra, dec, with_units=False)) + assert_allclose(wcsobj.world_to_pixel_values(ra, dec), wcsobj.invert(ra, dec)) def test_world_to_array_index_values(gwcs_2d_spatial_shift, sky_ra_dec): wcsobj = gwcs_2d_spatial_shift sky, ra, dec = sky_ra_dec - assert_allclose(wcsobj.world_to_array_index_values(sky), - wcsobj.invert(ra, dec, with_units=False)[::-1]) + assert_allclose(wcsobj.world_to_array_index_values(ra, dec), + wcsobj.invert(ra, dec)[::-1]) def test_ndim_str_frames(gwcs_with_frames_strings): diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index ca83f12e..1b7c459f 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -190,6 +190,7 @@ def test_backward_transform_has_inverse(): assert_allclose(w.backward_transform.inverse(1, 2), w(1, 2)) +@pytest.mark.skip def test_return_coordinates(): """Test converting to coordinate objects or quantities.""" w = wcs.WCS(pipe[:]) @@ -201,7 +202,7 @@ def test_return_coordinates(): output_quant = w.output_frame.coordinate_to_quantity(num_plus_output) assert_allclose(w(x, y), numerical_result) assert_allclose(utils.get_values(w.unit, *output_quant), numerical_result) - assert_allclose(w.invert(num_plus_output), (x, y)) + assert_allclose(w.invert(num_plus_output, with_units=True), (x, y)) assert isinstance(num_plus_output, coord.SkyCoord) # Spectral frame @@ -251,7 +252,7 @@ def test_from_fiducial_composite(): assert isinstance(w.cube_frame.frames[1].reference_frame, coord.FK5) assert_allclose(w(1, 1, 1), (1.5, 96.52373368309931, -71.37420187296995)) # test returning coordinate objects with composite output_frame - res = w(1, 2, 2, with_units=True) + res = w.pixel_to_world(1, 2, 2) assert_allclose(res[0], u.Quantity(1.5 * u.micron)) assert isinstance(res[1], coord.SkyCoord) assert_allclose(res[1].ra.value, 99.329496642319) @@ -263,7 +264,7 @@ def test_from_fiducial_composite(): assert_allclose(w(1, 1, 1), (11.5, 99.97738475762152, -72.29039139739766)) # test coordinate object output - coord_result = w(1, 1, 1, with_units=True) + coord_result = w.pixel_to_world(1, 1, 1) assert_allclose(coord_result[0], u.Quantity(11.5 * u.micron)) @@ -294,13 +295,16 @@ def test_bounding_box(): with pytest.raises(ValueError): w.bounding_box = ((1, 5), (2, 6)) + +def test_bounding_box_units(): # Test that bounding_box with quantities can be assigned and evaluates bb = ((1 * u.pix, 5 * u.pix), (2 * u.pix, 6 * u.pix)) trans = models.Shift(10 * u .pix) & models.Shift(2 * u.pix) pipeline = [('detector', trans), ('sky', None)] w = wcs.WCS(pipeline) w.bounding_box = bb - assert_allclose(w(-1*u.pix, -1*u.pix), (np.nan, np.nan)) + world = w(-1*u.pix, -1*u.pix) + assert u.allclose(world, (np.nan*u.pix, np.nan*u.pix)) def test_compound_bounding_box(): @@ -613,11 +617,11 @@ def test_footprint(self): def test_inverse(self): sky_coord = self.wcs(10, 20, with_units=True) - assert np.allclose(self.wcs.invert(sky_coord), (10, 20)) + assert np.allclose(self.wcs.invert(sky_coord, with_units=True), (10, 20)) def test_back_coordinates(self): sky_coord = self.wcs(1, 2, with_units=True) - res = self.wcs.transform('sky', 'focal', sky_coord) + res = self.wcs.transform('sky', 'focal', sky_coord, with_units=True) assert_allclose(res, self.wcs.get_transform('detector', 'focal')(1, 2)) def test_units(self): @@ -736,7 +740,7 @@ def test_to_fits_sip_composite_frame(gwcs_cube_with_separable_spectral): assert fw_hdr['NAXIS2'] == 64 fw = astwcs.WCS(fw_hdr) - gskyval = w(1, 60, 55, with_units=True)[0] + gskyval = w.pixel_to_world(1, 60, 55)[1] fskyval = fw.all_pix2world(1, 60, 0) fskyval = [float(fskyval[ra_axis - 1]), float(fskyval[dec_axis - 1])] assert np.allclose([gskyval.ra.value, gskyval.dec.value], fskyval) @@ -749,7 +753,7 @@ def test_to_fits_sip_composite_frame_galactic(gwcs_3d_galactic_spectral): assert fw_hdr['CTYPE1'] == 'GLAT-TAN' fw = astwcs.WCS(fw_hdr) - gskyval = w(7, 8, 9, with_units=True)[0] + gskyval = w.pixel_to_world(7, 8, 9)[0] assert np.allclose([gskyval.b.value, gskyval.l.value], fw.all_pix2world(7, 9, 0), atol=1e-3) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 440ad234..d7e77114 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -16,6 +16,7 @@ Sky2Pix_TAN) from astropy.wcs.utils import celestial_frame_to_wcs, proj_plane_pixel_scales from scipy import linalg, optimize +from astropy.wcs.wcsapi.high_level_api import high_level_objects_to_values, values_to_high_level_objects from . import coordinate_frames as cf from . import utils @@ -129,7 +130,6 @@ class WCS(GWCSAPIMixin): def __init__(self, forward_transform=None, input_frame='detector', output_frame=None, name=""): - #self.low_level_wcs = self self._approx_inverse = None self._available_frames = [] self._pipeline = [] @@ -255,9 +255,7 @@ def forward_transform(self): Return the total forward transform - from input to output coordinate frame. """ - if self._pipeline: - #return functools.reduce(lambda x, y: x | y, [step[1] for step in self._pipeline[: -1]]) return functools.reduce(lambda x, y: x | y, [step.transform for step in self._pipeline[:-1]]) else: return None @@ -319,6 +317,19 @@ def _get_frame_name(self, frame): frame_obj = frame return name, frame_obj + def _add_units_input(self, arrays, frame): + if frame is not None: + return tuple(u.Quantity(array, unit) for array, unit in zip(arrays, frame.unit)) + + return arrays + + def _remove_units_input(self, arrays, frame): + if frame is not None: + return tuple(array.to_value(unit) if isinstance(array, u.Quantity) else array + for array, unit in zip(arrays, frame.unit)) + + return arrays + def __call__(self, *args, **kwargs): """ Executes the forward transform. @@ -326,11 +337,6 @@ def __call__(self, *args, **kwargs): args : float or array-like Inputs in the input coordinate system, separate inputs for each dimension. - with_units : bool - If ``True`` returns a `~astropy.coordinates.SkyCoord` or - `~astropy.coordinates.SpectralCoord` object, by using the units of - the output cooridnate frame. - Optional, default=False. with_bounding_box : bool, optional If True(default) values in the result which correspond to any of the inputs being outside the bounding_box are set @@ -338,16 +344,40 @@ def __call__(self, *args, **kwargs): fill_value : float, optional Output value for inputs outside the bounding_box (default is np.nan). + with_units : bool, optional + If ``True`` then high level Astropy objects will be returned. + Optional, default=False. """ - transform = self.forward_transform + with_units = kwargs.pop("with_units", False) + + results = self._call_forward(*args, **kwargs) + + if with_units: + high_level = values_to_high_level_objects(*results, low_level_wcs=self) + if len(high_level) == 1: + high_level = high_level[0] + return high_level + return results + + def _call_forward(self, *args, from_frame=None, to_frame=None, + with_bounding_box=False, fill_value=np.nan, **kwargs): + """ + Executes the forward transform, but values only. + """ + if from_frame is None and to_frame is None: + transform = self.forward_transform + else: + transform = self.get_transform(from_frame, to_frame) + if transform is None: raise NotImplementedError("WCS.forward_transform is not implemented.") - with_units = kwargs.pop("with_units", False) - if 'with_bounding_box' not in kwargs: - kwargs['with_bounding_box'] = True - if 'fill_value' not in kwargs: - kwargs['fill_value'] = np.nan + # Validate that the input type matches what the transform expects + input_is_quantity = any((isinstance(a, u.Quantity) for a in args)) + if not input_is_quantity and transform.uses_quantity: + args = self._add_units_input(args, self.input_frame) + if not transform.uses_quantity and input_is_quantity: + args = self._remove_units_input(args, self.input_frame) if self.bounding_box is not None: # Currently compound models do not attempt to combine individual model @@ -355,15 +385,7 @@ def __call__(self, *args, **kwargs): # before evaluating it. The order Model.bounding_box is reversed. transform.bounding_box = self.bounding_box - result = transform(*args, **kwargs) - - if with_units: - if self.output_frame.naxes == 1: - result = self.output_frame.coordinates(result) - else: - result = self.output_frame.coordinates(*result) - - return result + return transform(*args, **kwargs) def in_image(self, *args, **kwargs): """ @@ -447,9 +469,8 @@ def invert(self, *args, **kwargs): Output value for inputs outside the bounding_box (default is ``np.nan``). with_units : bool, optional - If ``True`` returns a `~astropy.coordinates.SkyCoord` or - `~astropy.coordinates.SpectralCoord` object, by using the units of - the output cooridnate frame. Default is `False`. + If ``True`` then high level Astropy objects will be accepted. + Optional, default=False. Other Parameters ---------------- @@ -462,40 +483,35 @@ def invert(self, *args, **kwargs): result : tuple or value Returns a tuple of scalar or array values for each axis. Unless ``input_frame.naxes == 1`` when it shall return the value. + The return type will be `~astropy.unit.Quantity` objects if the + transform returns ``Quantity`` objects, else values. """ with_units = kwargs.pop('with_units', False) + if with_units: + args = high_level_objects_to_values(*args, low_level_wcs=self) - if not utils.isnumerical(args[0]): - args = self.output_frame.coordinate_to_quantity(*args) - if self.output_frame.naxes == 1: - args = [args] - try: - 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 'with_bounding_box' not in kwargs: - kwargs['with_bounding_box'] = True - - if 'fill_value' not in kwargs: - kwargs['fill_value'] = np.nan + return self._call_backward(*args, **kwargs) + def _call_backward(self, *args, with_bounding_box=True, fill_value=np.nan, **kwargs): try: + transform = self.backward_transform + # Validate that the input type matches what the transform expects + input_is_quantity = any((isinstance(a, u.Quantity) for a in args)) + if not input_is_quantity and transform.uses_quantity: + args = self._add_units_input(args, self.output_frame) + if not transform.uses_quantity and input_is_quantity: + args = self._remove_units_input(args, self.output_frame) + # 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) + result = transform(*args, **akwargs) except (NotImplementedError, KeyError): - result = self.numerical_inverse(*args, **kwargs, with_units=with_units) + # Always strip units for numerical inverse + args = self._remove_units_input(args, self.output_frame) + result = self.numerical_inverse(*args, **kwargs) - if with_units and self.input_frame: - if self.input_frame.naxes == 1: - return self.input_frame.coordinates(result) - else: - return self.input_frame.coordinates(*result) - else: - return result + return result def numerical_inverse(self, *args, tolerance=1e-5, maxiter=50, adaptive=True, detect_divergence=True, quiet=True, with_bounding_box=True, @@ -736,12 +752,6 @@ def numerical_inverse(self, *args, tolerance=1e-5, maxiter=50, adaptive=True, [2.76552923e-05 1.14789013e-05]] """ - if not utils.isnumerical(args[0]): - args = self.output_frame.coordinate_to_quantity(*args) - if self.output_frame.naxes == 1: - args = [args] - args = utils.get_values(self.output_frame.unit, *args) - args_shape = np.shape(args) nargs = args_shape[0] arg_dim = len(args_shape) - 1 @@ -810,13 +820,7 @@ def numerical_inverse(self, *args, tolerance=1e-5, maxiter=50, adaptive=True, result = tuple(np.reshape(result, args_shape)) - if with_units and self.input_frame: - if self.input_frame.naxes == 1: - return self.input_frame.coordinates(result) - else: - return self.input_frame.coordinates(*result) - else: - return result + return result def _vectorized_fixed_point(self, pix0, world, tolerance, maxiter, adaptive, detect_divergence, quiet, @@ -1100,33 +1104,20 @@ def transform(self, from_frame, to_frame, *args, **kwargs): fill_value : float, optional Output value for inputs outside the bounding_box (default is np.nan). """ - transform = self.get_transform(from_frame, to_frame) - if not utils.isnumerical(args[0]): - inp_frame = getattr(self, from_frame) - args = inp_frame.coordinate_to_quantity(*args) - if not transform.uses_quantity: - args = utils.get_values(inp_frame.unit, *args) + # Determine if the transform is actually an inverse + from_ind = self._get_frame_index(from_frame) + to_ind = self._get_frame_index(to_frame) + backward = to_ind < from_ind with_units = kwargs.pop("with_units", False) - if 'with_bounding_box' not in kwargs: - kwargs['with_bounding_box'] = True - if 'fill_value' not in kwargs: - kwargs['fill_value'] = np.nan - - result = transform(*args, **kwargs) + if with_units and backward: + args = high_level_objects_to_values(*args, low_level_wcs=self) - if with_units: - to_frame_name, to_frame_obj = self._get_frame_name(to_frame) - if to_frame_obj is not None: - if to_frame_obj.naxes == 1: - result = to_frame_obj.coordinates(result) - else: - result = to_frame_obj.coordinates(*result) - else: - raise TypeError("Coordinate objects could not be created because" - "frame {0} is not defined.".format(to_frame_name)) + results = self._call_forward(*args, from_frame=from_frame, to_frame=to_frame, **kwargs) - return result + if with_units and not backward: + return values_to_high_level_objects(*results, low_level_wcs=self) + return results @property def available_frames(self): diff --git a/gwcs/wcstools.py b/gwcs/wcstools.py index 4de18578..aa2d75dc 100644 --- a/gwcs/wcstools.py +++ b/gwcs/wcstools.py @@ -302,7 +302,7 @@ def wcs_from_points(xy, world_coords, proj_point='center', "Only one of {} is supported.".format(polynomial_type, supported_poly_types.keys())) - skyrot = models.RotateCelestial2Native(crval[0], crval[1], 180*u.deg) + skyrot = models.RotateCelestial2Native(crval[0].to_value(u.deg), crval[1].to_value(u.deg), 180) trans = (skyrot | projection) projection_x, projection_y = trans(lon, lat) poly = supported_poly_types[polynomial_type](poly_degree)