diff --git a/CHANGES.rst b/CHANGES.rst index 47b109e1..d9155f8d 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,18 @@ +0.15.0 (2020-11-13) +------------------- +New Features +^^^^^^^^^^^^ + +- Added ``insert_frame`` method to modify the pipeline of a ``WCS`` object. [#299] + +- Added ``to_fits_tab`` method to generate FITS header and binary table + extension following FITS WCS ``-TAB`` convension. [#295] + +- Added ``in_image`` function for testing whether a point in world coordinates + maps back to the domain of definition of the forward transformation. [#322] + +- Implemented iterative inverse for some imaging WCS. [#324] + 0.14.0 (2020-08-19) ------------------- New Features @@ -19,7 +34,6 @@ Bug Fixes - Add an optional parameter ``input_frame`` to ``wcstools.wcs_from_fiducial`. [#312] - 0.13.0 (2020-03-26) ------------------- New Features diff --git a/gwcs/coordinate_frames.py b/gwcs/coordinate_frames.py index 6a955024..b15f32ff 100644 --- a/gwcs/coordinate_frames.py +++ b/gwcs/coordinate_frames.py @@ -270,7 +270,6 @@ def axes_type(self): def coordinates(self, *args): """ Create world coordinates object""" - args = [args[i] for i in self.axes_order] coo = tuple([arg * un if not hasattr(arg, "to") else arg.to(un) for arg, un in zip(args, self.unit)]) return coo @@ -279,6 +278,10 @@ def coordinate_to_quantity(self, *coords): Given a rich coordinate object return an astropy quantity object. """ # NoOp leaves it to the model to handle + # If coords is a 1-tuple of quantity then return the element of the tuple + # This aligns the behavior with the other implementations + if not hasattr(coords, 'unit') and len(coords) == 1: + return coords[0] return coords @property @@ -286,12 +289,15 @@ def axis_physical_types(self): return self._axis_physical_types @property - def _world_axis_object_components(self): - raise NotImplementedError(f"This method is not implemented for {type(self)}") + def _world_axis_object_classes(self): + return {self._axes_type[0]: ( + u.Quantity, + (), + {'unit': self.unit[0]})} @property - def _world_axis_object_classes(self): - raise NotImplementedError(f"This method is not implemented for {type(self)}") + def _world_axis_object_components(self): + return [(self._axes_type[0], 0, 'value')] class CelestialFrame(CoordinateFrame): @@ -501,7 +507,12 @@ def _world_axis_object_classes(self): @property def _world_axis_object_components(self): - return [('temporal', 0, 'value')] + if isinstance(self.reference_frame.value, np.ndarray): + return [('temporal', 0, 'value')] + + def offset_from_time_and_reference(time): + return (time - self.reference_frame).sec + return [('temporal', 0, offset_from_time_and_reference)] def coordinates(self, *args): if np.isscalar(args): @@ -512,14 +523,16 @@ def coordinates(self, *args): return self._convert_to_time(dt, unit=self.unit[0], **self._attrs) def _convert_to_time(self, dt, *, unit, **kwargs): - if not isinstance(self.reference_frame.value, np.ndarray): - if not hasattr(dt, 'unit'): - dt = dt * unit - return self.reference_frame + dt - - else: + if (not isinstance(dt, time.TimeDelta) and + isinstance(dt, time.Time) or + isinstance(self.reference_frame.value, np.ndarray)): return time.Time(dt, **kwargs) + if not hasattr(dt, 'unit'): + dt = dt * unit + + return self.reference_frame + dt + def coordinate_to_quantity(self, *coords): if isinstance(coords[0], time.Time): ref_value = self.reference_frame.value @@ -532,6 +545,8 @@ def coordinate_to_quantity(self, *coords): # Is already a quantity elif hasattr(coords[0], 'unit'): return coords[0] + if isinstance(coords[0], np.ndarray): + return coords[0] * self.unit[0] else: raise ValueError("Can not convert {} to Quantity".format(coords[0])) diff --git a/gwcs/tests/conftest.py b/gwcs/tests/conftest.py index e0fa8085..e5efe34f 100644 --- a/gwcs/tests/conftest.py +++ b/gwcs/tests/conftest.py @@ -235,7 +235,7 @@ def sellmeier_zemax(): E_coef=E_coef) -@pytest.fixture +@pytest.fixture(scope="function") def gwcs_3d_galactic_spectral(): """ This fixture has the axes ordered as lat, spectral, lon. diff --git a/gwcs/tests/data/miri_lrs_wcs.asdf b/gwcs/tests/data/miri_lrs_wcs.asdf new file mode 100644 index 00000000..48d0e44d --- /dev/null +++ b/gwcs/tests/data/miri_lrs_wcs.asdf @@ -0,0 +1,955 @@ +#ASDF 1.0.0 +#ASDF_STANDARD 1.5.0 +%YAML 1.1 +%TAG ! tag:stsci.edu:asdf/ +--- !core/asdf-1.1.0 +asdf_library: !core/software-1.0.0 {author: Space Telescope Science Institute, homepage: 'http://github.com/spacetelescope/asdf', + name: asdf, version: 2.6.0} +history: + extensions: + - !core/extension_metadata-1.0.0 + extension_class: asdf.extension.BuiltinExtension + software: !core/software-1.0.0 {name: asdf, version: 2.6.0} + - !core/extension_metadata-1.0.0 + extension_class: gwcs.extension.GWCSExtension + software: !core/software-1.0.0 {name: gwcs, version: 0.13.1.dev16+g71f9b60} +wcs: ! + name: '' + steps: + - ! + frame: ! + axes_names: [x, y] + axes_order: [0, 1] + axis_physical_types: ['custom:x', 'custom:y'] + name: detector + unit: [!unit/unit-1.0.0 'pixel', !unit/unit-1.0.0 'pixel'] + transform: !transform/compose-1.2.0 + bounding_box: + - [6.5, 396.5] + - [302.5, 346.5] + bounds: + amplitude_7: [null, null] + angle_22: [null, null] + angle_3: [null, null] + c0_0_13: [null, null] + c0_0_14: [null, null] + c0_0_16: [null, null] + c0_0_17: [null, null] + c0_10: [null, null] + c0_11: [null, null] + c0_1_13: [null, null] + c0_1_14: [null, null] + c0_1_16: [null, null] + c0_1_17: [null, null] + c0_2_13: [null, null] + c0_2_14: [null, null] + c0_3_13: [null, null] + c0_3_14: [null, null] + c0_4_13: [null, null] + c0_4_14: [null, null] + c1_0_13: [null, null] + c1_0_14: [null, null] + c1_0_16: [null, null] + c1_0_17: [null, null] + c1_10: [null, null] + c1_11: [null, null] + c1_1_13: [null, null] + c1_1_14: [null, null] + c1_2_13: [null, null] + c1_2_14: [null, null] + c1_3_13: [null, null] + c1_3_14: [null, null] + c2_0_13: [null, null] + c2_0_14: [null, null] + c2_1_13: [null, null] + c2_1_14: [null, null] + c2_2_13: [null, null] + c2_2_14: [null, null] + c3_0_13: [null, null] + c3_0_14: [null, null] + c3_1_13: [null, null] + c3_1_14: [null, null] + c4_0_13: [null, null] + c4_0_14: [null, null] + offset_1: [null, null] + offset_2: [null, null] + offset_20: [null, null] + offset_21: [null, null] + offset_4: [null, null] + offset_5: [null, null] + offset_8: [null, null] + fixed: {amplitude_7: false, angle_22: false, angle_3: false, c0_0_13: false, + c0_0_14: false, c0_0_16: false, c0_0_17: false, c0_10: false, c0_11: false, + c0_1_13: false, c0_1_14: false, c0_1_16: false, c0_1_17: false, c0_2_13: false, + c0_2_14: false, c0_3_13: false, c0_3_14: false, c0_4_13: false, c0_4_14: false, + c1_0_13: false, c1_0_14: false, c1_0_16: false, c1_0_17: false, c1_10: false, + c1_11: false, c1_1_13: false, c1_1_14: false, c1_2_13: false, c1_2_14: false, + c1_3_13: false, c1_3_14: false, c2_0_13: false, c2_0_14: false, c2_1_13: false, + c2_1_14: false, c2_2_13: false, c2_2_14: false, c3_0_13: false, c3_0_14: false, + c3_1_13: false, c3_1_14: false, c4_0_13: false, c4_0_14: false, offset_1: false, + offset_2: false, offset_20: false, offset_21: false, offset_4: false, offset_5: false, + offset_8: false} + forward: + - !transform/remap_axes-1.2.0 + bounds: {} + fixed: {} + mapping: [0, 1, 0, 1] + - !transform/concatenate-1.2.0 + forward: + - !transform/compose-1.2.0 + forward: + - !transform/compose-1.2.0 + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/shift-1.2.0 + bounds: + offset: &id001 [null, null] + fixed: {offset: false} + offset: -325.13 + - !transform/shift-1.2.0 + bounds: + offset: *id001 + fixed: {offset: false} + offset: -299.7 + - !transform/rotate2d-1.2.0 + angle: 0.24155757363306068 + bounds: + angle: &id002 [null, null] + fixed: {angle: false} + - !transform/compose-1.2.0 + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/shift-1.2.0 + bounds: + offset: *id001 + fixed: {offset: false} + offset: 325.13 + - !transform/shift-1.2.0 + bounds: + offset: *id001 + fixed: {offset: false} + offset: 299.7 + - !transform/compose-1.2.0 + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/identity-1.2.0 + bounds: {} + fixed: {} + - !transform/constant-1.2.0 + bounds: + amplitude: [null, null] + fixed: {amplitude: false} + inverse: !transform/constant-1.2.0 + bounds: + amplitude: [null, null] + fixed: {amplitude: false} + value: 299.7 + value: 299.7 + - !transform/compose-1.2.0 + forward: + - !transform/compose-1.2.0 + forward: + - !transform/compose-1.2.0 + forward: + - !transform/compose-1.2.0 + forward: + - !transform/compose-1.2.0 + forward: + - !transform/compose-1.2.0 + forward: + - !transform/compose-1.2.0 + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/shift-1.2.0 + bounds: + offset: *id001 + fixed: {offset: false} + offset: -4.0 + - !transform/identity-1.2.0 + bounds: {} + fixed: {} + - !transform/concatenate-1.2.0 + forward: + - !transform/polynomial-1.2.0 + bounds: + c0: [null, null] + c1: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 0 + datatype: float64 + byteorder: little + shape: [2] + domain: [-1, 1] + fixed: {c0: false, c1: false} + inverse: !transform/polynomial-1.2.0 + bounds: + c0: [null, null] + c1: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 1 + datatype: float64 + byteorder: little + shape: [2] + domain: [-1, 1] + fixed: {c0: false, c1: false} + window: [-1, 1] + name: M_column_correction + window: [-1, 1] + - !transform/polynomial-1.2.0 + bounds: + c0: [null, null] + c1: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 2 + datatype: float64 + byteorder: little + shape: [2] + domain: [-1, 1] + fixed: {c0: false, c1: false} + inverse: !transform/polynomial-1.2.0 + bounds: + c0: [null, null] + c1: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 3 + datatype: float64 + byteorder: little + shape: [2] + domain: [-1, 1] + fixed: {c0: false, c1: false} + window: [-1, 1] + name: M_row_correction + window: [-1, 1] + - !transform/remap_axes-1.2.0 + bounds: {} + fixed: {} + inverse: !transform/identity-1.2.0 + bounds: {} + fixed: {} + n_dims: 2 + mapping: [0, 1, 0, 1] + - !transform/concatenate-1.2.0 + forward: + - !transform/polynomial-1.2.0 + bounds: + c0_0: [null, null] + c0_1: [null, null] + c0_2: [null, null] + c0_3: [null, null] + c0_4: [null, null] + c1_0: [null, null] + c1_1: [null, null] + c1_2: [null, null] + c1_3: [null, null] + c2_0: [null, null] + c2_1: [null, null] + c2_2: [null, null] + c3_0: [null, null] + c3_1: [null, null] + c4_0: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 4 + datatype: float64 + byteorder: little + shape: [5, 5] + domain: + - [-1, 1] + - [-1, 1] + fixed: {c0_0: false, c0_1: false, c0_2: false, c0_3: false, + c0_4: false, c1_0: false, c1_1: false, c1_2: false, + c1_3: false, c2_0: false, c2_1: false, c2_2: false, + c3_0: false, c3_1: false, c4_0: false} + inverse: !transform/polynomial-1.2.0 + bounds: + c0_0: [null, null] + c0_1: [null, null] + c0_2: [null, null] + c0_3: [null, null] + c0_4: [null, null] + c1_0: [null, null] + c1_1: [null, null] + c1_2: [null, null] + c1_3: [null, null] + c2_0: [null, null] + c2_1: [null, null] + c2_2: [null, null] + c3_0: [null, null] + c3_1: [null, null] + c4_0: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 5 + datatype: float64 + byteorder: little + shape: [5, 5] + domain: + - [-1, 1] + - [-1, 1] + fixed: {c0_0: false, c0_1: false, c0_2: false, c0_3: false, + c0_4: false, c1_0: false, c1_1: false, c1_2: false, + c1_3: false, c2_0: false, c2_1: false, c2_2: false, + c3_0: false, c3_1: false, c4_0: false} + window: + - [-1, 1] + - [-1, 1] + name: B_correction + window: + - [-1, 1] + - [-1, 1] + - !transform/polynomial-1.2.0 + bounds: + c0_0: [null, null] + c0_1: [null, null] + c0_2: [null, null] + c0_3: [null, null] + c0_4: [null, null] + c1_0: [null, null] + c1_1: [null, null] + c1_2: [null, null] + c1_3: [null, null] + c2_0: [null, null] + c2_1: [null, null] + c2_2: [null, null] + c3_0: [null, null] + c3_1: [null, null] + c4_0: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 6 + datatype: float64 + byteorder: little + shape: [5, 5] + domain: + - [-1, 1] + - [-1, 1] + fixed: {c0_0: false, c0_1: false, c0_2: false, c0_3: false, + c0_4: false, c1_0: false, c1_1: false, c1_2: false, + c1_3: false, c2_0: false, c2_1: false, c2_2: false, + c3_0: false, c3_1: false, c4_0: false} + inverse: !transform/polynomial-1.2.0 + bounds: + c0_0: [null, null] + c0_1: [null, null] + c0_2: [null, null] + c0_3: [null, null] + c0_4: [null, null] + c1_0: [null, null] + c1_1: [null, null] + c1_2: [null, null] + c1_3: [null, null] + c2_0: [null, null] + c2_1: [null, null] + c2_2: [null, null] + c3_0: [null, null] + c3_1: [null, null] + c4_0: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 7 + datatype: float64 + byteorder: little + shape: [5, 5] + domain: + - [-1, 1] + - [-1, 1] + fixed: {c0_0: false, c0_1: false, c0_2: false, c0_3: false, + c0_4: false, c1_0: false, c1_1: false, c1_2: false, + c1_3: false, c2_0: false, c2_1: false, c2_2: false, + c3_0: false, c3_1: false, c4_0: false} + window: + - [-1, 1] + - [-1, 1] + name: A_correction + window: + - [-1, 1] + - [-1, 1] + - !transform/remap_axes-1.2.0 + bounds: {} + fixed: {} + inverse: !transform/remap_axes-1.2.0 + bounds: {} + fixed: {} + mapping: [0, 1, 0, 1] + mapping: [0, 1, 0, 1] + - !transform/concatenate-1.2.0 + forward: + - !transform/polynomial-1.2.0 + bounds: + c0_0: [null, null] + c0_1: [null, null] + c1_0: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 8 + datatype: float64 + byteorder: little + shape: [2, 2] + domain: + - [-1, 1] + - [-1, 1] + fixed: {c0_0: false, c0_1: false, c1_0: false} + name: TI_row_correction + window: + - [-1, 1] + - [-1, 1] + - !transform/polynomial-1.2.0 + bounds: + c0_0: [null, null] + c0_1: [null, null] + c1_0: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 9 + datatype: float64 + byteorder: little + shape: [2, 2] + domain: + - [-1, 1] + - [-1, 1] + fixed: {c0_0: false, c0_1: false, c1_0: false} + name: TI_column_correction + window: + - [-1, 1] + - [-1, 1] + - !transform/identity-1.2.0 + bounds: {} + fixed: {} + inverse: !transform/remap_axes-1.2.0 + bounds: {} + fixed: {} + mapping: [0, 1, 0, 1] + n_dims: 2 + - !transform/remap_axes-1.2.0 + bounds: {} + fixed: {} + mapping: [1, 0] + - !transform/compose-1.2.0 + forward: + - !transform/compose-1.2.0 + forward: + - !transform/compose-1.2.0 + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/shift-1.2.0 + bounds: + offset: *id001 + fixed: {offset: false} + offset: -325.13 + - !transform/shift-1.2.0 + bounds: + offset: *id001 + fixed: {offset: false} + offset: -299.7 + - !transform/rotate2d-1.2.0 + angle: 0.24155757363306068 + bounds: + angle: *id002 + fixed: {angle: false} + - !transform/remap_axes-1.2.0 + bounds: {} + fixed: {} + mapping: [1] + - !transform/tabular-1.2.0 + bounding_box: [-291.60259000740217, 95.69787960086536] + bounds: {} + bounds_error: false + fill_value: .nan + fixed: {} + inverse: !transform/tabular-1.2.0 + bounding_box: [3.596040329703134, 14.387876373781149] + bounds: {} + bounds_error: false + fill_value: .nan + fixed: {} + lookup_table: !core/ndarray-1.0.0 + source: 12 + datatype: float64 + byteorder: little + shape: [699] + name: waverefinv + points: + - !core/ndarray-1.0.0 + source: 13 + datatype: float64 + byteorder: little + shape: [699] + lookup_table: !core/ndarray-1.0.0 + source: 10 + datatype: float64 + byteorder: little + shape: [388] + name: waveref + points: + - !core/ndarray-1.0.0 + source: 11 + datatype: float64 + byteorder: little + shape: [388] + inverse: !transform/compose-1.2.0 + bounds: + angle_14: [null, null] + angle_17: [null, null] + c0_0_2: [null, null] + c0_0_3: [null, null] + c0_0_5: [null, null] + c0_0_6: [null, null] + c0_1_2: [null, null] + c0_1_3: [null, null] + c0_1_5: [null, null] + c0_1_6: [null, null] + c0_2_5: [null, null] + c0_2_6: [null, null] + c0_3_5: [null, null] + c0_3_6: [null, null] + c0_4_5: [null, null] + c0_4_6: [null, null] + c0_8: [null, null] + c0_9: [null, null] + c1_0_2: [null, null] + c1_0_3: [null, null] + c1_0_5: [null, null] + c1_0_6: [null, null] + c1_1_5: [null, null] + c1_1_6: [null, null] + c1_2_5: [null, null] + c1_2_6: [null, null] + c1_3_5: [null, null] + c1_3_6: [null, null] + c1_8: [null, null] + c1_9: [null, null] + c2_0_5: [null, null] + c2_0_6: [null, null] + c2_1_5: [null, null] + c2_1_6: [null, null] + c2_2_5: [null, null] + c2_2_6: [null, null] + c3_0_5: [null, null] + c3_0_6: [null, null] + c3_1_5: [null, null] + c3_1_6: [null, null] + c4_0_5: [null, null] + c4_0_6: [null, null] + offset_10: [null, null] + offset_12: [null, null] + offset_13: [null, null] + offset_18: [null, null] + offset_19: [null, null] + fixed: {angle_14: false, angle_17: false, c0_0_2: false, c0_0_3: false, c0_0_5: false, + c0_0_6: false, c0_1_2: false, c0_1_3: false, c0_1_5: false, c0_1_6: false, + c0_2_5: false, c0_2_6: false, c0_3_5: false, c0_3_6: false, c0_4_5: false, + c0_4_6: false, c0_8: false, c0_9: false, c1_0_2: false, c1_0_3: false, c1_0_5: false, + c1_0_6: false, c1_1_5: false, c1_1_6: false, c1_2_5: false, c1_2_6: false, + c1_3_5: false, c1_3_6: false, c1_8: false, c1_9: false, c2_0_5: false, c2_0_6: false, + c2_1_5: false, c2_1_6: false, c2_2_5: false, c2_2_6: false, c3_0_5: false, + c3_0_6: false, c3_1_5: false, c3_1_6: false, c4_0_5: false, c4_0_6: false, + offset_10: false, offset_12: false, offset_13: false, offset_18: false, + offset_19: false} + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/compose-1.2.0 + forward: + - !transform/compose-1.2.0 + forward: + - !transform/compose-1.2.0 + forward: + - !transform/remap_axes-1.2.0 + bounds: {} + fixed: {} + mapping: [1, 0] + - !transform/compose-1.2.0 + forward: + - !transform/remap_axes-1.2.0 + bounds: {} + fixed: {} + mapping: [0, 1, 0, 1] + - !transform/compose-1.2.0 + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/polynomial-1.2.0 + bounds: + c0_0: [null, null] + c0_1: [null, null] + c1_0: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 14 + datatype: float64 + byteorder: little + shape: [2, 2] + domain: + - [-1, 1] + - [-1, 1] + fixed: {c0_0: false, c0_1: false, c1_0: false} + name: T_row_correction + window: + - [-1, 1] + - [-1, 1] + - !transform/polynomial-1.2.0 + bounds: + c0_0: [null, null] + c0_1: [null, null] + c1_0: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 15 + datatype: float64 + byteorder: little + shape: [2, 2] + domain: + - [-1, 1] + - [-1, 1] + fixed: {c0_0: false, c0_1: false, c1_0: false} + name: T_column_correction + window: + - [-1, 1] + - [-1, 1] + - !transform/compose-1.2.0 + forward: + - !transform/remap_axes-1.2.0 + bounds: {} + fixed: {} + mapping: [0, 1, 0, 1] + - !transform/compose-1.2.0 + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/polynomial-1.2.0 + bounds: + c0_0: [null, null] + c0_1: [null, null] + c0_2: [null, null] + c0_3: [null, null] + c0_4: [null, null] + c1_0: [null, null] + c1_1: [null, null] + c1_2: [null, null] + c1_3: [null, null] + c2_0: [null, null] + c2_1: [null, null] + c2_2: [null, null] + c3_0: [null, null] + c3_1: [null, null] + c4_0: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 16 + datatype: float64 + byteorder: little + shape: [5, 5] + domain: + - [-1, 1] + - [-1, 1] + fixed: {c0_0: false, c0_1: false, c0_2: false, c0_3: false, + c0_4: false, c1_0: false, c1_1: false, c1_2: false, + c1_3: false, c2_0: false, c2_1: false, c2_2: false, + c3_0: false, c3_1: false, c4_0: false} + window: + - [-1, 1] + - [-1, 1] + - !transform/polynomial-1.2.0 + bounds: + c0_0: [null, null] + c0_1: [null, null] + c0_2: [null, null] + c0_3: [null, null] + c0_4: [null, null] + c1_0: [null, null] + c1_1: [null, null] + c1_2: [null, null] + c1_3: [null, null] + c2_0: [null, null] + c2_1: [null, null] + c2_2: [null, null] + c3_0: [null, null] + c3_1: [null, null] + c4_0: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 17 + datatype: float64 + byteorder: little + shape: [5, 5] + domain: + - [-1, 1] + - [-1, 1] + fixed: {c0_0: false, c0_1: false, c0_2: false, c0_3: false, + c0_4: false, c1_0: false, c1_1: false, c1_2: false, + c1_3: false, c2_0: false, c2_1: false, c2_2: false, + c3_0: false, c3_1: false, c4_0: false} + window: + - [-1, 1] + - [-1, 1] + - !transform/compose-1.2.0 + forward: + - !transform/identity-1.2.0 + bounds: {} + fixed: {} + n_dims: 2 + - !transform/compose-1.2.0 + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/polynomial-1.2.0 + bounds: + c0: [null, null] + c1: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 18 + datatype: float64 + byteorder: little + shape: [2] + domain: [-1, 1] + fixed: {c0: false, c1: false} + window: [-1, 1] + - !transform/polynomial-1.2.0 + bounds: + c0: [null, null] + c1: [null, null] + coefficients: !core/ndarray-1.0.0 + source: 19 + datatype: float64 + byteorder: little + shape: [2] + domain: [-1, 1] + fixed: {c0: false, c1: false} + window: [-1, 1] + - !transform/concatenate-1.2.0 + forward: + - !transform/shift-1.2.0 + bounds: + offset: &id003 [null, null] + fixed: {offset: false} + offset: 4.0 + - !transform/identity-1.2.0 + bounds: {} + fixed: {} + - !transform/compose-1.2.0 + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/shift-1.2.0 + bounds: + offset: *id003 + fixed: {offset: false} + offset: -325.13 + - !transform/shift-1.2.0 + bounds: + offset: *id003 + fixed: {offset: false} + offset: -299.7 + - !transform/rotate2d-1.2.0 + angle: 0.24155757363306068 + bounds: + angle: &id004 [null, null] + fixed: {angle: false} + - !transform/remap_axes-1.2.0 + bounds: {} + fixed: {} + mapping: [0] + n_inputs: 2 + - !transform/tabular-1.2.0 + bounding_box: [3.596040329703134, 14.387876373781149] + bounds: {} + bounds_error: false + fill_value: .nan + fixed: {} + lookup_table: !core/ndarray-1.0.0 + source: 20 + datatype: float64 + byteorder: little + shape: [699] + name: waverefinv + points: + - !core/ndarray-1.0.0 + source: 21 + datatype: float64 + byteorder: little + shape: [699] + - !transform/compose-1.2.0 + forward: + - !transform/rotate2d-1.2.0 + angle: -0.24155757363306068 + bounds: + angle: *id004 + fixed: {angle: false} + - !transform/concatenate-1.2.0 + forward: + - !transform/shift-1.2.0 + bounds: + offset: *id003 + fixed: {offset: false} + offset: 325.13 + - !transform/shift-1.2.0 + bounds: + offset: *id003 + fixed: {offset: false} + offset: 299.7 + - ! + frame: ! + frames: + - ! + axes_names: [v2, v3] + axes_order: [0, 1] + axis_physical_types: ['custom:v2', 'custom:v3'] + name: v2v3_spatial + unit: [!unit/unit-1.0.0 'arcsec', !unit/unit-1.0.0 'arcsec'] + - ! + axes_names: [lambda] + axes_order: [2] + axis_physical_types: [em.wl] + name: spec + unit: [!unit/unit-1.0.0 'um'] + name: v2v3 + transform: !transform/concatenate-1.2.0 + bounds: + angles_2: [null, null] + factor_0: [null, null] + factor_1: [null, null] + fixed: {angles_2: false, factor_0: false, factor_1: false} + forward: + - !transform/compose-1.2.0 + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/scale-1.2.0 + bounds: + factor: &id005 [null, null] + factor: 0.0002777777777777778 + fixed: {factor: false} + - !transform/scale-1.2.0 + bounds: + factor: *id005 + factor: 0.0002777777777777778 + fixed: {factor: false} + - !transform/rotate_sequence_3d-1.0.0 + angles: [-0.12597594444444443, 0.10374517305555556, 0.0, -72.0545718, -5.630568] + axes_order: zyxyz + bounds: + angles: [null, null] + fixed: {angles: false} + name: v23tosky + rotation_type: spherical + - !transform/identity-1.2.0 + bounds: {} + fixed: {} + - ! + frame: ! + frames: + - ! + axes_names: [RA, DEC] + axes_order: [0, 1] + axis_physical_types: [pos.eq.ra, pos.eq.dec] + name: icrs + reference_frame: ! + frame_attributes: {} + unit: [!unit/unit-1.0.0 'deg', !unit/unit-1.0.0 'deg'] + - ! + axes_names: [lambda] + axes_order: [2] + axis_physical_types: [em.wl] + name: spec + unit: [!unit/unit-1.0.0 'um'] + name: world +... +BLK0۴ke-.u33333)?BLK0Eտ5Nh{@D@BLK0۴ke-.u33333)?BLK0Eտ5Nh{@D@BLK0 Z惂<LCjX?_@ r+9:mM)3><ǹ?"KmL-?@n.>TF>6 \UnҪx->Fa>V9>`f>BLK0*W˲VN{C&IiujBCƞ>@ Ո8x:?Z4~\W?rfƾd;P>6X=#>`hg=eLD-BLK0ȓr)ǸW`ӷi?ϝ_2*0IAn d>([RѨ@H`f0'9YE (p]>o-L'? t>3MX/Egy>\` ^>BLK09UHw o+u F_|0W4?Wʽxн{D=_=mR?w8j{?h  +X,,R[ތ>Iİ؃>%Mu,i"qF!Z9 +ؐ=BLK0 bdإZ*O^ +ףp=vwcCnqN?g$.?BLK0 ҍZoJ鄊zG!{g$.cCnqN?BLK0 U Lb~xR0h9-,@"\B,@RV,@ gބ,@DLh,@`c♐,@CF,@*,@a^x,@Fp,@BIzh,@Xw`,@mX,@`P,@tWH,@ <ځ@,@B7s8,@`a=a0,@a MN(,@AeW8 ,@V,@-j,@L),@ )+@r+@aƺ+@}AY+@&//+@}}+@y%+@n@!+@2jl+@HZ,4+@$+@Bٽ+@<~+@;+@s +@m٭+@*b|+@t+@|}k+@Coc+@:[+@B1R+@َaJ+@KuB+@LIE9+@'ޒ71+@{(+@%a +@'+@$)+@ZԬ+@-{s*@^?ծ*@r{*@,3*@[o*@J*@n p*@ +k*@ǘTK*@CRk*@X*@Sky*@%J֗*@b0*@mՆ*@W1}*@V,˜(u*@sl*@WIc*@Z*@m?R*@X7I|I*@H@*@d7*@/*@ǙG&*@kq*@Jd*@4G *@ 1&*@:)@L)@J)@")@t+)@('0)@?1)@0.)@ +.')@8oݑ)@fϭ )@w9 )@B8ߍ)@fs„)@^G|{)@Xm|r)@cxSRi)@V%$`)@/V)@EyM)@@~D)@=M>;)@IP1)@M/()@O~{c)@-)@${P? )@e^)@T-i(@Yt(@|F/(@IZ(@^z'N(@Hg((@I/XY(@8m׷(@'Q(@:f Ť(@b5(@Ȏ(@](@$ (@`a (@T~'@8_Q'@B'@;Vҷ'@h'@Y'@[B'@'@?` +'@?ah˖'@b}'@R,_'@oh;'@+5}'@Sr'@{a׭h'@5s^'@l2T'@I'@ݠ?'@wo0O5'@Q*'@TX '@V,{!7'@0/ + '@,+^'@3&@S#o&@z"&@/d.g&@ן/&@G&@>d&@-M&@NLi&@p޺&@qD &@#FT&@Z˕w&@MRl&@ﮑb&@y5W&@*^L&@AXA&@p36&@צ+&@/? &@:&@} O +&@4]%@ڪ %@tij%@ʅS%@k5s%@Q9PN%@<T"%@}G%@_"%@)tu%@`+-%@l߄%@Cy%@zA+n%@b%@"ZW%@oiK%@j@%@ 4%@{[)%@t%@>,%@l%@z$|$@&$@(x i$@$@$@;ta$@$@`i1$@>$@SB$@&sh7=$@91.x$@ԉl$@9}X_$@RS$@yoG$@Q;$@!o/$@}"$@Z5U$@{ $@IUIy#@D#@s#@#@A#@"*~#@j#@VV%#@;Z#@QdՄ#@ k~#@r*oq#@xŽd#@0VW#@ bJ#@A=#@!_d0#@/##@h#@ѡ#@%H"@!q["@n"@ =?"@T{b"@,ɷ"@Y)o""@ n"@"-|"@"Cހ"@)s"@yҎqe"@ZQ2"W"@_*I"@W@ +;"@bY,"@"]"@41}"@r1"@Nj!@ӆ؎o!@[|!@%s!@ ޹!@׫Y;!@0o!@+]2Ǎ!@;(~!@2Op!@3H(a!@X)R!@kC!@3!@=5$!@"2!@.C!@> @ x @m @k @m`˸ @.[P @&TW @_}!聉 @`Xy @MMi @5Y @ Q{iI @o09 @S( @c| @$E @@\iG@MEq@5R@6wkh@F~QF@h8$@&2@X՘@|"ջ@=n@eJu@&ͭQ@8-@W @%p@@@RJw@3Q@hѭ,@E@VA*@ +Pɹ@?!@&U2l@qWD@qԗ:pb@.*-@!S@]S@ދ@=FT@ZH@P@ZU˳ħ@; l@R|/@\ü@h@Ho@mŪ,@/o@w@V(X@|6$0 @/Dؿ@0Kn@E@4@h[@J5@Dr@nna@ 5oc@K`ʰ @BLK0 LSGiVSkk59r.L)r8,ٝrе rjy|q=Nq q5ĐqχqĊqhKRq3hq:q4 ~yqYziqfwYqᖁtIqwSq9q3hX%n)q+9jqfg qdpvlap1:>^p[pd}Wp^TpH?Qp1 XNp*KypdGipVDYpAIp/ރq>9pȡdC;)paeE8p(&4 p( rcoZ`]o琹VonR]PoJro%}դCRoWH=2oX6o0n3*n#!#nU^{n/ rnRn=f +2n!d +nSL&mQmZmjmi,=qmQQmw1mp(m2l olNlwl8[ql)Ql1lJ8}Fl~>kFk1kUՄkGDy~qk{cxQkq1kqdkkJejD ^jxOXjQjQKqj;EQjB$>1jt8j2W&2iڹ+i Am%i@țirO]qiYQi] 1i +塠i=lcDho$hzh/hiph;+wPhm0hhӥobg-1g7򩿰gk;MgupgI7Pg80g5Xܟgi{ff=$fnjfukf1pfcCyPf +Ws0fɑlffe/IB`ea' YẻSe5-MpeOFPe-Du@0e_:eR3dU`-d`'d*ا d\oKpd[Pd} 0d6d'cYb~c#"c!cioc$0h OcW)/c>TcŬbLnb"/@bU[㺏bⲇobit+Ob5ϧ/bxrbQaza <^aaoaOI{Oa*Bt/ana84h`a`LGH|[` U`UNo`܌gHO`dN B/`I;`䢥j_]%]_4Q^_&,|D_7^,2 +^X;R^^I8^"X]f>)]tp]]RD]\JGŝ\͎]\P֫\M[Ve[٬][}\x[;lZHc_Z!R]Z0iFZw>Y9YLo,YA[]Yiu0Y xwXr{X֔\XIP yObӛ_6O-6*FNR߹,6NoHM5MXfL5LK5K +zJJ92a5JVGIrP.5I"Hu(n4HA.G 44G;FA4Fk90f=E5e<(RSe;K^qe:Sjde9v1e8~d72d6kd5@%fd4޳C3d3kad2Qc1؝c0Rw.yi,p*,(V+h8&#C$R=\l"yt OGB-K)w|GI#) ? 9 ؿ0j?~T?r9a@I @욬G@bj4@ 9ރ@9D@B!@6ө C#@_úmC%@1P1C'@ +݉QD)@iqD+@X}E-@@AE/@00@·d'1@{FZ2@q[o(3@!c +4@GV5@J%6@u>X7@;28@&t9@ZV:@ǎ 8$;@/UW<@=@ݼ>@rܿ?@n7PS@@*@@2DSA@գ]A@CwSB@xɅB@SC@mgC@QSD@O4ID@TE@+)E@ݞBTF@$ \F@Y~uTG@G@j`TH@MH@/1zBTI@dtI@m$UJ@g'J@aAUK@:[wZK@pUsUL@gOYL@JIʦUM@.C;M@F=UN@{6N@0 VO@*%O@ O+P@@p,kP@B2 )9P@#EP@yR+Q@R_kQ@ +lQ@H|xQ@{+R@}v4kR@잫R@pR@N]+S@jkS@ѫS@wdS@h?+T@RZ^kT@K۰T@=XiT@#/!,U@ R*lU@Xϒ7U@LKDU@Q,V@)E]lV@tjV@^?-wV@,W@9lW@.VW@ɐ3W@BLK0؄tbc\ &j(W@eW@dW@ƾW@,`W@1W@x5W@3TW@ӑf/W@˗W@B+̚W@W@4y`W@!W@eH̨W@$^W@W)ٝW@-s:W@Q$W@UW@y:†W@5{`W@V4nzW@QM(RtW@ mW@&vgW@6q1`W@p>g ZW@qUESW@6?JLW@/|.EW@=Y=W@W6W@#$ /W@*Zf'W@6ZW@zW@R*9W@wZW@TV@stWoV@&lV@V@tV@v}SV@˯fV@{FV@WNV@F퉬V@<7V@{=}R@jR@UMP{R@5*DiR@WR@$)ER@y~3R@`!R@s^R@Q@ztBQ@#Q@)gkQ@hQ@/V@Q@DQ@-{Q@[8hQ@JUQ@YHxBQ@ʹ/Q@hӠQ@AuM Q@n(P@[}P@P@eXP@P@tHޓP@aXP@>rUlP@iXP@ZCP@t/P@ /nP@YP@`99O@Z5'O@ +O@WZÀiO@y?O@Ԉb)O@pt=N@2m%N@ܖN@.=!rmN@窈BN@+N@ +CM@ rM@lM@X7kM@mw?M@ M@|L@r3PܛL@rrGL@{'bL@ 6L@XqG L@@\$GK@KmK@܍ǁK@vBHTK@#d&K@M6 J@ιJ@%]J@x9anJ@?J@zMJ@LyWI@UI@fkI@0]FUI@?i%I@>WH@K]H@.͖H@GRwefH@hũ6H@AG\_H@NWVeG@'ZG@mtG@|}CG@ }G@lF@㽎EF@,F~F@~:MF@v}F@"ʭE@cbaE@@tdE@_qSE@|I E@ǬD@J2D@[b&FD@!R.4UD@n1!D@>JC@FC@yikC@qpSC@坥gC@lO}B@ 5B@$o*łB@.NB@%pB@0ϢA@f5"Y|A@)FzA@r`:DA@XcA@o@@Sߣ@@ m@@ v7@@B k@@?@5uSң(?@I>@9)M>@)=@fp=@32a@=@vX<@"<@;@Чy B;@l)T:@d܁a:@+NZ9@,?p~9@|> 9@*8@-,o(8@T+7@pzB7@Ii6@Eq.[6@bP5@Ar5@%V4@brF4@<4@=PD3@]l'3@2@92@1@r +GK1@0@ܢh[0@ +/@>Z~.@: e-@ ,@(%D+@)@+@6̛*@}V)@z)/q((@~pq1'@2 /:&@e_C%@%KJ$@1e}Q#@jX"@$W]!@'- +c @*=@!1@= @s@|1P@V"@'7@N5@w ` @ @gy@%m@ʠ%n?=h?? \L?pu$ӿRwOy3W~К-Ai/WiD! +pb)DuXp[CFu _Z +?#[˱<j6v+ 8,9!Q79H"@TW#pGg$w Lw%i)&i ᝙'Z^t(L)w*[+&u),M<.٫j$/[ppQ00jp41t1:L24e2?FRf3Tw3EV4o5FRt5/+6y!O64AcI7vUU7jnxh8kX{88͈9p:b:vϥ;;C9;`"_<SL#ܩ>=o(xW=?j$?l2@xq|@3N<@ТAʀf|\A;ZAscbA=Br*Bؗ)tBĞ# Ch kCCZ,D[PD`D,~ZTF؞lkF'=F=GZUGGGϡGjL'@HyhU؎H>Hhf,IK{I"QI>nJ jJփtxJe +KCZKzK)K2*KK~SR T[V7TkcT#tTX7ٺT`p-UT$>UY>Uu.jU쨄זU5yUG=UfUV GHV҆tVxV8/V/BVk)'WT$TWӁWW(1W)^ Xt6X>dXJZX=K )X&X\AYoLHYn)vYȷYmn0Y0E7ZB$c.Zxl\ZO՟ZpZ<\Z7BZqq]e[ &`E[Ҿs[enQA[?g[适;Z\5`/\aɣ{^\ ʛ\06\ۯHE\ ]6S81K]z]8ſn]wh7,]hx4 ^&9^ xUi^oø^b7^rafqM^cp*_@Z_0Ju_K{_-_3;`s'`ț;D?`+W`Vp` `ce ` `S`H`.n޽`3~&aOq&Saf-6axDcOaXƪha3oa)a.zea3~1a/Wa5aoQ9bu~ 1b-^JbzrcbIG2$}bK/ubW?bVb~{3b +bv}cA30cQIc cc}c+#xcDTec^cqs/ac\pckdU/V3d!2$Md8Xrhd&wad0ղdUdT!r#xd\dqje& eԉ;e +r*VeB>u^pee%Pe<e+[XgR.tgVg##4gJl'g-<נgfKgWAh\M/e5hapQhX`mhzpp#Mp*b]p lp*'\pm|p7pɄpR +!p氤p^Zp t(pIpxS6ZTp Ό q$q>1v(q=58qG[GqaMWq75gq R^wq(}V3qum q +A@q3̶qͧ)q$]lqيq j}qOtr$nr<6dm&rƈp6rBLK0j:)O=!.K`ʰ @MZ @Ӿ @5# @[YC @Lc @{' @'7W @kzF @ U @8ceU@7ht @{K@@ƿS`@3@G%@TQ@σ@@WxO@=@@a]@#pN}@gI.@ν=@1ML@3-\@<. 6 @E=@]Ez-@"u;MO=@ČT$M@f\\@idl@#l|@Lsx@{M@R"@22 @Iƒ@va@y:v@K@\ @h@"*@Bȟ:@tJ@QIZ@(6 j@My@leȉ@}9@r@RH@g@!@8@ +%@|"P-q@: +5F@Q<(@bi~D7@8LG@SW@H[og@fcDw@ k@.r@zæ@r&O@> m@UÑB@Xm}@7@@>@ek%@@5@$E@T@h*N׿d@ +Bߔt@Yi@Nq|>@6@@4@d @xh@=@$@^.M,@F4"@];2@Du{CfB@5K;R@Rb@*Zq@cb@nj@qd@y9@T2L@I@a@:yz@ܐ4b@~7@ @b/@d?@ƋO@`_@J6K5o@M +@eߎ@0}y@Ҕ3@t^@ħ4@a @Z@ +@"# @@:J+]@Q32-@i:=@&xBL@Ș2J\@jQl@ ȦY[|@`a0@Pi@pګ@&x@6>I@UY@zm.@w@1 +@`릭@̥*@_W:@F,J@Z@*i@,BHիy@Y݀@pqU@v*@0@V@Ϥ@^ @<T@)@."@"FG*'@]27@fu9}G@uARW@/I'g@LPv@ӣXц@]`@2h{@oP@v2w%@JF@a@\y@ty@.N@B#%@ע4@\D@(T@wd@l6Lt@NE!@e@R}ˣ@s렳@-u@8J@ۡ @|[ +@ @"@b:!t"@RD)I2@i0B@H8Q@r@a@,Hq@.Or@ߠWG@rZ_@g@&n@X>v@UC~p@mE@>@q@+@$夙/@㟬n?@hYCO@ +_@*n@NB~@YBӗ@ql@4A@֠p@x*@@@^X k@@ @.@DF ,@]A(<@u/L@*7i\@̤o?>l@n)G|@N@V@TX^@fg@2m<@:Ju@a@}@~y@ +@¨ne@d(:*@:@I@JWY@i@6cy@0N8@e? @t}@᷸@m@Z'a@6@ @@ V@"@:'@&R`7@i>'5G@j. +W@ 6f@l>v@P&F@M^@U3@6U]@&e@z>l@Vt@m=|\@`1@@k$@F%4@ߢD@ZT@,T/d@*t@pBك@Zɮ@q<у@VX@-@j@<$@@@ @+ @b @3c @@ @պ( @0 @wך8 @H޴@ @ꑢH @n&jP @LTX @ ).?` @])h @.%5p @0w @<= @HzӇ @rTWE @C`4ɨ @lM @w} @Th @R @X\= @)b' @?d @˾ @k @m @>s @ @m{ @K{!@(f!@SQ!@$);!@4&&!@@.!@Ly5!@hXV=!@9d3E!@ +pM!@{%U!@ʩ]!@}-ze!@Ndm!@a5Ou!@>9}!@=$!@!@cD!@4!@LΤ!@lи!@ JT!@x'؍!@I!\x!@-b!@8cM!@D7!@Pxk"!@^\U !@/h2s!@t!@z"@ "@s"@D"@`v#"@=a+"@K3"@6;"@Yԙ C"@* K"@R"@l%Z"@ Ib"@n&-j"@?%r"@14z"@$@Z)$@\$@9b$@$@`i$@1$@q%@ %@hy}%@uEg%@F)"R#%@5=+%@@܈'3%@L ;%@XB%@[dsJ%@,pPR%@{-Z%@· +b%@#j%@pħ{r%@A+fz%@~P%@[3;%@8%%@;%@V%@'B%@ϱ%@ J%@gΤ%@k!DR%@<-!y%@ 9Yd%@DN%@Pa9%@\#%@Qhri%@"tO&@,q&@ċ &@x&@f &@7(&@}x0&@Zb8&@7 M@&@{7H&@L"P&@Η X&@_&@ g&@f#o&@a%Cw&@21 +&@=&@H2v&@T`&@v`:K&@Glq5&@xNB &@+ +&@J&@&@\Q&@-մ&@|Y&@Y݉&@6at&@q^&@BhI&@3'@p'@'@ex'@W)B%'@(5-'@@5'@Lه='@X E'@ldrM'@=pp]U'@|MG]'@߇*2e'@m'@"u'@R|'@#*܄'@{ƌ'@X2'@5'@g:'@8p'@ A['@ E'@I0'@|!d'@M-AQ'@9'@DX'@P'@\`'@bh'@3toh(@Ln (@Ջ)pY(@C(@ww.#(@H+(@3(@z:(@WB(@4 J(@]R(@.Z(@̖b(@ mj(@Wr(@r%c"Bz(@C1@,(@=*(@H(@T1(@`֡(@Xl9(@)xn(@KA(@ˏ(ŀ(@Ik(@mU(@>P@(@*(@yX(@V(@3`(@S)@$g)@)@o)@~ )@h)bwi()@95?S0)@ +A>8)@L)@)@XֆH)@}d +O)@NpW)@|m_)@Jg)@'o)@w)@c!})@4g)@Û)R)@x<)@U1')@x2)@I9)@)@ @Ѿ)@Ļ)@!H)@^-a̐)@/9>P{)@Ee)@PWP)@\:)@sh_%)@Dt*@lg *@I*@&o*@%*@Yv-*@*5*@ƚ~y=*@wdE*@TNM*@n1 +9U*@?#]*@e*@ ɕl*@t*@%|*@T1`!*@%==*@H)*@Tw*@`0b*@ilL*@:x87*@ k!*@܏H@ *@%*@~H*@O*@ O*@ʙӠ*@vW*@Su*@d0_`+@5 J +@f5+@+@n +#+@y)*+@J5_v2+@A<:+@L~B+@XJ+@dӅR+@_p tZ+@0|^b+@jIj+@ғG3r+@$z+@t+@E +@ûݑ+@Θ(ș+@u+@R0+@Z/+@+ 8r+@ \+@?G+@!1+@o-G+@@9^+@E;O+@P+@\V+@hڰ,@Ut^,@&,@ifp,@ȗFZ ,@#nE(,@j/0,@;u8,@ Ǻ@,@җ}G,@tO,@QW,@P. _,@! g,@ o,@Ɣnw,@%Y,@e1C,@6=] .,@I:,@T(,@`,@zl/خ,@Kx¶,@7,@h,@BLK0 UϜ4q&+C- Fy@eCnqN?h$.?BLK0 Maxd!&޺]o>/.yh$.cCnqN?BLK09UHw o+u F_|0W4?Wʽxн{D=_=mR?w8j{?h  +X,,R[ތ>Iİ؃>%Mu,i"qF!Z9 +ؐ=BLK0*W˲VN{C&IiujBCƞ>@ Ո8x:?Z4~\W?rfƾd;P>6X=#>`hg=eLD-BLK0Eտ5Nh{@D@BLK0Eտ5Nh{@D@BLK0؄tbc\ &j(W@eW@dW@ƾW@,`W@1W@x5W@3TW@ӑf/W@˗W@B+̚W@W@4y`W@!W@eH̨W@$^W@W)ٝW@-s:W@Q$W@UW@y:†W@5{`W@V4nzW@QM(RtW@ mW@&vgW@6q1`W@p>g ZW@qUESW@6?JLW@/|.EW@=Y=W@W6W@#$ /W@*Zf'W@6ZW@zW@R*9W@wZW@TV@stWoV@&lV@V@tV@v}SV@˯fV@{FV@WNV@F퉬V@<7V@{=}R@jR@UMP{R@5*DiR@WR@$)ER@y~3R@`!R@s^R@Q@ztBQ@#Q@)gkQ@hQ@/V@Q@DQ@-{Q@[8hQ@JUQ@YHxBQ@ʹ/Q@hӠQ@AuM Q@n(P@[}P@P@eXP@P@tHޓP@aXP@>rUlP@iXP@ZCP@t/P@ /nP@YP@`99O@Z5'O@ +O@WZÀiO@y?O@Ԉb)O@pt=N@2m%N@ܖN@.=!rmN@窈BN@+N@ +CM@ rM@lM@X7kM@mw?M@ M@|L@r3PܛL@rrGL@{'bL@ 6L@XqG L@@\$GK@KmK@܍ǁK@vBHTK@#d&K@M6 J@ιJ@%]J@x9anJ@?J@zMJ@LyWI@UI@fkI@0]FUI@?i%I@>WH@K]H@.͖H@GRwefH@hũ6H@AG\_H@NWVeG@'ZG@mtG@|}CG@ }G@lF@㽎EF@,F~F@~:MF@v}F@"ʭE@cbaE@@tdE@_qSE@|I E@ǬD@J2D@[b&FD@!R.4UD@n1!D@>JC@FC@yikC@qpSC@坥gC@lO}B@ 5B@$o*łB@.NB@%pB@0ϢA@f5"Y|A@)FzA@r`:DA@XcA@o@@Sߣ@@ m@@ v7@@B k@@?@5uSң(?@I>@9)M>@)=@fp=@32a@=@vX<@"<@;@Чy B;@l)T:@d܁a:@+NZ9@,?p~9@|> 9@*8@-,o(8@T+7@pzB7@Ii6@Eq.[6@bP5@Ar5@%V4@brF4@<4@=PD3@]l'3@2@92@1@r +GK1@0@ܢh[0@ +/@>Z~.@: e-@ ,@(%D+@)@+@6̛*@}V)@z)/q((@~pq1'@2 /:&@e_C%@%KJ$@1e}Q#@jX"@$W]!@'- +c @*=@!1@= @s@|1P@V"@'7@N5@w ` @ @gy@%m@ʠ%n?=h?? \L?pu$ӿRwOy3W~К-Ai/WiD! +pb)DuXp[CFu _Z +?#[˱<j6v+ 8,9!Q79H"@TW#pGg$w Lw%i)&i ᝙'Z^t(L)w*[+&u),M<.٫j$/[ppQ00jp41t1:L24e2?FRf3Tw3EV4o5FRt5/+6y!O64AcI7vUU7jnxh8kX{88͈9p:b:vϥ;;C9;`"_<SL#ܩ>=o(xW=?j$?l2@xq|@3N<@ТAʀf|\A;ZAscbA=Br*Bؗ)tBĞ# Ch kCCZ,D[PD`D,~ZTF؞lkF'=F=GZUGGGϡGjL'@HyhU؎H>Hhf,IK{I"QI>nJ jJփtxJe +KCZKzK)K2*KK~SR T[V7TkcT#tTX7ٺT`p-UT$>UY>Uu.jU쨄זU5yUG=UfUV GHV҆tVxV8/V/BVk)'WT$TWӁWW(1W)^ Xt6X>dXJZX=K )X&X\AYoLHYn)vYȷYmn0Y0E7ZB$c.Zxl\ZO՟ZpZ<\Z7BZqq]e[ &`E[Ҿs[enQA[?g[适;Z\5`/\aɣ{^\ ʛ\06\ۯHE\ ]6S81K]z]8ſn]wh7,]hx4 ^&9^ xUi^oø^b7^rafqM^cp*_@Z_0Ju_K{_-_3;`s'`ț;D?`+W`Vp` `ce ` `S`H`.n޽`3~&aOq&Saf-6axDcOaXƪha3oa)a.zea3~1a/Wa5aoQ9bu~ 1b-^JbzrcbIG2$}bK/ubW?bVb~{3b +bv}cA30cQIc cc}c+#xcDTec^cqs/ac\pckdU/V3d!2$Md8Xrhd&wad0ղdUdT!r#xd\dqje& eԉ;e +r*VeB>u^pee%Pe<e+[XgR.tgVg##4gJl'g-<נgfKgWAh\M/e5hapQhX`mhzpp#Mp*b]p lp*'\pm|p7pɄpR +!p氤p^Zp t(pIpxS6ZTp Ό q$q>1v(q=58qG[GqaMWq75gq R^wq(}V3qum q +A@q3̶qͧ)q$]lqيq j}qOtr$nr<6dm&rƈp6rBLK0j:)O=!.K`ʰ @MZ @Ӿ @5# @[YC @Lc @{' @'7W @kzF @ U @8ceU@7ht @{K@@ƿS`@3@G%@TQ@σ@@WxO@=@@a]@#pN}@gI.@ν=@1ML@3-\@<. 6 @E=@]Ez-@"u;MO=@ČT$M@f\\@idl@#l|@Lsx@{M@R"@22 @Iƒ@va@y:v@K@\ @h@"*@Bȟ:@tJ@QIZ@(6 j@My@leȉ@}9@r@RH@g@!@8@ +%@|"P-q@: +5F@Q<(@bi~D7@8LG@SW@H[og@fcDw@ k@.r@zæ@r&O@> m@UÑB@Xm}@7@@>@ek%@@5@$E@T@h*N׿d@ +Bߔt@Yi@Nq|>@6@@4@d @xh@=@$@^.M,@F4"@];2@Du{CfB@5K;R@Rb@*Zq@cb@nj@qd@y9@T2L@I@a@:yz@ܐ4b@~7@ @b/@d?@ƋO@`_@J6K5o@M +@eߎ@0}y@Ҕ3@t^@ħ4@a @Z@ +@"# @@:J+]@Q32-@i:=@&xBL@Ș2J\@jQl@ ȦY[|@`a0@Pi@pګ@&x@6>I@UY@zm.@w@1 +@`릭@̥*@_W:@F,J@Z@*i@,BHիy@Y݀@pqU@v*@0@V@Ϥ@^ @<T@)@."@"FG*'@]27@fu9}G@uARW@/I'g@LPv@ӣXц@]`@2h{@oP@v2w%@JF@a@\y@ty@.N@B#%@ע4@\D@(T@wd@l6Lt@NE!@e@R}ˣ@s렳@-u@8J@ۡ @|[ +@ @"@b:!t"@RD)I2@i0B@H8Q@r@a@,Hq@.Or@ߠWG@rZ_@g@&n@X>v@UC~p@mE@>@q@+@$夙/@㟬n?@hYCO@ +_@*n@NB~@YBӗ@ql@4A@֠p@x*@@@^X k@@ @.@DF ,@]A(<@u/L@*7i\@̤o?>l@n)G|@N@V@TX^@fg@2m<@:Ju@a@}@~y@ +@¨ne@d(:*@:@I@JWY@i@6cy@0N8@e? @t}@᷸@m@Z'a@6@ @@ V@"@:'@&R`7@i>'5G@j. +W@ 6f@l>v@P&F@M^@U3@6U]@&e@z>l@Vt@m=|\@`1@@k$@F%4@ߢD@ZT@,T/d@*t@pBك@Zɮ@q<у@VX@-@j@<$@@@ @+ @b @3c @@ @պ( @0 @wך8 @H޴@ @ꑢH @n&jP @LTX @ ).?` @])h @.%5p @0w @<= @HzӇ @rTWE @C`4ɨ @lM @w} @Th @R @X\= @)b' @?d @˾ @k @m @>s @ @m{ @K{!@(f!@SQ!@$);!@4&&!@@.!@Ly5!@hXV=!@9d3E!@ +pM!@{%U!@ʩ]!@}-ze!@Ndm!@a5Ou!@>9}!@=$!@!@cD!@4!@LΤ!@lи!@ JT!@x'؍!@I!\x!@-b!@8cM!@D7!@Pxk"!@^\U !@/h2s!@t!@z"@ "@s"@D"@`v#"@=a+"@K3"@6;"@Yԙ C"@* K"@R"@l%Z"@ Ib"@n&-j"@?%r"@14z"@$@Z)$@\$@9b$@$@`i$@1$@q%@ %@hy}%@uEg%@F)"R#%@5=+%@@܈'3%@L ;%@XB%@[dsJ%@,pPR%@{-Z%@· +b%@#j%@pħ{r%@A+fz%@~P%@[3;%@8%%@;%@V%@'B%@ϱ%@ J%@gΤ%@k!DR%@<-!y%@ 9Yd%@DN%@Pa9%@\#%@Qhri%@"tO&@,q&@ċ &@x&@f &@7(&@}x0&@Zb8&@7 M@&@{7H&@L"P&@Η X&@_&@ g&@f#o&@a%Cw&@21 +&@=&@H2v&@T`&@v`:K&@Glq5&@xNB &@+ +&@J&@&@\Q&@-մ&@|Y&@Y݉&@6at&@q^&@BhI&@3'@p'@'@ex'@W)B%'@(5-'@@5'@Lه='@X E'@ldrM'@=pp]U'@|MG]'@߇*2e'@m'@"u'@R|'@#*܄'@{ƌ'@X2'@5'@g:'@8p'@ A['@ E'@I0'@|!d'@M-AQ'@9'@DX'@P'@\`'@bh'@3toh(@Ln (@Ջ)pY(@C(@ww.#(@H+(@3(@z:(@WB(@4 J(@]R(@.Z(@̖b(@ mj(@Wr(@r%c"Bz(@C1@,(@=*(@H(@T1(@`֡(@Xl9(@)xn(@KA(@ˏ(ŀ(@Ik(@mU(@>P@(@*(@yX(@V(@3`(@S)@$g)@)@o)@~ )@h)bwi()@95?S0)@ +A>8)@L)@)@XֆH)@}d +O)@NpW)@|m_)@Jg)@'o)@w)@c!})@4g)@Û)R)@x<)@U1')@x2)@I9)@)@ @Ѿ)@Ļ)@!H)@^-a̐)@/9>P{)@Ee)@PWP)@\:)@sh_%)@Dt*@lg *@I*@&o*@%*@Yv-*@*5*@ƚ~y=*@wdE*@TNM*@n1 +9U*@?#]*@e*@ ɕl*@t*@%|*@T1`!*@%==*@H)*@Tw*@`0b*@ilL*@:x87*@ k!*@܏H@ *@%*@~H*@O*@ O*@ʙӠ*@vW*@Su*@d0_`+@5 J +@f5+@+@n +#+@y)*+@J5_v2+@A<:+@L~B+@XJ+@dӅR+@_p tZ+@0|^b+@jIj+@ғG3r+@$z+@t+@E +@ûݑ+@Θ(ș+@u+@R0+@Z/+@+ 8r+@ \+@?G+@!1+@o-G+@@9^+@E;O+@P+@\V+@hڰ,@Ut^,@&,@ifp,@ȗFZ ,@#nE(,@j/0,@;u8,@ Ǻ@,@җ}G,@tO,@QW,@P. _,@! g,@ o,@Ɣnw,@%Y,@e1C,@6=] .,@I:,@T(,@`,@zl/خ,@Kx¶,@7,@h,@#ASDF BLOCK INDEX +%YAML 1.1 +--- +- 35600 +- 35670 +- 35740 +- 35810 +- 35880 +- 36134 +- 36388 +- 36642 +- 36896 +- 36982 +- 37068 +- 40226 +- 43384 +- 49030 +- 54676 +- 54762 +- 54848 +- 55102 +- 55356 +- 55426 +- 55496 +- 61142 +... diff --git a/gwcs/tests/data/nircamwcs.asdf b/gwcs/tests/data/nircamwcs.asdf new file mode 100644 index 00000000..f0e2c14b Binary files /dev/null and b/gwcs/tests/data/nircamwcs.asdf differ diff --git a/gwcs/tests/test_api.py b/gwcs/tests/test_api.py index 4942a8b3..3635a35d 100644 --- a/gwcs/tests/test_api.py +++ b/gwcs/tests/test_api.py @@ -173,10 +173,10 @@ def test_world_axis_object_components_1d(gwcs_1d_freq): def test_world_axis_object_components_4d(gwcs_4d_identity_units): waoc = gwcs_4d_identity_units.world_axis_object_components - assert waoc == [('celestial', 0, 'spherical.lon'), - ('celestial', 1, 'spherical.lat'), - ('spectral', 0, 'value'), - ('temporal', 0, 'value')] + assert waoc[0:3] == [('celestial', 0, 'spherical.lon'), + ('celestial', 1, 'spherical.lat'), + ('spectral', 0, 'value')] + assert waoc[3][0:2] == ('temporal', 0) def test_world_axis_object_classes_2d(gwcs_2d_spatial_shift): diff --git a/gwcs/tests/test_coordinate_systems.py b/gwcs/tests/test_coordinate_systems.py index 1f5df866..056aea1b 100644 --- a/gwcs/tests/test_coordinate_systems.py +++ b/gwcs/tests/test_coordinate_systems.py @@ -4,10 +4,12 @@ from numpy.testing import assert_allclose import astropy.units as u -from astropy.time import Time +from astropy.time import Time, TimeDelta from astropy import coordinates as coord from astropy.tests.helper import assert_quantity_allclose +from astropy.modeling import models as m +from .. import WCS from .. import coordinate_frames as cf import astropy @@ -86,6 +88,30 @@ def test_coordinates_composite(inputs): assert_allclose(result[1].value, inputs[2]) +def test_coordinates_composite_order(): + time = cf.TemporalFrame(Time("2011-01-01T00:00:00"), name='time', unit=[u.s, ], axes_order=(0, )) + dist = cf.CoordinateFrame(name='distance', naxes=1, + axes_type=["SPATIAL"], unit=[u.m, ], axes_order=(1, )) + frame = cf.CompositeFrame([time, dist]) + result = frame.coordinates(0, 0) + assert result[0] == Time("2011-01-01T00:00:00") + assert u.allclose(result[1], 0*u.m) + + +def test_bare_baseframe(): + # This is a regression test for the following call: + frame = cf.CoordinateFrame(1, "SPATIAL", (0,), unit=(u.km,)) + assert u.allclose(frame.coordinate_to_quantity((1*u.m,)), 1*u.m) + + # Now also setup the same situation through the whole call stack to be safe. + w = WCS(forward_transform=m.Tabular1D(points=np.arange(10)*u.pix, + lookup_table=np.arange(10)*u.km), + output_frame=frame, + input_frame=cf.CoordinateFrame(1, "PIXEL", (0,), unit=(u.pix,), name="detector_frame") + ) + assert u.allclose(w.world_to_pixel(0*u.km), 0) + + @pytest.mark.parametrize(('frame'), coord_frames) def test_celestial_attributes_length(frame): """ @@ -138,6 +164,10 @@ def test_base_coordinate(): assert_quantity_allclose(q1, 12 * u.deg) assert_quantity_allclose(q2, 3 * u.arcsec) + q1, q2 = frame.coordinate_to_quantity((12 * u.deg, 3 * u.arcsec)) + assert_quantity_allclose(q1, 12 * u.deg) + assert_quantity_allclose(q2, 3 * u.arcsec) + def test_temporal_relative(): t = cf.TemporalFrame(reference_frame=Time("2018-01-01T00:00:00"), unit=u.s) @@ -150,6 +180,7 @@ def test_temporal_relative(): t = cf.TemporalFrame(reference_frame=Time("2018-01-01T00:00:00")) assert t.coordinates(10 * u.s) == Time("2018-01-01T00:00:00") + 10 * u.s + assert t.coordinates(TimeDelta(10, format='sec')) == Time("2018-01-01T00:00:00") + 10 * u.s a = t.coordinates((10, 20) * u.s) assert a[0] == Time("2018-01-01T00:00:00") + 10 * u.s diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index f53a6417..e72f1137 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -9,6 +9,8 @@ import pytest from astropy.utils.data import get_pkg_data_filename from astropy import wcs as astwcs +from astropy.wcs import wcsapi +from astropy.time import Time from .. import wcs from ..wcstools import (wcs_from_fiducial, grid_from_bounding_box, wcs_from_points) @@ -17,6 +19,7 @@ from ..utils import CoordinateFrameError import asdf + m1 = models.Shift(12.4) & models.Shift(-2) m2 = models.Scale(2) & models.Scale(-2) m = m1 | m2 @@ -25,6 +28,7 @@ detector = cf.Frame2D(name='detector', axes_order=(0, 1)) focal = cf.Frame2D(name='focal', axes_order=(0, 1), unit=(u.m, u.m)) spec = cf.SpectralFrame(name='wave', unit=[u.m, ], axes_order=(2, ), axes_names=('lambda', )) +time = cf.TemporalFrame(name='time', unit=[u.s, ], axes_order=(3, ), axes_names=('time', ), reference_frame=Time("2020-01-01")) pipe = [(detector, m1), (focal, m2), @@ -99,6 +103,36 @@ def test_insert_transform(): assert_allclose(gw.forward_transform(1, 2), (m1 | m2)(1, 2)) +def test_insert_frame(): + """ Test inserting a frame into an existing pipeline """ + w = wcs.WCS(pipe[:]) + original_result = w(1, 2) + mnew = models.Shift(1) & models.Shift(1) + new_frame = cf.Frame2D(name='new') + + # Insert at the beginning + w.insert_frame(new_frame, mnew, w.input_frame) + assert_allclose(w(0, 1), original_result) + + tr = w.get_transform('detector', w.output_frame) + assert_allclose(tr(1, 2), original_result) + + # Insert at the end + w = wcs.WCS(pipe[:]) + with pytest.raises(ValueError, match=r"New coordinate frame.*"): + w.insert_frame('not a frame', mnew, new_frame) + + w.insert_frame('icrs', mnew, new_frame) + assert_allclose([x - 1 for x in w(1, 2)], original_result) + + tr = w.get_transform('detector', 'icrs') + assert_allclose(tr(1, 2), original_result) + + # Force error by trying same operation + with pytest.raises(ValueError, match=r".*both frames.*"): + w.insert_frame('icrs', mnew, new_frame) + + def test_set_transform(): """ Test setting a transform between two frames in the pipeline.""" w = wcs.WCS(forward_transform=pipe[:]) @@ -374,25 +408,41 @@ def test_high_level_api(): """ Test WCS high level API. """ - output_frame = cf.CompositeFrame(frames=[icrs, spec]) - transform = m1 & models.Scale(1.5) - det = cf.CoordinateFrame(naxes=3, unit=(u.pix, u.pix, u.pix), - axes_order=(0, 1, 2), - axes_type=('length', 'length', 'length')) + output_frame = cf.CompositeFrame(frames=[icrs, spec, time]) + transform = m1 & models.Scale(1.5) & models.Scale(2) + det = cf.CoordinateFrame(naxes=4, unit=(u.pix, u.pix, u.pix, u.pix), + axes_order=(0, 1, 2, 3), + axes_type=('length', 'length', 'length', 'length')) w = wcs.WCS(forward_transform=transform, output_frame=output_frame, input_frame=det) + wrapped = wcsapi.HighLevelWCSWrapper(w) - r, d, lam = w(xv, yv, xv) - world_coord = w.pixel_to_world(xv, yv, xv) + r, d, lam, t = w(xv, yv, xv, xv) + world_coord = w.pixel_to_world(xv, yv, xv, xv) assert isinstance(world_coord[0], coord.SkyCoord) assert isinstance(world_coord[1], u.Quantity) + assert isinstance(world_coord[2], Time) assert_allclose(world_coord[0].data.lon.value, r) assert_allclose(world_coord[0].data.lat.value, d) assert_allclose(world_coord[1].value, lam) + assert_allclose((world_coord[2] - time.reference_frame).to(u.s).value, t) + + wrapped_world_coord = wrapped.pixel_to_world(xv, yv, xv, xv) + assert_allclose(wrapped_world_coord[0].data.lon.value, r) + assert_allclose(wrapped_world_coord[0].data.lat.value, d) + assert_allclose(wrapped_world_coord[1].value, lam) + assert_allclose((world_coord[2] - time.reference_frame).to(u.s).value, t) + + x1, y1, z1, k1 = w.world_to_pixel(*world_coord) + assert_allclose(x1, xv) + assert_allclose(y1, yv) + assert_allclose(z1, xv) + assert_allclose(k1, xv) - x1, y1, z1 = w.world_to_pixel(*world_coord) + x1, y1, z1, k1 = wrapped.world_to_pixel(*world_coord) assert_allclose(x1, xv) assert_allclose(y1, yv) assert_allclose(z1, xv) + assert_allclose(k1, xv) @pytest.mark.remote_data @@ -472,8 +522,7 @@ def test_footprint(self): def test_inverse(self): sky_coord = self.wcs(1, 2, with_units=True) - with pytest.raises(NotImplementedError): - self.wcs.invert(sky_coord) + assert np.allclose(self.wcs.invert(sky_coord), (1, 2)) def test_back_coordinates(self): sky_coord = self.wcs(1, 2, with_units=True) @@ -522,3 +571,229 @@ def test_to_fits_sip(): with pytest.raises(ValueError): miriwcs.bounding_box = None mirisip = miriwcs.to_fits_sip(bounding_box=None, max_inv_pix_error=0.1) + + +def test_to_fits_tab_no_bb(gwcs_3d_galactic_spectral): + # gWCS: + w = gwcs_3d_galactic_spectral + w.bounding_box = None + + # FITS WCS -TAB: + with pytest.raises(ValueError): + hdr, bt = w.to_fits_tab() + + +def test_to_fits_tab_cube(gwcs_3d_galactic_spectral): + # gWCS: + w = gwcs_3d_galactic_spectral + + # FITS WCS -TAB: + hdr, bt = w.to_fits_tab() + hdulist = fits.HDUList( + [fits.PrimaryHDU(np.ones(w.pixel_n_dim * (2, )), hdr), bt] + ) + fits_wcs = astwcs.WCS(hdulist[0].header, hdulist) + + hdr, bt = w.to_fits_tab(bounding_box=w.bounding_box) + hdulist = fits.HDUList( + [fits.PrimaryHDU(np.ones(w.pixel_n_dim * (2, )), hdr), bt] + ) + fits_wcs_user_bb = astwcs.WCS(hdulist[0].header, hdulist) + + # test points: + (xmin, xmax), (ymin, ymax), (zmin, zmax) = w.bounding_box + np.random.seed(1) + x = xmin + (xmax - xmin) * np.random.random(100) + y = ymin + (ymax - ymin) * np.random.random(100) + z = zmin + (zmax - zmin) * np.random.random(100) + + # test: + assert np.allclose(w(x, y, z), fits_wcs.wcs_pix2world(x, y, z, 0), + rtol=1e-6, atol=1e-7) + + assert np.allclose(w(x, y, z), fits_wcs_user_bb.wcs_pix2world(x, y, z, 0), + rtol=1e-6, atol=1e-7) + + +def test_to_fits_tab_miri_image(): + # gWCS: + af = asdf.open(get_pkg_data_filename('data/miriwcs.asdf')) + w = af.tree['wcs'] + + # FITS WCS -TAB: + hdr, bt = w.to_fits_tab(sampling=0.5) + hdulist = fits.HDUList( + [fits.PrimaryHDU(np.ones(w.pixel_n_dim * (2, )), hdr), bt] + ) + fits_wcs = astwcs.WCS(hdulist[0].header, hdulist) + + # test points: + (xmin, xmax), (ymin, ymax) = w.bounding_box + np.random.seed(1) + x = xmin + (xmax - xmin) * np.random.random(100) + y = ymin + (ymax - ymin) * np.random.random(100) + + # test: + assert np.allclose(w(x, y), fits_wcs.wcs_pix2world(x, y, 0), + rtol=1e-6, atol=1e-7) + + +def test_to_fits_tab_miri_lrs(): + af = asdf.open(get_pkg_data_filename('data/miri_lrs_wcs.asdf')) + w = af.tree['wcs'] + + # FITS WCS -TAB: + hdr, bt = w.to_fits_tab(sampling=0.25) + hdulist = fits.HDUList( + [fits.PrimaryHDU(np.ones(w.pixel_n_dim * (2, )), hdr), bt] + ) + fits_wcs = astwcs.WCS(hdulist[0].header, hdulist) + + # test points: + (xmin, xmax), (ymin, ymax) = w.bounding_box + np.random.seed(1) + x = xmin + (xmax - xmin) * np.random.random(100) + y = ymin + (ymax - ymin) * np.random.random(100) + + # test: + ref = np.array(w(x, y)) + tab = np.array(fits_wcs.wcs_pix2world(x, y, 0, 0)) + m = np.cumprod(np.isfinite(ref), dtype=np.bool_, axis=0) + assert np.allclose(ref[m], tab[m], rtol=5e-6, atol=5e-6, equal_nan=True) + + +def test_in_image(): + # create a 1-dim WCS: + w1 = wcs.WCS( + [ + (cf.SpectralFrame(name='input', axes_names=('x',), unit=(u.pix,)), models.Scale(2)), + (cf.SpectralFrame(name='output', axes_names=('x'), unit=(u.pix,)), None) + ] + ) + w1.bounding_box = (1, 5) + + assert np.isscalar(w1.in_image(4)) + assert w1.in_image(4) + assert not w1.in_image(14) + assert np.array_equal( + w1.in_image([[-1, 4, 11], [2, 3, 12]]), + [[False, True, False], [True, True, False]], + ) + + # create a 2-dim WCS: + w2 = wcs.WCS([(cf.Frame2D(name='input', axes_names=('x', 'y'), unit=(u.pix, u.pix)), + models.Scale(2) & models.Scale(1.5)), + (cf.Frame2D(name='output', axes_names=('x', 'y'), unit=(u.pix, u.pix)), + None)]) + w2.bounding_box = [(1, 100), (2, 20)] + + assert np.isscalar(w2.in_image(2, 6)) + assert not np.isscalar(w2.in_image([2], [6])) + 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]], + [[9, np.nan, 11, 15], [8, 9, 67, 98], [1, 1, np.nan, -10]] + ), + [[ True, False, True, True], [ True, True, False, False], [False, False, False, False]], + ) + + +def test_iter_inv(): + w = asdf.open(get_pkg_data_filename('data/nircamwcs.asdf')).tree['wcs'] + # remove analytic/user-supplied inverse: + w.pipeline[0].transform.inverse = None + w.bounding_box = None + + # test single point + assert np.allclose((1, 2), w.invert(*w(1, 2))) + assert np.allclose( + (np.nan, np.nan), + w.numerical_inverse(*w(np.nan, 2)), + equal_nan=True + ) + + # prepare to test a vector of points: + np.random.seed(10) + x, y = 2047 * np.random.random((2, 10000)) # "truth" + + # test adaptive: + xp, yp = w.invert( + *w(x, y), + adaptive=True, + detect_divergence=True, + quiet=False + ) + assert np.allclose((x, y), (xp, yp)) + + w = asdf.open(get_pkg_data_filename('data/nircamwcs.asdf')).tree['wcs'] + + # test single point + assert np.allclose((1, 2), w.numerical_inverse(*w(1, 2))) + assert np.allclose( + (np.nan, np.nan), + w.numerical_inverse(*w(np.nan, 2)), + equal_nan=True + ) + + # don't detect devergence + xp, yp = w.numerical_inverse( + *w(x, y), + adaptive=True, + detect_divergence=False, + quiet=False + ) + assert np.allclose((x, y), (xp, yp)) + + with pytest.raises(wcs.NoConvergence) as e: + w.numerical_inverse( + *w([1, 20, 200, 2000], [200, 1000, 2000, 5]), + adaptive=True, + detect_divergence=True, + maxiter=2, # force not reaching requested accuracy + quiet=False + ) + + xp, yp = e.value.best_solution.T + assert e.value.slow_conv.size == 4 + assert np.all(np.sort(e.value.slow_conv) == np.arange(4)) + + # test non-adaptive: + xp, yp = w.numerical_inverse( + *w(x, y, with_bounding_box=False), + adaptive=False, + detect_divergence=True, + quiet=False, + with_bounding_box=False + ) + assert np.allclose((x, y), (xp, yp)) + + # test non-adaptive: + x[0] = 3000 + y[0] = 10000 + xp, yp = w.numerical_inverse( + *w(x, y, with_bounding_box=False), + adaptive=False, + detect_divergence=True, + quiet=False, + with_bounding_box=False + ) + assert np.allclose((x, y), (xp, yp)) + + # test non-adaptive with non-recoverable divergence: + x[0] = 300000 + y[0] = 1000000 + with pytest.raises(wcs.NoConvergence) as e: + xp, yp = w.numerical_inverse( + *w(x, y, with_bounding_box=False), + adaptive=False, + detect_divergence=True, + quiet=False, + with_bounding_box=False + ) + assert np.allclose((x, y), (xp, yp)) + + xp, yp = e.value.best_solution.T + assert np.allclose((x[1:], y[1:]), (xp[1:], yp[1:])) + assert e.value.divergent[0] == 0 diff --git a/gwcs/utils.py b/gwcs/utils.py index 0d80f06a..706be3fc 100644 --- a/gwcs/utils.py +++ b/gwcs/utils.py @@ -11,6 +11,7 @@ from astropy.io import fits from astropy import coordinates as coords from astropy import units as u +from astropy.time import Time, TimeDelta # these ctype values do not include yzLN and yzLT pairs @@ -461,6 +462,8 @@ def isnumerical(val): isnum = False elif isinstance(val, u.Quantity): isnum = False + elif isinstance(val, (Time, TimeDelta)): + isnum = False elif (isinstance(val, np.ndarray) and not np.issubdtype(val.dtype, np.floating) and not np.issubdtype(val.dtype, np.integer)): diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 42fdf079..9297d839 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -4,10 +4,11 @@ import warnings import numpy as np import numpy.linalg as npla +from scipy import optimize from astropy.modeling.core import Model # , fix_inputs from astropy.modeling import utils as mutils from astropy.modeling.models import (Shift, Polynomial2D, Sky2Pix_TAN, - RotateCelestial2Native) + RotateCelestial2Native, Mapping) from astropy.modeling.fitting import LinearLSQFitter import astropy.io.fits as fits @@ -19,15 +20,61 @@ from gwcs import coordinate_frames as cf -HAS_FIX_INPUTS = True - try: from astropy.modeling.core import fix_inputs + HAS_FIX_INPUTS = True except ImportError: HAS_FIX_INPUTS = False -__all__ = ['WCS'] +__all__ = ['WCS', 'NoConvergence'] + +_ITER_INV_KWARGS = ['tolerance', 'maxiter', 'adaptive', 'detect_divergence', 'quiet'] + + +class NoConvergence(Exception): + """ + An error class used to report non-convergence and/or divergence + of numerical methods. It is used to report errors in the + iterative solution used by + the :py:meth:`~astropy.wcs.WCS.all_world2pix`. + + Attributes + ---------- + + best_solution : `numpy.ndarray` + Best solution achieved by the numerical method. + + accuracy : `numpy.ndarray` + Estimate of the accuracy of the ``best_solution``. + + niter : `int` + Number of iterations performed by the numerical method + to compute ``best_solution``. + + divergent : None, `numpy.ndarray` + Indices of the points in ``best_solution`` array + for which the solution appears to be divergent. If the + solution does not diverge, ``divergent`` will be set to `None`. + + slow_conv : None, `numpy.ndarray` + Indices of the solutions in ``best_solution`` array + for which the solution failed to converge within the + specified maximum number of iterations. If there are no + non-converging solutions (i.e., if the required accuracy + has been achieved for all input data points) + then ``slow_conv`` will be set to `None`. + + """ + def __init__(self, *args, best_solution=None, accuracy=None, niter=None, + divergent=None, slow_conv=None): + super().__init__(*args) + + self.best_solution = best_solution + self.accuracy = accuracy + self.niter = niter + self.divergent = divergent + self.slow_conv = slow_conv class WCS(GWCSAPIMixin): @@ -53,6 +100,7 @@ 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 = [] self._name = name @@ -283,42 +331,128 @@ def __call__(self, *args, **kwargs): return result + def in_image(self, *args, **kwargs): + """ + This method tests if one or more of the input world coordinates are + contained within forward transformation's image and that it maps to + the domain of definition of the forward transformation. + In practical terms, this function tests + that input world coordinate(s) can be converted to input frame and that + it is within the forward transformation's ``bounding_box`` when + defined. + + Parameters + ---------- + args : float, array like, `~astropy.coordinates.SkyCoord` or + `~astropy.units.Unit` coordinates to be inverted + + kwargs : dict + keyword arguments to be passed either to ``backward_transform`` + (when defined) or to the iterative invert method. + + Returns + ------- + result : bool, numpy.ndarray + A single boolean value or an array of boolean values with `True` + indicating that the WCS footprint contains the coordinate + and `False` if input is outside the footprint. + + """ + kwargs['with_bounding_box'] = True + kwargs['fill_value'] = np.nan + + 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 + + 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): """ - Invert coordinates. + Invert coordinates from output frame to input frame using analytical or + user-supplied inverse. When neither analytical nor user-supplied + inverses are defined, a numerical solution will be attempted using + :py:meth:`numerical_inverse`. - The analytical inverse of the forward transform is used, if available. - If not an iterative method is used. + .. note:: + Currently numerical inverse is implemented only for 2D imaging WCS. Parameters ---------- args : float, array like, `~astropy.coordinates.SkyCoord` or `~astropy.units.Unit` - coordinates to be inverted - kwargs : dict - keyword arguments to be passed to the iterative invert method. + Coordinates to be inverted. The number of arguments must be equal + to the number of world coordinates given by ``world_n_dim``. + 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 to ``fill_value``. + If `True` (default) values in the result which correspond to any + of the inputs being outside the bounding_box are set to + ``fill_value``. + fill_value : float, optional - Output value for inputs outside the bounding_box (default is np.nan). + 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.units.Quantity` object, by using the units of + the output cooridnate frame. Default is `False`. + + Other Parameters + ---------------- + kwargs : dict + Keyword arguments to be passed to :py:meth:`numerical_inverse` + (when defined) or to the iterative invert method. + + Returns + ------- + result : tuple + Returns a tuple of scalar or array values for each axis. + """ + with_units = kwargs.pop('with_units', False) + if not utils.isnumerical(args[0]): args = self.output_frame.coordinate_to_quantity(*args) if self.output_frame.naxes == 1: args = [args] - if not self.backward_transform.uses_quantity: + 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) - 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 try: - result = self.backward_transform(*args, **kwargs) + # 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) except (NotImplementedError, KeyError): - result = self._invert(*args, **kwargs) + result = self.numerical_inverse(*args, **kwargs, with_units=with_units) if with_units and self.input_frame: if self.input_frame.naxes == 1: @@ -328,11 +462,534 @@ def invert(self, *args, **kwargs): else: return result - def _invert(self, *args, **kwargs): + def numerical_inverse(self, *args, **kwargs): """ - Implement iterative inverse here. + Invert coordinates from output frame to input frame using numerical + inverse. + + .. note:: + Currently numerical inverse is implemented only for 2D imaging WCS. + + .. note:: + This method uses a combination of vectorized fixed-point + iterations algorithm and `scipy.optimize.root`. The later is used + for input coordinates for which vectorized algorithm diverges. + + Parameters + ---------- + args : float, array like, `~astropy.coordinates.SkyCoord` or `~astropy.units.Unit` + Coordinates to be inverted. The number of arguments must be equal + to the number of world coordinates given by ``world_n_dim``. + + 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 to + ``fill_value``. + + fill_value : float, optional + 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.units.Quantity` object, by using the units of + the output cooridnate frame. Default is `False`. + + tolerance : float, optional + *Absolute tolerance* of solution. Iteration terminates when the + iterative solver estimates that the "true solution" is + within this many pixels current estimate, more + specifically, when the correction to the solution found + during the previous iteration is smaller + (in the sense of the L2 norm) than ``tolerance``. + Default ``tolerance`` is 1.0e-5. + + maxiter : int, optional + Maximum number of iterations allowed to reach a solution. + Default is 50. + + quiet : bool, optional + Do not throw :py:class:`NoConvergence` exceptions when + the method does not converge to a solution with the + required accuracy within a specified number of maximum + iterations set by ``maxiter`` parameter. Instead, + simply return the found solution. Default is `True`. + + Other Parameters + ---------------- + adaptive : bool, optional + Specifies whether to adaptively select only points that + did not converge to a solution within the required + accuracy for the next iteration. Default (`True`) is recommended. + + .. note:: + The :py:meth:`numerical_inverse` uses a vectorized + implementation of the method of consecutive + approximations (see ``Notes`` section below) in which it + iterates over *all* input points *regardless* until + the required accuracy has been reached for *all* input + points. In some cases it may be possible that + *almost all* points have reached the required accuracy + but there are only a few of input data points for + which additional iterations may be needed (this + depends mostly on the characteristics of the geometric + distortions for a given instrument). In this situation + it may be advantageous to set ``adaptive`` = `True` in + which case :py:meth:`numerical_inverse` will continue + iterating *only* over the points that have not yet + converged to the required accuracy. + + .. note:: + When ``detect_divergence`` is `True`, + :py:meth:`numerical_inverse` will automatically switch + to the adaptive algorithm once divergence has been + detected. + + detect_divergence : bool, optional + Specifies whether to perform a more detailed analysis + of the convergence to a solution. Normally + :py:meth:`numerical_inverse` may not achieve the required + accuracy if either the ``tolerance`` or ``maxiter`` arguments + are too low. However, it may happen that for some + geometric distortions the conditions of convergence for + the the method of consecutive approximations used by + :py:meth:`numerical_inverse` may not be satisfied, in which + case consecutive approximations to the solution will + diverge regardless of the ``tolerance`` or ``maxiter`` + settings. + + When ``detect_divergence`` is `False`, these divergent + points will be detected as not having achieved the + required accuracy (without further details). In addition, + if ``adaptive`` is `False` then the algorithm will not + know that the solution (for specific points) is diverging + and will continue iterating and trying to "improve" + diverging solutions. This may result in ``NaN`` or + ``Inf`` values in the return results (in addition to a + performance penalties). Even when ``detect_divergence`` + is `False`, :py:meth:`numerical_inverse`, at the end of the + iterative process, will identify invalid results + (``NaN`` or ``Inf``) as "diverging" solutions and will + raise :py:class:`NoConvergence` unless the ``quiet`` + parameter is set to `True`. + + When ``detect_divergence`` is `True` (default), + :py:meth:`numerical_inverse` will detect points for which + current correction to the coordinates is larger than + the correction applied during the previous iteration + **if** the requested accuracy **has not yet been + achieved**. In this case, if ``adaptive`` is `True`, + these points will be excluded from further iterations and + if ``adaptive`` is `False`, :py:meth:`numerical_inverse` will + automatically switch to the adaptive algorithm. Thus, the + reported divergent solution will be the latest converging + solution computed immediately *before* divergence + has been detected. + + .. note:: + When accuracy has been achieved, small increases in + current corrections may be possible due to rounding + errors (when ``adaptive`` is `False`) and such + increases will be ignored. + + .. note:: + Based on our testing using JWST NIRCAM images, setting + ``detect_divergence`` to `True` will incur about 5-10% + performance penalty with the larger penalty + corresponding to ``adaptive`` set to `True`. + Because the benefits of enabling this + feature outweigh the small performance penalty, + especially when ``adaptive`` = `False`, it is + recommended to set ``detect_divergence`` to `True`, + unless extensive testing of the distortion models for + images from specific instruments show a good stability + of the numerical method for a wide range of + coordinates (even outside the image itself). + + .. note:: + Indices of the diverging inverse solutions will be + reported in the ``divergent`` attribute of the + raised :py:class:`NoConvergence` exception object. + + Returns + ------- + result : tuple + Returns a tuple of scalar or array values for each axis. + + Raises + ------ + NoConvergence + The iterative method did not converge to a + solution to the required accuracy within a specified + number of maximum iterations set by the ``maxiter`` + parameter. To turn off this exception, set ``quiet`` to + `True`. Indices of the points for which the requested + accuracy was not achieved (if any) will be listed in the + ``slow_conv`` attribute of the + raised :py:class:`NoConvergence` exception object. + + See :py:class:`NoConvergence` documentation for + more details. + + NotImplementedError + Numerical inverse has not been implemented for this WCS. + + ValueError + Invalid argument values. + """ - raise NotImplementedError + tolerance = kwargs.get('tolerance', 1e-5) + maxiter = kwargs.get('maxiter', 50) + adaptive = kwargs.get('adaptive', True) + detect_divergence = kwargs.get('detect_divergence', True) + quiet = kwargs.get('quiet', True) + with_bounding_box = kwargs.get('with_bounding_box', True) + fill_value = kwargs.get('fill_value', np.nan) + with_units = kwargs.pop('with_units', False) + + 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 + + if nargs != self.world_n_dim: + raise ValueError("Number of input coordinates is different from " + "the number of defined world coordinates in the " + f"WCS ({self.world_n_dim:d})") + + if self.world_n_dim != self.pixel_n_dim: + raise NotImplementedError( + "Support for iterative inverse for transformations with " + "different number of inputs and outputs was not implemented." + ) + + # initial guess: + if nargs == 2 and self._approx_inverse is None: + self._calc_approx_inv(max_inv_pix_error=5, inv_degree=None) + + if self._approx_inverse is None: + if self.bounding_box is None: + x0 = np.ones(self.pixel_n_dim) + else: + x0 = np.mean(self.bounding_box, axis=-1) + + if arg_dim == 0: + argsi = args + + if nargs == 2 and self._approx_inverse is not None: + x0 = self._approx_inverse(*argsi) + if not np.all(np.isfinite(x0)): + return [np.array(np.nan) for _ in range(nargs)] + + result = tuple(self._vectorized_fixed_point( + x0, argsi, + tolerance=tolerance, + maxiter=maxiter, + adaptive=adaptive, + detect_divergence=detect_divergence, + quiet=quiet, + with_bounding_box=with_bounding_box, + fill_value=fill_value + ).T.ravel().tolist()) + + else: + arg_shape = args_shape[1:] + nelem = np.prod(arg_shape) + + args = np.reshape(args, (nargs, nelem)) + + if self._approx_inverse is None: + x0 = np.full((nelem, nargs), x0) + else: + x0 = np.array(self._approx_inverse(*args)).T + + result = self._vectorized_fixed_point( + x0, args.T, + tolerance=tolerance, + maxiter=maxiter, + adaptive=adaptive, + detect_divergence=detect_divergence, + quiet=quiet, + with_bounding_box=with_bounding_box, + fill_value=fill_value + ).T + + 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 + + def _vectorized_fixed_point(self, pix0, world, tolerance, maxiter, + adaptive, detect_divergence, quiet, + with_bounding_box, fill_value): + # ############################################################ + # # INITIALIZE ITERATIVE PROCESS: ## + # ############################################################ + + # make a copy of the initial approximation + pix0 = np.atleast_2d(np.array(pix0)) # 0-order solution + pix = np.array(pix0) + + world0 = np.atleast_2d(np.array(world)) + world = np.array(world0) + + # estimate pixel scale using approximate algorithm + # from https://trs.jpl.nasa.gov/handle/2014/40409 + if self.bounding_box is None: + crpix = np.ones(self.pixel_n_dim) + else: + crpix = np.mean(self.bounding_box, axis=-1) + + l1, phi1 = np.deg2rad(self.__call__(*(crpix - 0.5))) + l2, phi2 = np.deg2rad(self.__call__(*(crpix + [-0.5, 0.5]))) + l3, phi3 = np.deg2rad(self.__call__(*(crpix + 0.5))) + l4, phi4 = np.deg2rad(self.__call__(*(crpix + [0.5, -0.5]))) + area = np.abs(0.5 * ((l4 - l2) * np.sin(phi1) + + (l1 - l3) * np.sin(phi2) + + (l2 - l4) * np.sin(phi3) + + (l3 - l2) * np.sin(phi4))) + inv_pscale = 1 / np.rad2deg(np.sqrt(area)) + + # form equation: + def f(x): + w = np.array(self.__call__(*(x.T), with_bounding_box=False)).T + dw = np.mod(np.subtract(w, world) - 180.0, 360.0) - 180.0 + return np.add(inv_pscale * dw, x) + + def froot(x): + return np.mod(np.subtract(self.__call__(*x, with_bounding_box=False), worldi) - 180.0, 360.0) - 180.0 + + # compute correction: + def correction(pix): + p1 = f(pix) + p2 = f(p1) + d = p2 - 2.0 * p1 + pix + idx = np.where(d != 0) + corr = pix - p2 + corr[idx] = np.square(p1[idx] - pix[idx]) / d[idx] + return corr + + # initial iteration: + dpix = correction(pix) + + # Update initial solution: + pix -= dpix + + # Norm (L2) squared of the correction: + dn = np.sum(dpix * dpix, axis=1) + dnprev = dn.copy() # if adaptive else dn + tol2 = tolerance**2 + + # Prepare for iterative process + k = 1 + ind = None + inddiv = None + + # Turn off numpy runtime warnings for 'invalid' and 'over': + old_invalid = np.geterr()['invalid'] + old_over = np.geterr()['over'] + np.seterr(invalid='ignore', over='ignore') + + # ############################################################ + # # NON-ADAPTIVE ITERATIONS: ## + # ############################################################ + if not adaptive: + # Fixed-point iterations: + while (np.nanmax(dn) >= tol2 and k < maxiter): + # Find correction to the previous solution: + dpix = correction(pix) + + # Compute norm (L2) squared of the correction: + dn = np.sum(dpix * dpix, axis=1) + + # Check for divergence (we do this in two stages + # to optimize performance for the most common + # scenario when successive approximations converge): + + if detect_divergence: + divergent = (dn >= dnprev) + if np.any(divergent): + # Find solutions that have not yet converged: + slowconv = (dn >= tol2) + inddiv, = np.where(divergent & slowconv) + + if inddiv.shape[0] > 0: + # Update indices of elements that + # still need correction: + conv = (dn < dnprev) + iconv = np.where(conv) + + # Apply correction: + dpixgood = dpix[iconv] + pix[iconv] -= dpixgood + dpix[iconv] = dpixgood + + # For the next iteration choose + # non-divergent points that have not yet + # converged to the requested accuracy: + ind, = np.where(slowconv & conv) + world = world[ind] + dnprev[ind] = dn[ind] + k += 1 + + # Switch to adaptive iterations: + adaptive = True + break + + # Save current correction magnitudes for later: + dnprev = dn + + # Apply correction: + pix -= dpix + k += 1 + + # ############################################################ + # # ADAPTIVE ITERATIONS: ## + # ############################################################ + if adaptive: + if ind is None: + ind, = np.where(np.isfinite(pix).all(axis=1)) + world = world[ind] + + # "Adaptive" fixed-point iterations: + while (ind.shape[0] > 0 and k < maxiter): + # Find correction to the previous solution: + dpixnew = correction(pix[ind]) + + # Compute norm (L2) of the correction: + dnnew = np.sum(np.square(dpixnew), axis=1) + + # Bookkeeping of corrections: + dnprev[ind] = dn[ind].copy() + dn[ind] = dnnew + + if detect_divergence: + # Find indices of pixels that are converging: + conv = np.logical_or(dnnew < dnprev[ind], dnnew < tol2) + if not np.all(conv): + conv = np.ones_like(dnnew, dtype=np.bool) + iconv = np.where(conv) + iiconv = ind[iconv] + + # Apply correction: + dpixgood = dpixnew[iconv] + pix[iiconv] -= dpixgood + dpix[iiconv] = dpixgood + + # Find indices of solutions that have not yet + # converged to the requested accuracy + # AND that do not diverge: + subind, = np.where((dnnew >= tol2) & conv) + + else: + # Apply correction: + pix[ind] -= dpixnew + dpix[ind] = dpixnew + + # Find indices of solutions that have not yet + # converged to the requested accuracy: + subind, = np.where(dnnew >= tol2) + + # Choose solutions that need more iterations: + ind = ind[subind] + world = world[subind] + + k += 1 + + # ############################################################ + # # FINAL DETECTION OF INVALID, DIVERGING, ## + # # AND FAILED-TO-CONVERGE POINTS ## + # ############################################################ + # Identify diverging and/or invalid points: + invalid = ((~np.all(np.isfinite(pix), axis=1)) & + (np.all(np.isfinite(world0), axis=1))) + + # When detect_divergence is False, dnprev is outdated + # (it is the norm of the very first correction). + # Still better than nothing... + inddiv, = np.where(((dn >= tol2) & (dn >= dnprev)) | invalid) + if inddiv.shape[0] == 0: + inddiv = None + + # If there are divergent points, attempt to find a solution using + # scipy's 'hybr' method: + if detect_divergence and inddiv is not None and inddiv.size: + bad = [] + for idx in inddiv: + worldi = world0[idx] + result = optimize.root( + froot, + pix0[idx], + method='hybr', + tol=tolerance / (np.linalg.norm(pix0[idx]) + 1), + options={'maxfev': 2 * maxiter} + ) + + if result['success']: + pix[idx, :] = result['x'] + invalid[idx] = False + else: + bad.append(idx) + + if bad: + inddiv = np.array(bad, dtype=np.int) + else: + inddiv = None + + # Identify points that did not converge within 'maxiter' + # iterations: + if k >= maxiter: + ind, = np.where((dn >= tol2) & (dn < dnprev) & (~invalid)) + if ind.shape[0] == 0: + ind = None + else: + ind = None + + # Restore previous numpy error settings: + np.seterr(invalid=old_invalid, over=old_over) + + # ############################################################ + # # RAISE EXCEPTION IF DIVERGING OR TOO SLOWLY CONVERGING ## + # # DATA POINTS HAVE BEEN DETECTED: ## + # ############################################################ + if (ind is not None or inddiv is not None) and not quiet: + if inddiv is None: + raise NoConvergence( + "'WCS.invert' failed to " + "converge to the requested accuracy after {:d} " + "iterations.".format(k), best_solution=pix, + accuracy=np.abs(dpix), niter=k, + slow_conv=ind, divergent=None) + else: + raise NoConvergence( + "'WCS.invert' failed to " + "converge to the requested accuracy.\n" + "After {:d} iterations, the solution is diverging " + "at least for one input point." + .format(k), best_solution=pix, + accuracy=np.abs(dpix), niter=k, + slow_conv=ind, divergent=inddiv) + + if with_bounding_box and self.bounding_box is not None: + # find points outside the bounding box and replace their values + # with fill_value + valid = np.logical_not(invalid) + in_bb = np.ones_like(invalid, dtype=np.bool_) + + for c, (x1, x2) in zip(pix[valid].T, self.bounding_box): + in_bb[valid] &= (c >= x1) & (c <= x2) + pix[np.logical_not(in_bb)] = fill_value + + return pix def transform(self, from_frame, to_frame, *args, **kwargs): """ @@ -424,6 +1081,62 @@ def insert_transform(self, frame, transform, after=False): current_transform = self._pipeline[frame_ind].transform self._pipeline[frame_ind].transform = transform | current_transform + def insert_frame(self, input_frame, transform, output_frame): + """ + Insert a new frame into an existing pipeline. This frame must be + anchored to a frame already in the pipeline by a transform. This + existing frame is identified solely by its name, although an entire + `~gwcs.coordinate_frame.CoordinateFrame` can be passed (e.g., the + `input_frame` or `output_frame` attribute). This frame is never + modified. + + Parameters + ---------- + input_frame : str or `~gwcs.coordinate_frame.CoordinateFrame` + Coordinate frame at start of new transform + transform : `~astropy.modeling.Model` + New transform to be inserted in the pipeline + output_frame: str or `~gwcs.coordinate_frame.CoordinateFrame` + Coordinate frame at end of new transform + """ + input_name, input_frame_obj = self._get_frame_name(input_frame) + output_name, output_frame_obj = self._get_frame_name(output_frame) + try: + input_index = self._get_frame_index(input_frame) + except ValueError: + input_index = None + if input_frame_obj is None: + raise ValueError(f"New coordinate frame {input_name} must " + "be defined") + try: + output_index = self._get_frame_index(output_frame) + except ValueError: + output_index = None + if output_frame_obj is None: + raise ValueError(f"New coordinate frame {output_name} must " + "be defined") + + new_frames = [input_index, output_index].count(None) + if new_frames == 0: + raise ValueError("Could not insert frame as both frames " + f"{input_name} and {output_name} already exist") + elif new_frames == 2: + raise ValueError("Could not insert frame as neither frame " + f"{input_name} nor {output_name} exists") + + if input_index is None: + self._pipeline = (self._pipeline[:output_index] + + [Step(input_frame_obj, transform)] + + self._pipeline[output_index:]) + super(WCS, self).__setattr__(input_name, input_frame_obj) + else: + split_step = self._pipeline[input_index] + self._pipeline = (self._pipeline[:input_index] + + [Step(split_step.frame, transform), + Step(output_frame_obj, split_step.transform)] + + self._pipeline[input_index + 1:]) + super(WCS, self).__setattr__(output_name, output_frame_obj) + @property def unit(self): """The unit of the coordinates in the output coordinate system.""" @@ -440,7 +1153,6 @@ def unit(self): def output_frame(self): """Return the output coordinate frame.""" if self._pipeline: - #frame = self._pipeline[-1][0] frame = self._pipeline[-1].frame if not isinstance(frame, str): frame = frame.name @@ -452,7 +1164,6 @@ def output_frame(self): def input_frame(self): """Return the input coordinate frame.""" if self._pipeline: - #frame = self._pipeline[0][0] frame = self._pipeline[0].frame if not isinstance(frame, str): frame = frame.name @@ -806,6 +1517,201 @@ def to_fits_sip(self, bounding_box=None, max_pix_error=0.25, degree=None, hdr['cd2_2'] = cdmat[1][1] return hdr + def to_fits_tab(self, bounding_box=None, bin_ext_name='WCS-TABLE', + coord_col_name='coordinates', sampling=1): + """ + Construct a FITS WCS ``-TAB``-based approximation to the WCS + in the form of a FITS header and a binary table extension. For the + description of the FITS WCS ``-TAB`` convention, see + "Representations of spectral coordinates in FITS" in + `Greisen, E. W. et al. A&A 446 (2) 747-771 (2006) + `_ . + + Parameters + ---------- + bounding_box : tuple, optional + Specifies the range of acceptable values for each input axis. + The order of the axes is + `~gwcs.coordinate_frames.CoordinateFrame.axes_order`. + For two image axes ``bounding_box`` is of the form + ``((xmin, xmax), (ymin, ymax))``. + + bin_ext_name : str, optional + Extension name for the `~astropy.io.fits.BinTableHDU` extension. + + coord_col_name : str, optional + Field name of the coordinate array in the structured array + stored in `~astropy.io.fits.BinTableHDU` data. This corresponds to + ``TTYPEi`` field in the FITS header of the binary table extension. + + sampling : float, tuple, optional + The target "density" of grid nodes per pixel to be used when + creating the coordinate array for the ``-TAB`` FITS WCS convention. + It is equal to ``1/step`` where ``step`` is the distance between + grid nodes in pixels. ``sampling`` can be specified as a single + number to be used for all axes or as a `tuple` of numbers + that specify the sampling for each image axis. + + Returns + ------- + hdr : `~astropy.io.fits.Header` + Header with WCS-TAB information associated (to be used) with image + data. + + bin_table : `~astropy.io.fits.BinTableHDU` + Binary table extension containing the coordinate array. + + Raises + ------ + ValueError + When ``bounding_box`` is not defined either through the input + ``bounding_box`` parameter or this object's ``bounding_box`` + property. + + ValueError + When ``sampling`` is a `tuple` of length larger than 1 that + does not match the number of image axes. + + RuntimeError + If the number of image axes (`~gwcs.WCS.pixel_n_dim`) is larger + than the number of world axes (`~gwcs.WCS.world_n_dim`). + + """ + if bounding_box is None: + if self.bounding_box is None: + raise ValueError( + "Need a valid bounding_box to compute the footprint." + ) + bounding_box = self.bounding_box + + else: + # validate user-supplied bounding box: + frames = self.available_frames + transform_0 = self.get_transform(frames[0], frames[1]) + mutils._BoundingBox.validate(transform_0, bounding_box) + + if self.pixel_n_dim > self.world_n_dim: + raise RuntimeError( + "The case when the number of input axes is larger than the " + "number of output axes is not supported." + ) + + try: + sampling = np.broadcast_to(sampling, (self.pixel_n_dim, )) + except ValueError: + raise ValueError("Number of sampling values either must be 1 " + "or it must match the number of pixel axes.") + + # 1D grid coordinates: + gcrds = [] + cdelt = [] + for (xmin, xmax), s in zip(bounding_box, sampling): + npix = max(2, 1 + int(np.ceil(abs((xmax - xmin) / s)))) + gcrds.append(np.linspace(xmin, xmax, npix)) + cdelt.append((npix - 1) / (xmax - xmin) if xmin != xmax else 1) + + # n-dim coordinate arrays: + coord = np.stack( + self(*np.meshgrid(*gcrds[::-1], indexing='ij')[::-1]), + axis=-1 + ) + + # create header with WCS info: + hdr = fits.Header() + + for k in range(self.world_n_dim): + k1 = k + 1 + ct = cf.get_ctype_from_ucd(self.world_axis_physical_types[k]) + if len(ct) > 4: + raise ValueError("Axis type name too long.") + + hdr['CTYPE{:d}'.format(k1)] = ct + (4 - len(ct)) * '-' + '-TAB' + hdr['CUNIT{:d}'.format(k1)] = self.world_axis_units[k] + hdr['PS{:d}_0'.format(k1)] = bin_ext_name + hdr['PS{:d}_1'.format(k1)] = coord_col_name + hdr['PV{:d}_3'.format(k1)] = k1 + hdr['CRVAL{:d}'.format(k1)] = 1 + + if k < self.pixel_n_dim: + hdr['CRPIX{:d}'.format(k1)] = gcrds[k][0] + 1 + hdr['PC{0:d}_{0:d}'.format(k1)] = 1.0 + hdr['CDELT{:d}'.format(k1)] = cdelt[k] + else: + hdr['CRPIX{:d}'.format(k1)] = 1 + coord = coord[None, :] + + # structured array (data) for binary table HDU: + arr = np.array( + [(coord, )], + dtype=[ + (coord_col_name, np.float64, coord.shape), + ] + ) + + # create binary table HDU: + bin_tab = fits.BinTableHDU(arr) + bin_tab.header['EXTNAME'] = bin_ext_name + + return hdr, bin_tab + + def _calc_approx_inv(self, max_inv_pix_error=5, inv_degree=None, npoints=16): + """ + Compute polynomial fit for the inverse transformation to be used as + initial aproximation/guess for the iterative solution. + """ + self._approx_inverse = None + + try: + # try to use analytic inverse if available: + self._approx_inverse = functools.partial(self.backward_transform, + with_bounding_box=False) + return + except (NotImplementedError, KeyError): + pass + + if not isinstance(self.output_frame, cf.CelestialFrame): + # The _calc_approx_inv method only works with celestial frame transforms + return + + # Determine reference points. + if self.bounding_box is None: + # A bounding_box is needed to proceed. + return + + (xmin, xmax), (ymin, ymax) = self.bounding_box + crpix1 = (xmax - xmin) // 2 + crpix2 = (ymax - ymin) // 2 + crval1, crval2 = self.forward_transform(crpix1, crpix2) + + # Rotate to native system and deproject. Set center of the projection + # transformation to the middle of the bounding box ("image") in order + # to minimize projection effects across the entire image, + # thus the initial shift. + ntransform = ((Shift(crpix1) & Shift(crpix2)) | self.forward_transform + | RotateCelestial2Native(crval1, crval2, 180) + | Sky2Pix_TAN()) + + u, v = _make_sampling_grid(npoints, self.bounding_box) + undist_x, undist_y = ntransform(u, v) + + # Double sampling to check if sampling is sufficient. + ud, vd = _make_sampling_grid(2 * npoints, self.bounding_box) + undist_xd, undist_yd = ntransform(ud, vd) + + fit_inv_poly_u, fit_inv_poly_v, max_inv_resid = _fit_2D_poly( + ntransform, + npoints, None, + max_inv_pix_error, + undist_x, undist_y, u, v, + undist_xd, undist_yd, ud, vd, + verbose=True + ) + + self._approx_inverse = (RotateCelestial2Native(crval1, crval2, 180) | + Sky2Pix_TAN() | Mapping((0, 1, 0, 1)) | + (fit_inv_poly_u & fit_inv_poly_v) | + (Shift(crpix1) & Shift(crpix2))) + def _fit_2D_poly(ntransform, npoints, degree, max_error, xin, yin, xout, yout, diff --git a/setup.py b/setup.py index 8272e92c..b578bc97 100755 --- a/setup.py +++ b/setup.py @@ -67,8 +67,9 @@ def get_package_data(): setup_requires=['setuptools_scm'], description=DESCRIPTION, install_requires=[ - 'astropy', + 'astropy>=4.1', 'numpy', + 'scipy', 'asdf'], packages=find_packages(), extras_require={