From dbaadfcfb5a89656976d8f8929e8e581889a3db1 Mon Sep 17 00:00:00 2001 From: Stuart Mumford Date: Thu, 12 Oct 2023 17:04:44 -0400 Subject: [PATCH] Add a test for axes ordering with CelestialFrame fixes #269 --- gwcs/coordinate_frames.py | 17 ++++++++++++-- gwcs/tests/test_wcs.py | 49 ++++++++++++++++++++++++++++++++++++--- 2 files changed, 61 insertions(+), 5 deletions(-) diff --git a/gwcs/coordinate_frames.py b/gwcs/coordinate_frames.py index 6321b33d..961c591d 100644 --- a/gwcs/coordinate_frames.py +++ b/gwcs/coordinate_frames.py @@ -436,6 +436,11 @@ def world_axis_object_classes(self): def world_axis_object_components(self): return [(f"{at}{i}" if i != 0 else at, 0, 'value') for i, at in enumerate(self._axes_type)] + @property + def _native_world_axis_object_components(self): + """Defines the target component ordering (i.e. not taking into account axes_order)""" + return self.world_axis_object_components + class CelestialFrame(CoordinateFrame): """ @@ -509,10 +514,18 @@ def world_axis_object_classes(self): 'unit': self.unit})} @property - def world_axis_object_components(self): + def _native_world_axis_object_components(self): return [('celestial', 0, 'spherical.lon'), ('celestial', 1, 'spherical.lat')] + @property + def world_axis_object_components(self): + # Sort the native waoc by the axes order. The axes order may have jumps + # in it if there are other frames in between the components. + ordered = np.array(self._native_world_axis_object_components, + dtype=object)[np.argsort(self.axes_order)] + return list(map(tuple, ordered)) + class SpectralFrame(CoordinateFrame): """ @@ -709,7 +722,7 @@ def _wao_renamed_components_iter(self): mapper = self._wao_classes_rename_map for frame in self.frames: renamed_components = [] - for comp in frame.world_axis_object_components: + for comp in frame._native_world_axis_object_components: comp = list(comp) rename = mapper[frame].get(comp[0]) if rename: diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 5276b937..139e0872 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -1272,9 +1272,14 @@ def test_spatial_spectral_stokes(): def test_split_frame_wcs(): - # Setup a model where the pixel & world axes are (lat, wave, lon) - spatial = models.Multiply(10*u.arcsec/u.pix) & models.Multiply(15*u.arcsec/u.pix) # pretend this is a spatial model + # Setup a WCS where the pixel & world axes are (lat, wave, lon) + + # We setup a model which is pretending to be a celestial transform. Note + # that we are pretending that this model is ordered lon, lat because that's + # what the projections require in astropy. + spatial = models.Multiply(20*u.arcsec/u.pix) & models.Multiply(15*u.arcsec/u.pix) compound = models.Linear1D(intercept=0*u.nm, slope=10*u.nm/u.pix) & spatial + # This forward transforms uses mappings to be (lat, wave, lon) forward = models.Mapping((1, 2, 0)) | compound | models.Mapping((2, 0, 1)) # Setup the output frame @@ -1287,14 +1292,52 @@ def test_split_frame_wcs(): axes_order=list(range(3)), unit=[u.pix]*3) iwcs = wcs.WCS(forward, input_frame, output_frame) - input_pixel = [0*u.pix, 1*u.pix, 2*u.pix] + input_pixel = [1*u.pix, 2*u.pix, 3*u.pix] output_world = iwcs.pixel_to_world_values(*input_pixel) output_pixel = iwcs.world_to_pixel_values(*output_world) assert_allclose(output_pixel, u.Quantity(input_pixel).to_value(u.pix)) + expected_world = [15*u.arcsec, 20*u.nm, 60*u.arcsec] + for expected, output in zip(expected_world, output_world): + assert_allclose(output, expected.value) + world_obj = iwcs.pixel_to_world(*input_pixel) assert isinstance(world_obj[0], coord.SkyCoord) assert isinstance(world_obj[1], coord.SpectralCoord) + assert u.allclose(world_obj[0].spherical.lat, expected_world[0]) + assert u.allclose(world_obj[0].spherical.lon, expected_world[2]) + assert u.allclose(world_obj[1], expected_world[1]) + obj_pixel = iwcs.world_to_pixel(*world_obj) assert_allclose(obj_pixel, u.Quantity(input_pixel).to_value(u.pix)) + + +def test_reordered_celestial(): + # This is a spatial model which is ordered lat, lon for the purposes of this test. + spatial = models.Multiply(20*u.deg/u.pix) & models.Multiply(15*u.deg/u.pix) + + celestial_frame = cf.CelestialFrame(axes_order=(1, 0), unit=(u.deg, u.deg), + reference_frame=coord.ICRS()) + + input_frame = cf.CoordinateFrame(2, ["PIXEL"]*2, + axes_order=list(range(2)), unit=[u.pix]*2) + + iwcs = wcs.WCS(spatial, input_frame, celestial_frame) + + input_pixel = [1*u.pix, 3*u.pix] + output_world = iwcs.pixel_to_world_values(*input_pixel) + output_pixel = iwcs.world_to_pixel_values(*output_world) + assert_allclose(output_pixel, u.Quantity(input_pixel).to_value(u.pix)) + + expected_world = [20*u.deg, 45*u.deg] + assert_allclose(output_world, [e.value for e in expected_world]) + + world_obj = iwcs.pixel_to_world(*input_pixel) + assert isinstance(world_obj, coord.SkyCoord) + + assert u.allclose(world_obj.spherical.lat, expected_world[0]) + assert u.allclose(world_obj.spherical.lon, expected_world[1]) + + obj_pixel = iwcs.world_to_pixel(world_obj) + assert_allclose(obj_pixel, u.Quantity(input_pixel).to_value(u.pix))