diff --git a/CHANGES.rst b/CHANGES.rst index fc4bafc4..4c155b31 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,20 @@ +0.13.0 (2020-03-26) +------------------- +New Features +^^^^^^^^^^^^ + +- Added two new transforms - ``SphericalToCartesian`` and + ``CartesianToSpherical``. [#275, #284, #285] + +- Added ``to_fits_sip`` method to generate FITS header with SIP keywords [#286] + +- Added ``get_ctype_from_ucd`` function. [#288] + +Bug Fixes +^^^^^^^^^ + +- Fixed an off by one issue in ``utils.make_fitswcs_transform``. [#290] + 0.12.0 (2019-12-24) ------------------- New Features diff --git a/gwcs/coordinate_frames.py b/gwcs/coordinate_frames.py index a958abcc..6a955024 100644 --- a/gwcs/coordinate_frames.py +++ b/gwcs/coordinate_frames.py @@ -18,6 +18,34 @@ 'CoordinateFrame', 'TemporalFrame'] +UCD1_TO_CTYPE = { + 'pos.eq.ra': 'RA', + 'pos.eq.dec': 'DEC', + 'pos.galactic.lon': 'GLON', + 'pos.galactic.lat': 'GLAT', + 'pos.ecliptic.lon': 'ELON', + 'pos.ecliptic.lat': 'ELAT', + 'pos.bodyrc.lon': 'TLON', + 'pos.bodyrc.lat': 'TLAT', + 'custom:pos.helioprojective.lat': 'HPLT', + 'custom:pos.helioprojective.lon': 'HPLN', + 'custom:pos.heliographic.stonyhurst.lon': 'HGLN', + 'custom:pos.heliographic.stonyhurst.lat': 'HGLT', + 'custom:pos.heliographic.carrington.lon': 'CRLN', + 'custom:pos.heliographic.carrington.lat': 'CRLT', + 'em.freq': 'FREQ', + 'em.energy': 'ENER', + 'em.wavenumber': 'WAVN', + 'em.wl': 'WAVE', + 'spect.dopplerVeloc.radio': 'VRAD', + 'spect.dopplerVeloc.opt': 'VOPT', + 'src.redshift': 'ZOPT', + 'spect.dopplerVeloc': 'VELO', + 'custom:spect.doplerVeloc.beta': 'BETA', + 'time': 'TIME', + } + + STANDARD_REFERENCE_FRAMES = [frame.upper() for frame in coord.builtin_frames.__all__] STANDARD_REFERENCE_POSITION = ["GEOCENTER", "BARYCENTER", "HELIOCENTER", @@ -25,6 +53,23 @@ "GALACTIC_CENTER", "LOCAL_GROUP_CENTER"] +def get_ctype_from_ucd(ucd): + """ + Return the FITS ``CTYPE`` corresponding to a UCD1 value. + + Parameters + ---------- + ucd : str + UCD string, for example one of ```WCS.world_axis_physical_types``. + + Returns + ------- + CTYPE : str + The corresponding FITS ``CTYPE`` value or an empty string. + """ + return UCD1_TO_CTYPE.get(ucd, "") + + class CoordinateFrame: """ Base class for Coordinate Frames. @@ -475,7 +520,6 @@ def _convert_to_time(self, dt, *, unit, **kwargs): else: return time.Time(dt, **kwargs) - def coordinate_to_quantity(self, *coords): if isinstance(coords[0], time.Time): ref_value = self.reference_frame.value @@ -629,11 +673,11 @@ def from_index(cls, indexes): indexes = np.asanyarray(indexes, dtype=int) out = np.empty_like(indexes, dtype=object) - out[nans] = np.nan - for profile, index in cls.profiles.items(): out[indexes == index] = profile + out[nans] = np.nan + if out.size == 1 and not nans: return StokesProfile(out.item()) elif nans.all(): diff --git a/gwcs/geometry.py b/gwcs/geometry.py index 59071383..bcb9f79b 100644 --- a/gwcs/geometry.py +++ b/gwcs/geometry.py @@ -3,11 +3,13 @@ Spectroscopy related models. """ +import numbers import numpy as np from astropy.modeling.core import Model +from astropy import units as u - -__all__ = ['ToDirectionCosines', 'FromDirectionCosines'] +__all__ = ['ToDirectionCosines', 'FromDirectionCosines', + 'SphericalToCartesian', 'CartesianToSpherical'] class ToDirectionCosines(Model): @@ -55,3 +57,152 @@ def evaluate(self, cosa, cosb, cosc, length): def inverse(self): return ToDirectionCosines() + + +class SphericalToCartesian(Model): + """ + Convert spherical coordinates on a unit sphere to cartesian coordinates. + Spherical coordinates when not provided as ``Quantity`` are assumed + to be in degrees with ``lon`` being the *longitude (or azimuthal) angle* + ``[0, 360)`` (or ``[-180, 180)``) and angle ``lat`` is the *latitude* + (or *elevation angle*) in range ``[-90, 90]``. + + """ + _separable = False + + _input_units_allow_dimensionless = True + + n_inputs = 2 + n_outputs = 3 + + def __init__(self, wrap_lon_at=360, **kwargs): + """ + Parameters + ---------- + wrap_lon_at : {360, 180}, optional + An **integer number** that specifies the range of the longitude + (azimuthal) angle. When ``wrap_lon_at`` is 180, the longitude angle + will have a range of ``[-180, 180)`` and when ``wrap_lon_at`` + is 360 (default), the longitude angle will have a range of + ``[0, 360)``. + + """ + super().__init__(**kwargs) + self.inputs = ('lon', 'lat') + self.outputs = ('x', 'y', 'z') + self.wrap_lon_at = wrap_lon_at + + @property + def wrap_lon_at(self): + """ An **integer number** that specifies the range of the longitude + (azimuthal) angle. + + Allowed values are 180 and 360. When ``wrap_lon_at`` + is 180, the longitude angle will have a range of ``[-180, 180)`` and + when ``wrap_lon_at`` is 360 (default), the longitude angle will have a + range of ``[0, 360)``. + + """ + return self._wrap_lon_at + + @wrap_lon_at.setter + def wrap_lon_at(self, wrap_angle): + if not (isinstance(wrap_angle, numbers.Integral) and wrap_angle in [180, 360]): + raise ValueError("'wrap_lon_at' must be an integer number: 180 or 360") + self._wrap_lon_at = wrap_angle + + def evaluate(self, lon, lat): + if isinstance(lon, u.Quantity) != isinstance(lat, u.Quantity): + raise TypeError("All arguments must be of the same type " + "(i.e., quantity or not).") + + lon = np.deg2rad(lon) + lat = np.deg2rad(lat) + + cs = np.cos(lat) + x = cs * np.cos(lon) + y = cs * np.sin(lon) + z = np.sin(lat) + + return x, y, z + + def inverse(self): + return CartesianToSpherical(wrap_lon_at=self._wrap_lon_at) + + @property + def input_units(self): + return {'lon': u.deg, 'lat': u.deg} + + +class CartesianToSpherical(Model): + """ + Convert cartesian coordinates to spherical coordinates on a unit sphere. + Output spherical coordinates are in degrees. When input cartesian + coordinates are quantities (``Quantity`` objects), output angles + will also be quantities in degrees. Angle ``lon`` is the *longitude* + (or *azimuthal angle*) in range ``[0, 360)`` (or ``[-180, 180)``) and + angle ``lat`` is the *latitude* (or *elevation angle*) in the + range ``[-90, 90]``. + + """ + _separable = False + + _input_units_allow_dimensionless = True + + n_inputs = 3 + n_outputs = 2 + + def __init__(self, wrap_lon_at=360, **kwargs): + """ + Parameters + ---------- + wrap_lon_at : {360, 180}, optional + An **integer number** that specifies the range of the longitude + (azimuthal) angle. When ``wrap_lon_at`` is 180, the longitude angle + will have a range of ``[-180, 180)`` and when ``wrap_lon_at`` + is 360 (default), the longitude angle will have a range of + ``[0, 360)``. + + """ + super().__init__(**kwargs) + self.inputs = ('x', 'y', 'z') + self.outputs = ('lon', 'lat') + self.wrap_lon_at = wrap_lon_at + + @property + def wrap_lon_at(self): + """ An **integer number** that specifies the range of the longitude + (azimuthal) angle. + + Allowed values are 180 and 360. When ``wrap_lon_at`` + is 180, the longitude angle will have a range of ``[-180, 180)`` and + when ``wrap_lon_at`` is 360 (default), the longitude angle will have a + range of ``[0, 360)``. + + """ + return self._wrap_lon_at + + @wrap_lon_at.setter + def wrap_lon_at(self, wrap_angle): + if not (isinstance(wrap_angle, numbers.Integral) and wrap_angle in [180, 360]): + raise ValueError("'wrap_lon_at' must be an integer number: 180 or 360") + self._wrap_lon_at = wrap_angle + + def evaluate(self, x, y, z): + nquant = [isinstance(i, u.Quantity) for i in (x, y, z)].count(True) + if nquant in [1, 2]: + raise TypeError("All arguments must be of the same type " + "(i.e., quantity or not).") + + h = np.hypot(x, y) + lat = np.rad2deg(np.arctan2(z, h)) + lon = np.rad2deg(np.arctan2(y, x)) + lon[h == 0] *= 0 + + if self._wrap_lon_at != 180: + lon = np.mod(lon, 360.0 * u.deg if nquant else 360.0) + + return lon, lat + + def inverse(self): + return SphericalToCartesian(wrap_lon_at=self._wrap_lon_at) diff --git a/gwcs/schemas/stsci.edu/gwcs/spherical_cartesian-1.0.0.yaml b/gwcs/schemas/stsci.edu/gwcs/spherical_cartesian-1.0.0.yaml new file mode 100644 index 00000000..61f959eb --- /dev/null +++ b/gwcs/schemas/stsci.edu/gwcs/spherical_cartesian-1.0.0.yaml @@ -0,0 +1,44 @@ +%YAML 1.1 +--- +$schema: "http://stsci.edu/schemas/yaml-schema/draft-01" +id: "http://stsci.edu/schemas/gwcs/spherical_cartesian-1.0.0" +tag: "tag:stsci.edu:gwcs/spherical_cartesian-1.0.0" + +title: > + Convert coordinates between spherical and Cartesian coordinates. + +description: | + This schema is for transforms which convert between spherical coordinates + (on the unit sphere) and Cartesian coordinates. + +examples: + - + - Convert spherical coordinates to Cartesian coordinates. + + - | + + ! + transform_type: spherical_to_cartesian + + - + - Convert Cartesian coordinates to spherical coordinates. + + - | + ! + transform_type: cartesian_to_spherical + +allOf: + - $ref: "tag:stsci.edu:asdf/transform/transform-1.1.0" + - object: + properties: + wrap_lon_at: + description: Angle at which to wrap the longitude angle. + type: integer + enum: [180, 360] + default: 360 + transform_type: + description: The type of transform/class to initialize. + type: string + enum: [spherical_to_cartesian, cartesian_to_spherical] + + required: [transform_type, wrap_lon_at] diff --git a/gwcs/tags/geometry_models.py b/gwcs/tags/geometry_models.py index db3fa49d..7227826a 100644 --- a/gwcs/tags/geometry_models.py +++ b/gwcs/tags/geometry_models.py @@ -3,10 +3,11 @@ """ from asdf import yamlutil from ..gwcs_types import GWCSTransformType -from .. geometry import ToDirectionCosines, FromDirectionCosines +from .. geometry import (ToDirectionCosines, FromDirectionCosines, + SphericalToCartesian, CartesianToSpherical) -__all__ = ['DirectionCosinesType'] +__all__ = ['DirectionCosinesType', 'SphericalCartesianType'] class DirectionCosinesType(GWCSTransformType): @@ -34,3 +35,35 @@ def to_tree_transform(cls, model, ctx): raise TypeError(f"Model of type {model.__class__} is not supported.") node = {'transform_type': transform_type} return yamlutil.custom_tree_to_tagged_tree(node, ctx) + + +class SphericalCartesianType(GWCSTransformType): + name = "spherical_cartesian" + types = [SphericalToCartesian, CartesianToSpherical] + version = "1.0.0" + + @classmethod + def from_tree_transform(cls, node, ctx): + transform_type = node['transform_type'] + wrap_lon_at = node['wrap_lon_at'] + if transform_type == 'spherical_to_cartesian': + return SphericalToCartesian(wrap_lon_at=wrap_lon_at) + elif transform_type == 'cartesian_to_spherical': + return CartesianToSpherical(wrap_lon_at=wrap_lon_at) + else: + raise TypeError(f"Unknown model_type {transform_type}") + + @classmethod + def to_tree_transform(cls, model, ctx): + if isinstance(model, SphericalToCartesian): + transform_type = 'spherical_to_cartesian' + elif isinstance(model, CartesianToSpherical): + transform_type = 'cartesian_to_spherical' + else: + raise TypeError(f"Model of type {model.__class__} is not supported.") + + node = { + 'transform_type': transform_type, + 'wrap_lon_at': model.wrap_lon_at + } + return yamlutil.custom_tree_to_tagged_tree(node, ctx) diff --git a/gwcs/tags/spectroscopy_models.py b/gwcs/tags/spectroscopy_models.py index 80c54f97..afc3fe40 100644 --- a/gwcs/tags/spectroscopy_models.py +++ b/gwcs/tags/spectroscopy_models.py @@ -2,11 +2,10 @@ ASDF tags for spectroscopy related models. """ -import numpy as np -from numpy.testing import assert_array_equal +from numpy.testing import assert_allclose from astropy import units as u -from astropy.tests.helper import assert_quantity_allclose +from astropy.units.quantity import _unquantify_allclose_arguments from asdf import yamlutil from ..gwcs_types import GWCSTransformType @@ -122,3 +121,14 @@ def assert_equal(cls, a, b): assert isinstance(b, WavelengthFromGratingEquation) # nosec assert_quantity_allclose(a.groove_density, b.groove_density) # nosec assert a.spectral_order.value == b.spectral_order.value # nosec + + +def assert_quantity_allclose(actual, desired, rtol=1.e-7, atol=None, **kwargs): + """ + Raise an assertion if two objects are not equal up to desired tolerance. + + This is a :class:`~astropy.units.Quantity`-aware version of + :func:`numpy.testing.assert_allclose`. + """ + assert_allclose(*_unquantify_allclose_arguments( + actual, desired, rtol, atol), **kwargs) diff --git a/gwcs/tags/tests/test_transforms.py b/gwcs/tags/tests/test_transforms.py index 8489597a..49a0c9b9 100644 --- a/gwcs/tags/tests/test_transforms.py +++ b/gwcs/tags/tests/test_transforms.py @@ -19,9 +19,13 @@ snell = sp.Snell3D() todircos = geometry.ToDirectionCosines() fromdircos = geometry.FromDirectionCosines() +tocart = geometry.SphericalToCartesian() +tospher = geometry.CartesianToSpherical() transforms = [todircos, fromdircos, + tospher, + tocart, snell, sell_glass, sell_zemax, diff --git a/gwcs/tests/conftest.py b/gwcs/tests/conftest.py index fa5b345b..e0fa8085 100644 --- a/gwcs/tests/conftest.py +++ b/gwcs/tests/conftest.py @@ -13,6 +13,7 @@ from .. import coordinate_frames as cf from .. import spectroscopy as sp from .. import wcs +from .. import geometry # frames detector_1d = cf.CoordinateFrame(name='detector', axes_order=(0,), naxes=1, axes_type="detector") @@ -271,3 +272,13 @@ def gwcs_3d_galactic_spectral(): owcs.pixel_shape = (10, 20, 30) return owcs + + +@pytest.fixture +def spher_to_cart(): + return geometry.SphericalToCartesian() + + +@pytest.fixture +def cart_to_spher(): + return geometry.CartesianToSpherical() diff --git a/gwcs/tests/data/miriwcs.asdf b/gwcs/tests/data/miriwcs.asdf new file mode 100644 index 00000000..e70bc620 --- /dev/null +++ b/gwcs/tests/data/miriwcs.asdf @@ -0,0 +1,268 @@ +#ASDF 1.0.0 +#ASDF_STANDARD 1.4.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.5.2} +history: + extensions: + - !core/extension_metadata-1.0.0 + extension_class: astropy.io.misc.asdf.extension.AstropyAsdfExtension + software: {name: astropy, version: '4.0'} + - !core/extension_metadata-1.0.0 + extension_class: gwcs.extension.GWCSExtension + software: {name: gwcs, version: 0.12.0} + - !core/extension_metadata-1.0.0 + extension_class: asdf.extension.BuiltinExtension + software: {name: asdf, version: 2.5.2} +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.1.0 + bounding_box: + - [-0.5, 1023.5] + - [3.5, 1027.5] + forward: + - !transform/concatenate-1.1.0 + forward: + - !transform/shift-1.2.0 {offset: 0.15} + - !transform/shift-1.2.0 {offset: -0.59} + - !transform/compose-1.1.0 + forward: + - !transform/compose-1.1.0 + forward: + - !transform/compose-1.1.0 + forward: + - !transform/compose-1.1.0 + forward: + - !transform/compose-1.1.0 + forward: + - !transform/compose-1.1.0 + forward: + - !transform/compose-1.1.0 + forward: + - !transform/concatenate-1.1.0 + forward: + - !transform/shift-1.2.0 {offset: -4.0} + - !transform/identity-1.1.0 {} + - !transform/concatenate-1.1.0 + forward: + - !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 0 + datatype: float64 + byteorder: little + shape: [2] + inverse: !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 1 + datatype: float64 + byteorder: little + shape: [2] + name: M_column_correction + - !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 2 + datatype: float64 + byteorder: little + shape: [2] + inverse: !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 3 + datatype: float64 + byteorder: little + shape: [2] + name: M_row_correction + - !transform/remap_axes-1.1.0 + inverse: !transform/identity-1.1.0 {n_dims: 2} + mapping: [0, 1, 0, 1] + - !transform/concatenate-1.1.0 + forward: + - !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 4 + datatype: float64 + byteorder: little + shape: [5, 5] + inverse: !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 5 + datatype: float64 + byteorder: little + shape: [5, 5] + name: B_correction + - !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 6 + datatype: float64 + byteorder: little + shape: [5, 5] + inverse: !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 7 + datatype: float64 + byteorder: little + shape: [5, 5] + name: A_correction + - !transform/remap_axes-1.1.0 + inverse: !transform/remap_axes-1.1.0 + mapping: [0, 1, 0, 1] + mapping: [0, 1, 0, 1] + - !transform/concatenate-1.1.0 + forward: + - !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 8 + datatype: float64 + byteorder: little + shape: [2, 2] + name: TI_row_correction + - !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 9 + datatype: float64 + byteorder: little + shape: [2, 2] + name: TI_column_correction + - !transform/identity-1.1.0 + inverse: !transform/remap_axes-1.1.0 + mapping: [0, 1, 0, 1] + n_dims: 2 + - !transform/remap_axes-1.1.0 + mapping: [1, 0] + inverse: !transform/compose-1.1.0 + forward: + - !transform/compose-1.1.0 + forward: + - !transform/remap_axes-1.1.0 + mapping: [1, 0] + - !transform/compose-1.1.0 + forward: + - !transform/remap_axes-1.1.0 + mapping: [0, 1, 0, 1] + - !transform/compose-1.1.0 + forward: + - !transform/concatenate-1.1.0 + forward: + - !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 10 + datatype: float64 + byteorder: little + shape: [2, 2] + name: T_row_correction + - !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 11 + datatype: float64 + byteorder: little + shape: [2, 2] + name: T_column_correction + - !transform/compose-1.1.0 + forward: + - !transform/remap_axes-1.1.0 + mapping: [0, 1, 0, 1] + - !transform/compose-1.1.0 + forward: + - !transform/concatenate-1.1.0 + forward: + - !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 12 + datatype: float64 + byteorder: little + shape: [5, 5] + - !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 13 + datatype: float64 + byteorder: little + shape: [5, 5] + - !transform/compose-1.1.0 + forward: + - !transform/identity-1.1.0 {n_dims: 2} + - !transform/compose-1.1.0 + forward: + - !transform/concatenate-1.1.0 + forward: + - !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 14 + datatype: float64 + byteorder: little + shape: [2] + - !transform/polynomial-1.1.0 + coefficients: !core/ndarray-1.0.0 + source: 15 + datatype: float64 + byteorder: little + shape: [2] + - !transform/concatenate-1.1.0 + forward: + - !transform/shift-1.2.0 {offset: 4.0} + - !transform/identity-1.1.0 {} + - !transform/concatenate-1.1.0 + forward: + - !transform/shift-1.2.0 {offset: -0.15} + - !transform/shift-1.2.0 {offset: 0.59} + - ! + frame: ! + axes_names: [x, y] + axes_order: [0, 1] + axis_physical_types: ['custom:x', 'custom:y'] + name: v2v3 + unit: [!unit/unit-1.0.0 'arcsec', !unit/unit-1.0.0 'arcsec'] + transform: !transform/compose-1.1.0 + forward: + - !transform/concatenate-1.1.0 + forward: + - !transform/scale-1.2.0 {factor: 0.0002777777777777778} + - !transform/scale-1.2.0 {factor: 0.0002777777777777778} + - !transform/rotate_sequence_3d-1.0.0 + angles: [-0.12597594444444443, 0.10374517305555556, 0.0, -72.0545718, -5.630568] + axes_order: zyxyz + name: v23tosky + rotation_type: spherical + - ! + frame: ! + axes_names: [lon, lat] + axes_order: [0, 1] + axis_physical_types: [pos.eq.ra, pos.eq.dec] + name: world + reference_frame: ! + frame_attributes: {} + unit: [!unit/unit-1.0.0 'deg', !unit/unit-1.0.0 'deg'] +... +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Ϝ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@#ASDF BLOCK INDEX +%YAML 1.1 +--- +- 9730 +- 9800 +- 9870 +- 9940 +- 10010 +- 10264 +- 10518 +- 10772 +- 11026 +- 11112 +- 11198 +- 11284 +- 11370 +- 11624 +- 11878 +- 11948 +... diff --git a/gwcs/tests/test_geometry.py b/gwcs/tests/test_geometry.py new file mode 100644 index 00000000..f395f3ef --- /dev/null +++ b/gwcs/tests/test_geometry.py @@ -0,0 +1,233 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +from itertools import product, permutations +import io + +import pytest + +import asdf +from asdf.tests.helpers import assert_roundtrip_tree +import numpy as np +from astropy import units as u + +from .. import geometry + + +_INV_SQRT2 = 1.0 / np.sqrt(2.0) + + +def test_spherical_cartesian_inverse(): + t = geometry.SphericalToCartesian() + assert type(t.inverse) == geometry.CartesianToSpherical + + t = geometry.CartesianToSpherical() + assert type(t.inverse) == geometry.SphericalToCartesian + + +@pytest.mark.parametrize( + 'testval, unit, wrap_at', + product( + [ + (45.0, -90.0, (0.0, 0.0, -1.0)), + (45.0, -45.0, (0.5, 0.5, -_INV_SQRT2)), + (45, 0.0, (_INV_SQRT2, _INV_SQRT2, 0.0)), + (45.0, 45, (0.5, 0.5, _INV_SQRT2)), + (45.0, 90.0, (0.0, 0.0, 1.0)), + (135.0, -90.0, (0.0, 0.0, -1.0)), + (135.0, -45.0, (-0.5, 0.5, -_INV_SQRT2)), + (135.0, 0.0, (-_INV_SQRT2, _INV_SQRT2, 0.0)), + (135.0, 45.0, (-0.5, 0.5, _INV_SQRT2)), + (135.0, 90.0, (0.0, 0.0, 1.0)), + (225.0, -90.0, (0.0, 0.0, -1.0)), + (225.0, -45.0, (-0.5, -0.5, -_INV_SQRT2)), + (225.0, 0.0, (-_INV_SQRT2, -_INV_SQRT2, 0.0)), + (225.0, 45.0, (-0.5, -0.5, _INV_SQRT2)), + (225.0, 90.0, (0.0, 0.0, 1.0)), + (315.0, -90.0, (0.0, 0.0, -1.0)), + (315.0, -45.0, (0.5, -0.5, -_INV_SQRT2)), + (315.0, 0.0, (_INV_SQRT2, -_INV_SQRT2, 0.0)), + (315.0, 45.0, (0.5, -0.5, _INV_SQRT2)), + (315.0, 90.0, (0.0, 0.0, 1.0)), + ], + [1, 1 * u.deg, 3600.0 * u.arcsec, np.pi / 180.0 * u.rad], + [180, 360], + ) +) +def test_spherical_to_cartesian(testval, unit, wrap_at): + s2c = geometry.SphericalToCartesian(wrap_lon_at=wrap_at) + ounit = 1 if unit == 1 else u.dimensionless_unscaled + lon, lat, expected = testval + + if wrap_at == 180: + lon = np.mod(lon - 180.0, 360.0) - 180.0 + + xyz = s2c(lon * unit, lat * unit) + if unit != 1: + assert xyz[0].unit == u.dimensionless_unscaled + assert u.allclose(xyz, u.Quantity(expected, ounit), atol=1e-15 * ounit) + + +@pytest.mark.parametrize( + 'lon, lat, unit, wrap_at', + list(product( + [0, 45, 90, 135, 180, 225, 270, 315, 360], + [-90, -89, -55, 0, 25, 89, 90], + [1, 1 * u.deg, 3600.0 * u.arcsec, np.pi / 180.0 * u.rad], + [180, 360], + )) +) +def test_spher2cart_roundrip(lon, lat, unit, wrap_at): + s2c = geometry.SphericalToCartesian(wrap_lon_at=wrap_at) + c2s = geometry.CartesianToSpherical(wrap_lon_at=wrap_at) + ounit = 1 if unit == 1 else u.deg + + if wrap_at == 180: + lon = np.mod(lon - 180.0, 360.0) - 180.0 + + assert s2c.wrap_lon_at == wrap_at + assert c2s.wrap_lon_at == wrap_at + + assert u.allclose( + c2s(*s2c(lon * unit, lat * unit)), + (lon * ounit, lat * ounit), + atol=1e-15 * ounit + ) + + +def test_cart2spher_at_pole(cart_to_spher): + assert np.allclose(cart_to_spher(0, 0, 1), (0, 90), rtol=0, atol=1e-15) + + +@pytest.mark.parametrize( + 'lonlat, unit, wrap_at', + list(product( + [ + [[1], [-80]], + [[325], [-89]], + [[0, 1, 120, 180, 225, 325, 359], [-89, 0, 89, 10, -15, 45, -30]], + [np.array([0.0, 1, 120, 180, 225, 325, 359]), np.array([-89, 0.0, 89, 10, -15, 45, -30])] + ], + [None, 1 * u.deg], + [180, 360], + )) +) +def test_spher2cart_roundrip_arr(lonlat, unit, wrap_at): + lon, lat = lonlat + s2c = geometry.SphericalToCartesian(wrap_lon_at=wrap_at) + c2s = geometry.CartesianToSpherical(wrap_lon_at=wrap_at) + + if wrap_at == 180: + if isinstance(lon, np.ndarray): + lon = np.mod(lon - 180.0, 360.0) - 180.0 + else: + lon = [((l - 180.0) % 360.0) - 180.0 for l in lon] + + atol = 1e-15 + if unit is None: + olon = lon + olat = lat + else: + olon = lon * u.deg + olat = lat * u.deg + lon = lon * unit + lat = lat * unit + atol = atol * u.deg + + assert u.allclose( + c2s(*s2c(lon, lat)), + (olon, olat), + atol=atol + ) + + +@pytest.mark.parametrize('unit1, unit2', [(u.deg, 1), (1, u.deg)]) +def test_spherical_to_cartesian_mixed_Q(spher_to_cart, unit1, unit2): + with pytest.raises(TypeError) as arg_err: + spher_to_cart(135.0 * unit1, 45.0 * unit2) + assert (arg_err.value.args[0] == "All arguments must be of the same type " + "(i.e., quantity or not).") + + +@pytest.mark.parametrize( + 'x, y, z', + list(set( + tuple(permutations([1 * u.m, 1, 1])) + tuple(permutations([1 * u.m, 1 * u.m, 1])) + )) +) +def test_cartesian_to_spherical_mixed_Q(cart_to_spher, x, y, z): + with pytest.raises(TypeError) as arg_err: + cart_to_spher(x, y, z) + assert (arg_err.value.args[0] == "All arguments must be of the same type " + "(i.e., quantity or not).") + + +@pytest.mark.parametrize('wrap_at', ['1', 180., True, 180j, [180], -180, 0]) +def test_c2s2c_wrong_wrap_type(spher_to_cart, cart_to_spher, wrap_at): + err_msg = "'wrap_lon_at' must be an integer number: 180 or 360" + with pytest.raises(ValueError) as arg_err: + geometry.SphericalToCartesian(wrap_lon_at=wrap_at) + assert arg_err.value.args[0] == err_msg + + with pytest.raises(ValueError) as arg_err: + spher_to_cart.wrap_lon_at = wrap_at + assert arg_err.value.args[0] == err_msg + + with pytest.raises(ValueError) as arg_err: + geometry.CartesianToSpherical(wrap_lon_at=wrap_at) + assert arg_err.value.args[0] == err_msg + + with pytest.raises(ValueError) as arg_err: + cart_to_spher.wrap_lon_at = wrap_at + assert arg_err.value.args[0] == err_msg + + +def test_cartesian_spherical_asdf(tmpdir): + s2c0 = geometry.SphericalToCartesian(wrap_lon_at=360) + c2s0 = geometry.CartesianToSpherical(wrap_lon_at=180) + + # asdf round-trip test: + assert_roundtrip_tree({'c2s': c2s0, 's2c': s2c0}, tmpdir) + + # create file object + f = asdf.AsdfFile({'c2s': c2s0, 's2c': s2c0}) + + # write to... + buf = io.BytesIO() + f.write_to(buf) + + # read back: + buf.seek(0) + f = asdf.open(buf) + + # retrieve transformations: + c2s = f['c2s'] + s2c = f['s2c'] + + pcoords = [ + (45.0, -90.0), (45.0, -45.0), (45, 0.0), + (45.0, 45), (45.0, 90.0), (135.0, -90.0), + (135.0, -45.0), (135.0, 0.0), (135.0, 45.0), + (135.0, 90.0) + ] + + ncoords = [ + (225.0, -90.0), (225.0, -45.0), + (225.0, 0.0), (225.0, 45.0), (225.0, 90.0), + (315.0, -90.0), (315.0, -45.0), (315.0, 0.0), + (315.0, 45.0), (315.0, 90.0) + ] + + for lon, lat in pcoords: + xyz = s2c(lon, lat) + assert xyz == s2c0(lon, lat) + lon2, lat2 = c2s(*xyz) + assert lon2, lat2 == c2s0(*xyz) + assert np.allclose((lon, lat), (lon2, lat2)) + + for lon, lat in ncoords: + xyz = s2c(lon, lat) + assert xyz == s2c0(lon, lat) + lon2, lat2 = c2s(*xyz) + lon3, lat3 = s2c.inverse(*xyz) + assert lon2, lat2 == c2s0(*xyz) + assert np.allclose((lon, lat), (lon2 + 360, lat2)) + assert np.allclose((lon, lat), (lon3, lat2)) diff --git a/gwcs/tests/test_utils.py b/gwcs/tests/test_utils.py index 75b58344..729bded0 100644 --- a/gwcs/tests/test_utils.py +++ b/gwcs/tests/test_utils.py @@ -25,11 +25,12 @@ def test_wrong_projcode2(): gwutils.create_projection_transform("TAM") +@pytest.mark.remote_data def test_fits_transform(): hdr = fits.Header.fromfile(get_pkg_data_filename('data/simple_wcs2.hdr')) gw1 = gwutils.make_fitswcs_transform(hdr) w1 = fitswcs.WCS(hdr) - assert_allclose(gw1(1, 2), w1.wcs_pix2world(1, 2, 1), atol=10 ** -8) + assert_allclose(gw1(1, 2), w1.wcs_pix2world(1, 2, 0), atol=10 ** -8) def test_lon_pole(): @@ -60,12 +61,10 @@ def test_unknown_ctype(): transform = gwutils.make_fitswcs_transform(wcsinfo) x = np.linspace(-5, 7, 10) y = np.linspace(-5, 7, 10) - expected = (np.array([-0.00079444, -0.0007463, -0.00069815, -0.00065, -0.00060185, - -0.0005537, -0.00050556, -0.00045741, -0.00040926, -0.00036111] - ), - np.array([-0.00075833, -0.00071019, -0.00066204, -0.00061389, -0.00056574, - -0.00051759, -0.00046944, -0.0004213, -0.00037315, -0.000325] - ) + expected = (np.array([-0.00075833, -0.00071019, -0.00066204, -0.00061389, -0.00056574, + -0.00051759, -0.00046944, -0.0004213 , -0.00037315, -0.000325]), + np.array([-0.00072222, -0.00067407, -0.00062593, -0.00057778, -0.00052963, + -0.00048148, -0.00043333, -0.00038519, -0.00033704, -0.00028889]) ) a, b = transform(x, y) assert_allclose(a, expected[0], atol=10**-8) diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 73d8ad61..a8702ab8 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -15,7 +15,7 @@ from .. import coordinate_frames as cf from .. import utils from ..utils import CoordinateFrameError - +import asdf m1 = models.Shift(12.4) & models.Shift(-2) m2 = models.Scale(2) & models.Scale(-2) @@ -254,6 +254,7 @@ def test_grid_from_bounding_box_step(): grid_from_bounding_box(bb, step=(1, 2, 1)) +@pytest.mark.remote_data def test_wcs_from_points(): np.random.seed(0) hdr = fits.Header.fromtextfile(get_pkg_data_filename("data/acs.hdr"), endcard=False) @@ -386,6 +387,7 @@ def test_high_level_api(): assert_allclose(z1, xv) +@pytest.mark.remote_data class TestImaging(object): def setup_class(self): hdr = fits.Header.fromtextfile(get_pkg_data_filename("data/acs.hdr"), endcard=False) @@ -484,3 +486,31 @@ def test_pixel_to_world(self): assert isinstance(sky_coord, coord.SkyCoord) assert_allclose(sky_coord.data.lon.value, ra) assert_allclose(sky_coord.data.lat.value, dec) + + +def test_to_fits_sip(): + y, x = np.mgrid[:1024:10, :1024:10] + xflat = np.ravel(x[1:-1, 1:-1]) + yflat = np.ravel(y[1:-1, 1:-1]) + af = asdf.open(get_pkg_data_filename('data/miriwcs.asdf')) + miriwcs = af.tree['wcs'] + bounding_box = ((0, 1024), (0, 1024)) + mirisip = miriwcs.to_fits_sip(bounding_box, max_inv_pix_error=0.1) + fitssip = astwcs.WCS(mirisip) + fitsvalx, fitsvaly = fitssip.all_pix2world(xflat+1, yflat+1, 1) + gwcsvalx, gwcsvaly = miriwcs(xflat, yflat) + assert_allclose(gwcsvalx, fitsvalx, atol=1e-10, rtol=0) + assert_allclose(gwcsvaly, fitsvaly, atol=1e-10, rtol=0) + fits_inverse_valx, fits_inverse_valy = fitssip.all_world2pix(fitsvalx, fitsvaly, 1) + assert_allclose(xflat, fits_inverse_valx - 1, atol=0.1, rtol=0) + assert_allclose(yflat, fits_inverse_valy - 1, atol=0.1, rtol=0) + + mirisip = miriwcs.to_fits_sip(bounding_box=None, max_inv_pix_error=0.1) + fitssip = astwcs.WCS(mirisip) + fitsvalx, fitsvaly = fitssip.all_pix2world(xflat+1, yflat+1, 1) + assert_allclose(gwcsvalx, fitsvalx, atol=1e-10, rtol=0) + assert_allclose(gwcsvaly, fitsvaly, atol=1e-10, rtol=0) + + with pytest.raises(ValueError): + miriwcs.bounding_box = None + mirisip = miriwcs.to_fits_sip(bounding_box=None, max_inv_pix_error=0.1) diff --git a/gwcs/utils.py b/gwcs/utils.py index 849acff0..0d80f06a 100644 --- a/gwcs/utils.py +++ b/gwcs/utils.py @@ -375,7 +375,7 @@ def fitswcs_linear(header): # raise DimensionsError("WCSLinearTransform supports only 2 or 3 dimensions, " # "{0} given".format(wcsaxes)) - translation_models = [astmodels.Shift(-shift, name='crpix' + str(i + 1)) + translation_models = [astmodels.Shift(-(shift - 1), name='crpix' + str(i + 1)) for i, shift in enumerate(crpix)] translation = functools.reduce(lambda x, y: x & y, translation_models) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 0119e9ad..5bd66b88 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -2,14 +2,20 @@ import functools import itertools import numpy as np +import numpy.linalg as npla 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) +from astropy.modeling.fitting import LinearLSQFitter +import astropy.io.fits as fits from .api import GWCSAPIMixin from . import coordinate_frames from .utils import CoordinateFrameError from .utils import _toindex from . import utils +from gwcs import coordinate_frames as cf HAS_FIX_INPUTS = True @@ -634,3 +640,274 @@ def fix_inputs(self, fixed): new_pipeline.append((step0[0], new_transform)) new_pipeline.extend(self.pipeline[1:]) return self.__class__(new_pipeline) + + def to_fits_sip(self, bounding_box=None, max_pix_error=0.25, degree=None, + max_inv_pix_error=0.25, inv_degree=None, + npoints=32, verbose=False): + """ + Construct a SIP-based approximation to the WCS in the form of a FITS header + + This assumes a tangent projection. + + The default mode in using this attempts to achieve roughly 0.25 pixel + accuracy over the whole image. + + Parameters + ---------- + bounding_box : tuple, optional + A pair of tuples, each consisting of two numbers + Represents the range of pixel values in both dimensions + ((xmin, xmax), (ymin, ymax)) + max_pix_error : float, optional + Maximum allowed error over the domain of the pixel array. This + error is the equivalent pixel error that corresponds to the maximum + error in the output coordinate resulting from the fit based on + a nominal plate scale. + degree : int, optional + Degree of the SIP polynomial. If supplied, max_pixel_error is ignored. + max_inv_error : float, optional + Maximum allowed inverse error over the domain of the pixel array + in pixel units. If None, no inverse is generated. + inv_degree : int, optional + Degree of the inverse SIP polynomial. If supplied max_inv_pixel_error + is ignored. + npoints : int, optional + The number of points in each dimension to sample the bounding box + for use in the SIP fit. + verbose : bool, optional + print progress of fits + + Returns + ------- + FITS header with all SIP WCS keywords + + Raises + ------ + ValueError + If the WCS is not 2D, an exception will be raised. If the specified accuracy + (both forward and inverse, both rms and maximum) is not achieved an exception + will be raised. + + Notes + ----- + + Use of this requires a judicious choice of required accuracies. Attempts to use + higher degrees (~7 or higher) will typically fail due floating point problems + that arise with high powers. + + + """ + if not isinstance(self.output_frame, cf.CelestialFrame): + raise ValueError( + "The to_fits_sip method only works with celestial frame transforms") + + transform = self.forward_transform + # Determine reference points. + if bounding_box is None and self.bounding_box is None: + raise ValueError("A bounding_box is needed to proceed.") + if bounding_box is None: + bounding_box = self.bounding_box + + (xmin, xmax), (ymin, ymax) = bounding_box + crpix1 = (xmax - xmin) // 2 + crpix2 = (ymax - ymin) // 2 + crval1, crval2 = transform(crpix1, crpix2) + hdr = fits.Header() + hdr['naxis'] = 2 + hdr['naxis1'] = xmax + hdr['naxis2'] = ymax + hdr['ctype1'] = 'RA---TAN-SIP' + hdr['ctype2'] = 'DEC--TAN-SIP' + hdr['CRPIX1'] = crpix1 + 1 + hdr['CRPIX2'] = crpix2 + 1 + hdr['CRVAL1'] = crval1 + hdr['CRVAL2'] = crval2 + hdr['cd1_1'] = 1 # Placeholders for FITS card order, all will change. + hdr['cd1_2'] = 0 + hdr['cd2_1'] = 0 + hdr['cd2_2'] = 1 + # Now rotate to native system and deproject. Recall that transform + # expects pixels in the original coordinate system, but the SIP + # transform is relative to crpix coordinates, thus the initial shift. + ntransform = ((Shift(crpix1) & Shift(crpix2)) | transform + | RotateCelestial2Native(crval1, crval2, 180) + | Sky2Pix_TAN()) + u, v = _make_sampling_grid(npoints, bounding_box) + undist_x, undist_y = ntransform(u, v) + # Double sampling to check if sampling is sufficient. + ud, vd = _make_sampling_grid(2 * npoints, bounding_box) + undist_xd, undist_yd = ntransform(ud, vd) + # Determine approximate pixel scale in order to compute error threshold + # from the specified pixel error. Computed at the center of the array. + x0, y0 = ntransform(0, 0) + xx, xy = ntransform(1, 0) + yx, yy = ntransform(0, 1) + pixarea = np.abs((xx - x0) * (yy - y0) - (xy - y0) * (yx - x0)) + plate_scale = np.sqrt(pixarea) + max_error = max_pix_error * plate_scale + # The fitting section. + fit_poly_x, fit_poly_y, max_resid = _fit_2D_poly(ntransform, npoints, + degree, max_error, + u, v, undist_x, undist_y, + ud, vd, undist_xd, undist_yd, + verbose=verbose) + # The following is necessary to put the fit into the SIP formalism. + cdmat, sip_poly_x, sip_poly_y = _reform_poly_coefficients(fit_poly_x, fit_poly_y) + # cdmat = np.array([[fit_poly_x.c1_0.value, fit_poly_x.c0_1.value], + # [fit_poly_y.c1_0.value, fit_poly_y.c0_1.value]]) + det = cdmat[0][0] * cdmat[1][1] - cdmat[0][1] * cdmat[1][0] + U = ( cdmat[1][1] * undist_x - cdmat[0][1] * undist_y) / det + V = (-cdmat[1][0] * undist_x + cdmat[0][0] * undist_y) / det + detd = cdmat[0][0] * cdmat[1][1] - cdmat[0][1] * cdmat[1][0] + Ud = ( cdmat[1][1] * undist_xd - cdmat[0][1] * undist_yd) / detd + Vd = (-cdmat[1][0] * undist_xd + cdmat[0][0] * undist_yd) / detd + + if max_inv_pix_error: + fit_inv_poly_u, fit_inv_poly_v, max_inv_resid = _fit_2D_poly(ntransform, + npoints, None, + max_inv_pix_error, + U, V, u-U, v-V, + Ud, Vd, ud-Ud, vd-Vd, + verbose=verbose) + pdegree = fit_poly_x.degree + if pdegree > 1: + hdr['a_order'] = pdegree + hdr['b_order'] = pdegree + _store_2D_coefficients(hdr, sip_poly_x, 'A') + _store_2D_coefficients(hdr, sip_poly_y, 'B') + hdr['sipmxerr'] = (max_resid * plate_scale, 'Max diff from GWCS (equiv pix).') + if max_inv_pix_error: + hdr['sipiverr'] = (max_inv_resid, 'Max diff for inverse (pixels)') + _store_2D_coefficients(hdr, fit_inv_poly_u, 'AP', keeplinear=True) + _store_2D_coefficients(hdr, fit_inv_poly_v, 'BP', keeplinear=True) + if max_inv_pix_error: + ipdegree = fit_inv_poly_u.degree + hdr['ap_order'] = ipdegree + hdr['bp_order'] = ipdegree + else: + hdr['ctype1'] = 'RA---TAN' + hdr['ctype2'] = 'DEC--TAN' + + hdr['cd1_1'] = cdmat[0][0] + hdr['cd1_2'] = cdmat[0][1] + hdr['cd2_1'] = cdmat[1][0] + hdr['cd2_2'] = cdmat[1][1] + return hdr + + +def _fit_2D_poly(ntransform, npoints, degree, max_error, + xin, yin, xout, yout, + xind, yind, xoutd, youtd, + verbose=False): + """ + Fit a pair of ordinary 2D polynomials to the supplied transform. + + """ + llsqfitter = LinearLSQFitter() + + # The case of one pass with the specified polynomial degree + if degree: + deglist = [degree] + else: + deglist = range(10) + prev_max_error = float(np.inf) + if verbose: + print(f'maximum_specified_error: {max_error}') + for deg in deglist: + poly_x = Polynomial2D(degree=deg) + poly_y = Polynomial2D(degree=deg) + fit_poly_x = llsqfitter(poly_x, xin, yin, xout) + fit_poly_y = llsqfitter(poly_y, xin, yin, yout) + max_resid = _compute_distance_residual(xout, yout, + fit_poly_x(xin, yin), + fit_poly_y(xin, yin)) + if max_resid > prev_max_error: + raise RuntimeError('Failed to achieve required error tolerance') + if verbose: + print(f'Degree = {deg}, max_resid = {max_resid}') + if max_resid < max_error: + # Check to see if double sampling meets error requirement. + max_resid = _compute_distance_residual(xoutd, youtd, + fit_poly_x(xind, yind), + fit_poly_y(xind, yind)) + if verbose: + print(f'Double sampling check: maximum residual={max_resid}') + if max_resid < max_error: + if verbose: + print('terminating condition met') + break + return fit_poly_x, fit_poly_y, max_resid + + +def _make_sampling_grid(npoints, bounding_box): + (xmin, xmax), (ymin, ymax) = bounding_box + xsize = xmax - xmin + ysize = ymax - ymin + crpix1 = int(xsize / 2) + crpix2 = int(ysize / 2) + stepsize_x = int(xsize / npoints) + stepsize_y = int(ysize / npoints) + # Ensure last row and column are part of the evaluation grid. + y, x = np.mgrid[: ymax + 1: stepsize_y, : xmax + 1: stepsize_x] + x[:, -1] = xsize - 1 + y[-1, :] = ysize - 1 + u = x - crpix1 + v = y - crpix2 + return u, v + +def _compute_distance_residual(undist_x, undist_y, fit_poly_x, fit_poly_y): + """ + Compute the distance residuals and return the rms and maximum values. + """ + dist = np.sqrt((undist_x - fit_poly_x)**2 + (undist_y - fit_poly_y)**2) + max_resid = dist.max() + return max_resid + + +def _reform_poly_coefficients(fit_poly_x, fit_poly_y): + """ + The fit polynomials must be recombined to align with the SIP decomposition + + The result is the f(u,v) and g(u,v) polynomials, and the CD matrix. + """ + # Extract values for CD matrix and recombining + c11 = fit_poly_x.c1_0.value + c12 = fit_poly_x.c0_1.value + c21 = fit_poly_y.c1_0.value + c22 = fit_poly_y.c0_1.value + sip_poly_x = fit_poly_x.copy() + sip_poly_y = fit_poly_y.copy() + # Force low order coefficients to be 0 as defined in SIP + sip_poly_x.c0_0 = 0 + sip_poly_y.c0_0 = 0 + sip_poly_x.c1_0 = 0 + sip_poly_x.c0_1 = 0 + sip_poly_y.c1_0 = 0 + sip_poly_y.c0_1 = 0 + + cdmat = ((c11, c12), (c21, c22)) + invcdmat = npla.inv(np.array(cdmat)) + degree = fit_poly_x.degree + # Now loop through all remaining coefficients + for i in range(0, degree + 1): + for j in range(0, degree + 1): + if (i + j > 1) and (i + j < degree + 1): + old_x = getattr(fit_poly_x, f'c{i}_{j}').value + old_y = getattr(fit_poly_y, f'c{i}_{j}').value + newcoeff = np.dot(invcdmat, np.array([[old_x], [old_y]])) + setattr(sip_poly_x, f'c{i}_{j}', newcoeff[0, 0]) + setattr(sip_poly_y, f'c{i}_{j}', newcoeff[1, 0]) + + return cdmat, sip_poly_x, sip_poly_y + + +def _store_2D_coefficients(hdr, poly_model, coeff_prefix, keeplinear=False): + """ + Write the polynomial model coefficients to the header. + """ + mindeg = int(not keeplinear) + degree = poly_model.degree + for i in range(0, degree + 1): + for j in range(0, degree + 1): + if (i + j) > mindeg and (i + j < degree + 1): + hdr[f'{coeff_prefix}_{i}_{j}'] = getattr(poly_model, f'c{i}_{j}').value