From f1366a4d35d66a0db17a09b0d86cb64bbdd77797 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Tue, 24 Dec 2019 10:21:07 -0500 Subject: [PATCH 01/38] prepare ghange log for development --- CHANGES.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index bd3b1707..cd1d57fb 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,6 @@ +0.12.1 (Unreleased) +------------------ + 0.12.0 (2019-12-24) ------------------- New Features From ad951ec4f75572c692a670527a685c281afd69e5 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Thu, 26 Dec 2019 14:15:45 -0500 Subject: [PATCH 02/38] Add Spherical<->Cartesian coordinate transformation models --- CHANGES.rst | 7 +- gwcs/geometry.py | 88 ++++++++++++++++- .../gwcs/spherical_cartesian-1.0.0.yaml | 39 ++++++++ gwcs/tags/geometry_models.py | 32 ++++++- gwcs/tags/tests/test_transforms.py | 4 + gwcs/tests/conftest.py | 11 +++ gwcs/tests/test_geometry.py | 94 +++++++++++++++++++ 7 files changed, 270 insertions(+), 5 deletions(-) create mode 100644 gwcs/schemas/stsci.edu/gwcs/spherical_cartesian-1.0.0.yaml create mode 100644 gwcs/tests/test_geometry.py diff --git a/CHANGES.rst b/CHANGES.rst index cd1d57fb..42715849 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,5 +1,10 @@ 0.12.1 (Unreleased) ------------------- +------------------- +New Features +^^^^^^^^^^^^ + +- Added two new transforms - ``SphericalToCartesian`` and ``CartesianToSpherical``. [#275] + 0.12.0 (2019-12-24) ------------------- diff --git a/gwcs/geometry.py b/gwcs/geometry.py index 59071383..9494ffb3 100644 --- a/gwcs/geometry.py +++ b/gwcs/geometry.py @@ -5,9 +5,10 @@ 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 +56,86 @@ 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 ``phi`` being the *azimuthal* angle ``[0, 360)`` + and ``theta`` being the *elevation* angle ``[-90, 90]``. + + """ + _separable = False + + _input_units_strict = True + _input_units_allow_dimensionless = True + + n_inputs = 2 + n_outputs = 3 + + def __init__(self, **kwargs): + super().__init__(**kwargs) + self.inputs = ('phi', 'theta') + self.outputs = ('x', 'y', 'z') + + @staticmethod + def evaluate(phi, theta): + if isinstance(phi, u.Quantity) != isinstance(theta, u.Quantity): + raise TypeError("All arguments must be of the same type " + "(i.e., quantity or not).") + + phi = np.deg2rad(phi) + theta = np.deg2rad(theta) + + cs = np.cos(theta) + x = cs * np.cos(phi) + y = cs * np.sin(phi) + z = np.sin(theta) + + return x, y, z + + def inverse(self): + return CartesianToSpherical() + + @property + def input_units(self): + return {'phi': u.deg, 'theta': u.deg} + + +class CartesianToSpherical(Model): + """ + Convert cartesian coordinates to spherical coordinates on a unit sphere. + Spherical coordinates are assumed to be in degrees with ``phi`` being + the *azimuthal* angle ``[0, 360)`` and ``theta`` being the *elevation* + angle ``[-90, 90]``. + + """ + _separable = False + + n_inputs = 3 + n_outputs = 2 + + def __init__(self, **kwargs): + super().__init__(**kwargs) + self.inputs = ('x', 'y', 'z') + self.outputs = ('phi', 'theta') + + @staticmethod + def evaluate(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) + phi = np.mod( + np.rad2deg(np.arctan2(y, x)), + 360.0 * u.deg if nquant else 360.0 + ) + theta = np.rad2deg(np.arctan2(z, h)) + + return phi, theta + + def inverse(self): + return SphericalToCartesian() 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..5b7e72ca --- /dev/null +++ b/gwcs/schemas/stsci.edu/gwcs/spherical_cartesian-1.0.0.yaml @@ -0,0 +1,39 @@ +%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: + transform_type: + description: | + The type of transform/class to initialize. + type: string + enum: [spherical_to_cartesian, cartesian_to_spherical] + required: [transform_type] diff --git a/gwcs/tags/geometry_models.py b/gwcs/tags/geometry_models.py index db3fa49d..5217bc77 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,30 @@ 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'] + if transform_type == 'spherical_to_cartesian': + return SphericalToCartesian() + elif transform_type == 'cartesian_to_spherical': + return CartesianToSpherical() + 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} + return yamlutil.custom_tree_to_tagged_tree(node, ctx) 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/test_geometry.py b/gwcs/tests/test_geometry.py new file mode 100644 index 00000000..d59a8c3f --- /dev/null +++ b/gwcs/tests/test_geometry.py @@ -0,0 +1,94 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +from itertools import product, permutations + +import numpy as np +import pytest +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', + 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] + ) +) +def test_spherical_to_cartesian(spher_to_cart, testval, unit): + ounit = 1 if unit == 1 else u.dimensionless_unscaled + phi, theta, expected = testval + xyz = spher_to_cart(phi * unit, theta * 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( + 'phi,theta,unit', + list(product( + 45.0 * np.arange(8), + [-90, -55, 0, 25, 90], + [1, 1 * u.deg, 3600.0 * u.arcsec, np.pi / 180.0 * u.rad] + )) +) +def test_spher2cart_roundrip(spher_to_cart, cart_to_spher, phi, theta, unit): + ounit = 1 if unit == 1 else u.deg + assert u.allclose( + cart_to_spher(*spher_to_cart(phi * unit, theta * unit)), + (phi * ounit, theta * ounit), + atol=1e-15 * ounit + ) + + +@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).") From 3241620e9ce5420632e49ff08826d5f85e5a8eed Mon Sep 17 00:00:00 2001 From: Ole Streicher Date: Sun, 5 Jan 2020 21:05:07 +0100 Subject: [PATCH 03/38] Disable tests that require remote data --- gwcs/tests/test_utils.py | 1 + gwcs/tests/test_wcs.py | 2 ++ 2 files changed, 3 insertions(+) diff --git a/gwcs/tests/test_utils.py b/gwcs/tests/test_utils.py index 75b58344..b500b954 100644 --- a/gwcs/tests/test_utils.py +++ b/gwcs/tests/test_utils.py @@ -25,6 +25,7 @@ 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) diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 73d8ad61..bb726f50 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -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) From 0b9e6defc5c493b45760320264869e7c700e4f7a Mon Sep 17 00:00:00 2001 From: Ole Streicher Date: Mon, 13 Jan 2020 11:29:31 +0100 Subject: [PATCH 04/38] Fix NaN assignment in coordinate_frames --- gwcs/coordinate_frames.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gwcs/coordinate_frames.py b/gwcs/coordinate_frames.py index a958abcc..20dc95fb 100644 --- a/gwcs/coordinate_frames.py +++ b/gwcs/coordinate_frames.py @@ -629,11 +629,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(): From 130fe2a2380cbc128732f614babfadb6ff7d0859 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Thu, 30 Jan 2020 18:26:03 -0500 Subject: [PATCH 05/38] Add wrap_phi_at and theta_def parameters to Spherical <-> Cartesian models --- CHANGES.rst | 2 +- gwcs/geometry.py | 143 +++++++++++++-- .../gwcs/spherical_cartesian-1.0.0.yaml | 17 +- gwcs/tags/geometry_models.py | 13 +- gwcs/tests/test_geometry.py | 168 ++++++++++++++++-- 5 files changed, 308 insertions(+), 35 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 42715849..5c395697 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -3,7 +3,7 @@ New Features ^^^^^^^^^^^^ -- Added two new transforms - ``SphericalToCartesian`` and ``CartesianToSpherical``. [#275] +- Added two new transforms - ``SphericalToCartesian`` and ``CartesianToSpherical``. [#275, #284] 0.12.0 (2019-12-24) diff --git a/gwcs/geometry.py b/gwcs/geometry.py index 9494ffb3..54a87895 100644 --- a/gwcs/geometry.py +++ b/gwcs/geometry.py @@ -3,6 +3,7 @@ Spectroscopy related models. """ +import numbers import numpy as np from astropy.modeling.core import Model from astropy import units as u @@ -58,12 +59,18 @@ def inverse(self): return ToDirectionCosines() +def _is_int(n): + return isinstance(n, numbers.Integral) and not isinstance(n, bool) + + 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 ``phi`` being the *azimuthal* angle ``[0, 360)`` - and ``theta`` being the *elevation* angle ``[-90, 90]``. + to be in degrees with ``phi`` being the *azimuthal angle* ``[0, 360)`` + (or ``[-180, 180)``) and ``theta`` being the *elevation angle* + ``[-90, 90]`` or the *inclination (polar) angle* in range + ``[0, 180]`` depending on the used definition. """ _separable = False @@ -74,13 +81,54 @@ class SphericalToCartesian(Model): n_inputs = 2 n_outputs = 3 - def __init__(self, **kwargs): + def __init__(self, wrap_phi_at=360, theta_def='elevation', **kwargs): + """ + Parameters + ---------- + wrap_phi_at : {180, 360}, optional + Specifies the range of the azimuthal angle. When ``wrap_phi_at`` is + 180, azimuthal angle will have a range of ``[-180, 180)`` and + when ``wrap_phi_at`` is 360 (default), the azimuthal angle will have + a range of ``[0, 360)`` + + theta_def : {'elevation', 'inclination', 'polar'}, optional + Specifies the definition used for the ``theta``: 'elevation' from + the reference plane or 'inclination' (or 'polar') angle measured + from the pole. + + """ super().__init__(**kwargs) + self.wrap_phi_at = wrap_phi_at + self.theta_def = theta_def self.inputs = ('phi', 'theta') self.outputs = ('x', 'y', 'z') - @staticmethod - def evaluate(phi, theta): + @property + def wrap_phi_at(self): + return self._wrap_phi_at + + @wrap_phi_at.setter + def wrap_phi_at(self, wrap_angle): + if not _is_int(wrap_angle): + raise TypeError("'wrap_phi_at' must be an integer number: 180 or 360") + if wrap_angle not in [180, 360]: + raise ValueError("Allowed 'wrap_phi_at' values are 180 and 360") + self._wrap_phi_at = wrap_angle + + @property + def theta_def(self): + return self._theta_def + + @theta_def.setter + def theta_def(self, theta_def): + if not isinstance(theta_def, str): + raise TypeError("'theta_def' must be a string.") + if theta_def not in ['elevation', 'inclination', 'polar']: + raise ValueError("Allowed 'theta_def' values are: 'elevation', " + "'inclination', or 'polar'.") + self._theta_def = theta_def + + def evaluate(self, phi, theta): if isinstance(phi, u.Quantity) != isinstance(theta, u.Quantity): raise TypeError("All arguments must be of the same type " "(i.e., quantity or not).") @@ -88,15 +136,24 @@ def evaluate(phi, theta): phi = np.deg2rad(phi) theta = np.deg2rad(theta) - cs = np.cos(theta) + if self._theta_def == 'elevation': + cs = np.cos(theta) + si = np.sin(theta) + else: + cs = np.sin(theta) + si = np.cos(theta) + x = cs * np.cos(phi) y = cs * np.sin(phi) - z = np.sin(theta) + z = si return x, y, z def inverse(self): - return CartesianToSpherical() + return CartesianToSpherical( + wrap_phi_at=self._wrap_phi_at, + theta_def=self._theta_def + ) @property def input_units(self): @@ -107,8 +164,9 @@ class CartesianToSpherical(Model): """ Convert cartesian coordinates to spherical coordinates on a unit sphere. Spherical coordinates are assumed to be in degrees with ``phi`` being - the *azimuthal* angle ``[0, 360)`` and ``theta`` being the *elevation* - angle ``[-90, 90]``. + the *azimuthal angle* ``[0, 360)`` (or ``[-180, 180)``) and ``theta`` + being the *elevation angle* ``[-90, 90]`` or the *inclination (polar) + angle* in range ``[0, 180]`` depending on the used definition. """ _separable = False @@ -116,26 +174,75 @@ class CartesianToSpherical(Model): n_inputs = 3 n_outputs = 2 - def __init__(self, **kwargs): + def __init__(self, wrap_phi_at=360, theta_def='elevation', **kwargs): + """ + Parameters + ---------- + wrap_phi_at : {180, 360}, optional + Specifies the range of the azimuthal angle. When ``wrap_phi_at`` is + 180, azimuthal angle will have a range of ``[-180, 180)`` and + when ``wrap_phi_at`` is 360 (default), the azimuthal angle will have + a range of ``[0, 360)`` + + theta_def : {'elevation', 'inclination', 'polar'}, optional + Specifies the definition used for the ``theta``: 'elevation' from + the reference plane or 'inclination' (or 'polar') angle measured + from the pole. + + """ super().__init__(**kwargs) + self.wrap_phi_at = wrap_phi_at + self.theta_def = theta_def self.inputs = ('x', 'y', 'z') self.outputs = ('phi', 'theta') - @staticmethod - def evaluate(x, y, z): + @property + def wrap_phi_at(self): + return self._wrap_phi_at + + @wrap_phi_at.setter + def wrap_phi_at(self, wrap_angle): + if not _is_int(wrap_angle): + raise TypeError("'wrap_phi_at' must be an integer number: 180 or 360") + if wrap_angle not in [180, 360]: + raise ValueError("Allowed 'wrap_phi_at' values are 180 and 360") + self._wrap_phi_at = wrap_angle + + @property + def theta_def(self): + return self._theta_def + + @theta_def.setter + def theta_def(self, theta_def): + if not isinstance(theta_def, str): + raise TypeError("'theta_def' must be a string.") + if theta_def not in ['elevation', 'inclination', 'polar']: + raise ValueError("Allowed 'theta_def' values are: 'elevation', " + "'inclination', or 'polar'.") + self._theta_def = theta_def + + 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) - phi = np.mod( - np.rad2deg(np.arctan2(y, x)), - 360.0 * u.deg if nquant else 360.0 - ) + phi = np.rad2deg(np.arctan2(y, x)) + if h == 0.0: + phi *= 0.0 + + if self._wrap_phi_at != 180: + phi = np.mod(phi, 360.0 * u.deg if nquant else 360.0) + theta = np.rad2deg(np.arctan2(z, h)) + if self._theta_def != 'elevation': + theta = (90.0 * u.deg if nquant else 90.0) - theta return phi, theta def inverse(self): - return SphericalToCartesian() + return SphericalToCartesian( + wrap_phi_at=self._wrap_phi_at, + theta_def=self._theta_def + ) 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 index 5b7e72ca..5d60edc2 100644 --- a/gwcs/schemas/stsci.edu/gwcs/spherical_cartesian-1.0.0.yaml +++ b/gwcs/schemas/stsci.edu/gwcs/spherical_cartesian-1.0.0.yaml @@ -31,9 +31,20 @@ allOf: - $ref: "tag:stsci.edu:asdf/transform/transform-1.1.0" - object: properties: - transform_type: + wrap_phi_at: + description: Angle to which to wrap the azimuthal angle. + type: integer + enum: [180, 360] + default: 180 + theta_def: description: | - The type of transform/class to initialize. + Definition for the theta angle: polar/inclination or elevation. + type: string + enum: [elevation, inclination, polar] + default: elevation + transform_type: + description: The type of transform/class to initialize. type: string enum: [spherical_to_cartesian, cartesian_to_spherical] - required: [transform_type] + + required: [transform_type, wrap_phi_at, theta_def] diff --git a/gwcs/tags/geometry_models.py b/gwcs/tags/geometry_models.py index 5217bc77..eed6c314 100644 --- a/gwcs/tags/geometry_models.py +++ b/gwcs/tags/geometry_models.py @@ -45,10 +45,12 @@ class SphericalCartesianType(GWCSTransformType): @classmethod def from_tree_transform(cls, node, ctx): transform_type = node['transform_type'] + wrap_phi_at = node['wrap_phi_at'] + theta_def = node['theta_def'] if transform_type == 'spherical_to_cartesian': - return SphericalToCartesian() + return SphericalToCartesian(wrap_phi_at=wrap_phi_at, theta_def=theta_def) elif transform_type == 'cartesian_to_spherical': - return CartesianToSpherical() + return CartesianToSpherical(wrap_phi_at=wrap_phi_at, theta_def=theta_def) else: raise TypeError(f"Unknown model_type {transform_type}") @@ -60,5 +62,10 @@ def to_tree_transform(cls, model, ctx): transform_type = 'cartesian_to_spherical' else: raise TypeError(f"Model of type {model.__class__} is not supported.") - node = {'transform_type': transform_type} + + node = { + 'transform_type': transform_type, + 'wrap_phi_at': model.wrap_phi_at, + 'theta_def': model.theta_def + } return yamlutil.custom_tree_to_tagged_tree(node, ctx) diff --git a/gwcs/tests/test_geometry.py b/gwcs/tests/test_geometry.py index d59a8c3f..88398323 100644 --- a/gwcs/tests/test_geometry.py +++ b/gwcs/tests/test_geometry.py @@ -1,5 +1,7 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst from itertools import product, permutations +import io +import asdf import numpy as np import pytest @@ -20,7 +22,7 @@ def test_spherical_cartesian_inverse(): @pytest.mark.parametrize( - 'testval,unit', + 'testval,unit,wrap_at,theta_def', product( [ (45.0, -90.0, (0.0, 0.0, -1.0)), @@ -44,31 +46,59 @@ def test_spherical_cartesian_inverse(): (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] + [1, 1 * u.deg, 3600.0 * u.arcsec, np.pi / 180.0 * u.rad], + [180, 360], + ['elevation', 'polar'], ) ) -def test_spherical_to_cartesian(spher_to_cart, testval, unit): +def test_spherical_to_cartesian(testval, unit, wrap_at, theta_def): + s2c = geometry.SphericalToCartesian(wrap_phi_at=wrap_at, theta_def=theta_def) ounit = 1 if unit == 1 else u.dimensionless_unscaled phi, theta, expected = testval - xyz = spher_to_cart(phi * unit, theta * unit) + + if wrap_at == 180: + phi = np.mod(phi - 180.0, 360.0) - 180.0 + + if theta_def == 'polar': + theta = 90 - theta + + xyz = s2c(phi * unit, theta * 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( - 'phi,theta,unit', + 'phi,theta,unit,wrap_at,theta_def', list(product( 45.0 * np.arange(8), - [-90, -55, 0, 25, 90], - [1, 1 * u.deg, 3600.0 * u.arcsec, np.pi / 180.0 * u.rad] + [-90, -89, -55, 0, 25, 89, 90], + [1, 1 * u.deg, 3600.0 * u.arcsec, np.pi / 180.0 * u.rad], + [180, 360], + ['elevation', 'polar'], )) ) -def test_spher2cart_roundrip(spher_to_cart, cart_to_spher, phi, theta, unit): +def test_spher2cart_roundrip(phi, theta, unit, wrap_at, theta_def): + s2c = geometry.SphericalToCartesian(wrap_phi_at=wrap_at, theta_def=theta_def) + c2s = geometry.CartesianToSpherical(wrap_phi_at=wrap_at, theta_def=theta_def) ounit = 1 if unit == 1 else u.deg + + if wrap_at == 180: + phi = np.mod(phi - 180.0, 360.0) - 180.0 + + sl = slice(int(theta == 90 or theta == -90), 2) + + if theta_def == 'polar': + theta = 90 - theta + + assert s2c.wrap_phi_at == wrap_at + assert s2c.theta_def == theta_def + assert c2s.wrap_phi_at == wrap_at + assert c2s.theta_def == theta_def + assert u.allclose( - cart_to_spher(*spher_to_cart(phi * unit, theta * unit)), - (phi * ounit, theta * ounit), + c2s(*s2c(phi * unit, theta * unit))[sl], + (phi * ounit, theta * ounit)[sl], atol=1e-15 * ounit ) @@ -92,3 +122,121 @@ def test_cartesian_to_spherical_mixed_Q(cart_to_spher, x, y, z): 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]]) +def test_c2s2c_wrong_wrap_type(spher_to_cart, cart_to_spher, wrap_at): + with pytest.raises(TypeError) as arg_err: + geometry.SphericalToCartesian(wrap_phi_at=wrap_at) + assert (arg_err.value.args[0] == "'wrap_phi_at' must be an integer number: " + "180 or 360") + + with pytest.raises(TypeError) as arg_err: + spher_to_cart.wrap_phi_at = wrap_at + assert (arg_err.value.args[0] == "'wrap_phi_at' must be an integer number: " + "180 or 360") + + with pytest.raises(TypeError) as arg_err: + geometry.SphericalToCartesian(wrap_phi_at=wrap_at) + assert (arg_err.value.args[0] == "'wrap_phi_at' must be an integer number: " + "180 or 360") + + with pytest.raises(TypeError) as arg_err: + cart_to_spher.wrap_phi_at = wrap_at + assert (arg_err.value.args[0] == "'wrap_phi_at' must be an integer number: " + "180 or 360") + + +@pytest.mark.parametrize('wrap_at', [-180, 0, 1]) +def test_c2s2c_wrong_wrap_val(spher_to_cart, cart_to_spher, wrap_at): + with pytest.raises(ValueError) as arg_err: + geometry.SphericalToCartesian(wrap_phi_at=wrap_at) + assert arg_err.value.args[0] == "Allowed 'wrap_phi_at' values are 180 and 360" + + with pytest.raises(ValueError) as arg_err: + spher_to_cart.wrap_phi_at = wrap_at + assert arg_err.value.args[0] == "Allowed 'wrap_phi_at' values are 180 and 360" + + with pytest.raises(ValueError) as arg_err: + geometry.SphericalToCartesian(wrap_phi_at=wrap_at) + assert arg_err.value.args[0] == "Allowed 'wrap_phi_at' values are 180 and 360" + + with pytest.raises(ValueError) as arg_err: + cart_to_spher.wrap_phi_at = wrap_at + assert arg_err.value.args[0] == "Allowed 'wrap_phi_at' values are 180 and 360" + + +@pytest.mark.parametrize('theta_def', [180., True, 180j, [180], np.int(1)]) +def test_c2s2c_wrong_theta_type(spher_to_cart, cart_to_spher, theta_def): + with pytest.raises(TypeError) as arg_err: + geometry.SphericalToCartesian(theta_def=theta_def) + assert arg_err.value.args[0] == "'theta_def' must be a string." + + with pytest.raises(TypeError) as arg_err: + spher_to_cart.theta_def = theta_def + assert arg_err.value.args[0] == "'theta_def' must be a string." + + with pytest.raises(TypeError) as arg_err: + geometry.SphericalToCartesian(theta_def=theta_def) + assert arg_err.value.args[0] == "'theta_def' must be a string." + + with pytest.raises(TypeError) as arg_err: + cart_to_spher.theta_def = theta_def + assert arg_err.value.args[0] == "'theta_def' must be a string." + + +@pytest.mark.parametrize('theta_def', ['ele', '']) +def test_c2s2c_wrong_theta_value(spher_to_cart, cart_to_spher, theta_def): + with pytest.raises(ValueError) as arg_err: + geometry.SphericalToCartesian(theta_def=theta_def) + assert (arg_err.value.args[0] == "Allowed 'theta_def' values are: " + "'elevation', 'inclination', or 'polar'.") + + with pytest.raises(ValueError) as arg_err: + spher_to_cart.theta_def = theta_def + assert (arg_err.value.args[0] == "Allowed 'theta_def' values are: " + "'elevation', 'inclination', or 'polar'.") + + with pytest.raises(ValueError) as arg_err: + geometry.SphericalToCartesian(theta_def=theta_def) + assert (arg_err.value.args[0] == "Allowed 'theta_def' values are: " + "'elevation', 'inclination', or 'polar'.") + + with pytest.raises(ValueError) as arg_err: + cart_to_spher.theta_def = theta_def + assert (arg_err.value.args[0] == "Allowed 'theta_def' values are: " + "'elevation', 'inclination', or 'polar'.") + + +def test_cartesian_spherical_asdf(spher_to_cart, cart_to_spher): + # create file object + f = asdf.AsdfFile({'c2s': cart_to_spher, 's2c': spher_to_cart}) + + # write to... + buf = io.BytesIO() + f.write_to(buf) + + # read back: + buf.seek(0) + f = asdf.open(buf) + + # retreave transformations: + c2s = f['c2s'] + s2c = f['s2c'] + + coords = [ + (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), (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 phi, th in coords: + xyz = s2c(phi, th) + assert xyz == spher_to_cart(phi, th) + phi2, th2 = c2s(*xyz) + assert phi2, th2 == cart_to_spher(*xyz) + assert np.allclose((phi, th), (phi2, th2)) From ee4bfa757db6d32d6d26f47dfd57def073ab0c8b Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Sat, 1 Feb 2020 00:57:06 -0500 Subject: [PATCH 06/38] address reviewers comments --- gwcs/geometry.py | 161 +++++++++++------- .../gwcs/spherical_cartesian-1.0.0.yaml | 6 +- gwcs/tests/test_geometry.py | 86 +++------- 3 files changed, 129 insertions(+), 124 deletions(-) diff --git a/gwcs/geometry.py b/gwcs/geometry.py index 54a87895..9d053597 100644 --- a/gwcs/geometry.py +++ b/gwcs/geometry.py @@ -11,6 +11,14 @@ __all__ = ['ToDirectionCosines', 'FromDirectionCosines', 'SphericalToCartesian', 'CartesianToSpherical'] +# allowable values for angle theta definition in spherical<->Cartesian +# transformations: +_THETA_NAMES = ['latitude', 'colatitude', 'altitude', 'elevation', + 'inclination', 'polar'] + +# theta definitions for which angle is measured from the "reference" plane: +_LAT_LIKE = ['latitude', 'altitude', 'elevation'] + class ToDirectionCosines(Model): """ @@ -59,18 +67,14 @@ def inverse(self): return ToDirectionCosines() -def _is_int(n): - return isinstance(n, numbers.Integral) and not isinstance(n, bool) - - 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 ``phi`` being the *azimuthal angle* ``[0, 360)`` - (or ``[-180, 180)``) and ``theta`` being the *elevation angle* - ``[-90, 90]`` or the *inclination (polar) angle* in range - ``[0, 180]`` depending on the used definition. + to be in degrees with ``phi`` being the *azimuthal angle* (or *longitude*) + ``[0, 360)`` (or ``[-180, 180)``) and ``theta`` being the *elevation angle* + (or *latitude*) ``[-90, 90]`` or the *inclination (polar, colatitude) angle* + in range ``[0, 180]`` depending on the used definition. """ _separable = False @@ -81,52 +85,71 @@ class SphericalToCartesian(Model): n_inputs = 2 n_outputs = 3 - def __init__(self, wrap_phi_at=360, theta_def='elevation', **kwargs): + def __init__(self, wrap_phi_at=360, theta_def='latitude', **kwargs): """ Parameters ---------- - wrap_phi_at : {180, 360}, optional - Specifies the range of the azimuthal angle. When ``wrap_phi_at`` is - 180, azimuthal angle will have a range of ``[-180, 180)`` and - when ``wrap_phi_at`` is 360 (default), the azimuthal angle will have - a range of ``[0, 360)`` - - theta_def : {'elevation', 'inclination', 'polar'}, optional - Specifies the definition used for the ``theta``: 'elevation' from - the reference plane or 'inclination' (or 'polar') angle measured - from the pole. + wrap_phi_at : {360, 180}, optional + An **integer number** that specifies the range of the azimuthal + (longitude) angle. When ``wrap_phi_at`` is 180, azimuthal angle + will have a range of ``[-180, 180)`` and when ``wrap_phi_at`` + is 360 (default), the azimuthal angle will have a range of + ``[0, 360)``. + + theta_def : {'latitude', 'colatitude', 'altitude', 'elevation', \ + 'inclination', 'polar'}, optional + Specifies the definition used for the angle ``theta``: + either 'elevation' angle (synonyms: 'latitude', 'altitude') from + the reference plane or 'inclination' angle (synonyms: 'polar', + 'colatitude') measured from the pole (zenith direction). """ super().__init__(**kwargs) - self.wrap_phi_at = wrap_phi_at - self.theta_def = theta_def self.inputs = ('phi', 'theta') self.outputs = ('x', 'y', 'z') + self.wrap_phi_at = wrap_phi_at + self.theta_def = theta_def @property def wrap_phi_at(self): + """ An **integer number** that specifies the range of the azimuthal + (longitude) angle. + + Allowed values are 180 and 360. When ``wrap_phi_at`` + is 180, azimuthal angle will have a range of ``[-180, 180)`` and when + ``wrap_phi_at`` is 360 (default), the azimuthal angle will have a + range of ``[0, 360)``. + + """ return self._wrap_phi_at @wrap_phi_at.setter def wrap_phi_at(self, wrap_angle): - if not _is_int(wrap_angle): - raise TypeError("'wrap_phi_at' must be an integer number: 180 or 360") - if wrap_angle not in [180, 360]: - raise ValueError("Allowed 'wrap_phi_at' values are 180 and 360") + if not (isinstance(wrap_angle, numbers.Integral) and wrap_angle in [180, 360]): + raise ValueError("'wrap_phi_at' must be an integer number: 180 or 360") self._wrap_phi_at = wrap_angle @property def theta_def(self): + """ Definition used for the ``theta`` angle, i.e., latitude or colatitude. + + When ``theta_def`` is either 'elevation', or 'latitude', or 'altitude', + angle ``theta_def`` is measured from the reference plane. + When ``theta_def`` is either 'inclination', or 'polar', or 'colatitude', + angle ``theta_def`` is measured from the pole (zenith direction). + + """ return self._theta_def @theta_def.setter def theta_def(self, theta_def): - if not isinstance(theta_def, str): - raise TypeError("'theta_def' must be a string.") - if theta_def not in ['elevation', 'inclination', 'polar']: - raise ValueError("Allowed 'theta_def' values are: 'elevation', " - "'inclination', or 'polar'.") + if theta_def not in _THETA_NAMES: + raise ValueError( + "'theta_def' must be a string with one of the following " + "values: {:s}".format(','.join(map(repr, _THETA_NAMES))) + ) self._theta_def = theta_def + self._is_theta_latitude = theta_def in _LAT_LIKE def evaluate(self, phi, theta): if isinstance(phi, u.Quantity) != isinstance(theta, u.Quantity): @@ -136,7 +159,7 @@ def evaluate(self, phi, theta): phi = np.deg2rad(phi) theta = np.deg2rad(theta) - if self._theta_def == 'elevation': + if self._is_theta_latitude: cs = np.cos(theta) si = np.sin(theta) else: @@ -163,10 +186,13 @@ def input_units(self): class CartesianToSpherical(Model): """ Convert cartesian coordinates to spherical coordinates on a unit sphere. - Spherical coordinates are assumed to be in degrees with ``phi`` being - the *azimuthal angle* ``[0, 360)`` (or ``[-180, 180)``) and ``theta`` - being the *elevation angle* ``[-90, 90]`` or the *inclination (polar) - angle* in range ``[0, 180]`` depending on the used definition. + Output spherical coordinates are in degrees. When input cartesian + coordinates are quantities (``Quantity`` objects), output angles + will also be quantities in degrees. Angle ``phi`` is the *azimuthal angle* + (or *longitude*) ``[0, 360)`` (or ``[-180, 180)``) and angle ``theta`` + is the *elevation angle* (or *latitude*) ``[-90, 90]`` or the *inclination + (polar, colatitude) angle* in range ``[0, 180]`` depending on the used + definition (see documentation for ``theta_def``). """ _separable = False @@ -174,52 +200,71 @@ class CartesianToSpherical(Model): n_inputs = 3 n_outputs = 2 - def __init__(self, wrap_phi_at=360, theta_def='elevation', **kwargs): + def __init__(self, wrap_phi_at=360, theta_def='latitude', **kwargs): """ Parameters ---------- - wrap_phi_at : {180, 360}, optional - Specifies the range of the azimuthal angle. When ``wrap_phi_at`` is - 180, azimuthal angle will have a range of ``[-180, 180)`` and - when ``wrap_phi_at`` is 360 (default), the azimuthal angle will have - a range of ``[0, 360)`` - - theta_def : {'elevation', 'inclination', 'polar'}, optional - Specifies the definition used for the ``theta``: 'elevation' from - the reference plane or 'inclination' (or 'polar') angle measured - from the pole. + wrap_phi_at : {360, 180}, optional + An **integer number** that specifies the range of the azimuthal + (longitude) angle. When ``wrap_phi_at`` is 180, azimuthal angle + will have a range of ``[-180, 180)`` and when ``wrap_phi_at`` + is 360 (default), the azimuthal angle will have a range of + ``[0, 360)``. + + theta_def : {'latitude', 'colatitude', 'altitude', 'elevation', \ + 'inclination', 'polar'}, optional + Specifies the definition used for the angle ``theta``: + either 'elevation' angle (synonyms: 'latitude', 'altitude') from + the reference plane or 'inclination' angle (synonyms: 'polar', + 'colatitude') measured from the pole (zenith direction). """ super().__init__(**kwargs) - self.wrap_phi_at = wrap_phi_at - self.theta_def = theta_def self.inputs = ('x', 'y', 'z') self.outputs = ('phi', 'theta') + self.wrap_phi_at = wrap_phi_at + self.theta_def = theta_def @property def wrap_phi_at(self): + """ An **integer number** that specifies the range of the azimuthal + (longitude) angle. + + Allowed values are 180 and 360. When ``wrap_phi_at`` + is 180, azimuthal angle will have a range of ``[-180, 180)`` and when + ``wrap_phi_at`` is 360 (default), the azimuthal angle will have a + range of ``[0, 360)``. + + """ return self._wrap_phi_at @wrap_phi_at.setter def wrap_phi_at(self, wrap_angle): - if not _is_int(wrap_angle): - raise TypeError("'wrap_phi_at' must be an integer number: 180 or 360") - if wrap_angle not in [180, 360]: - raise ValueError("Allowed 'wrap_phi_at' values are 180 and 360") + if not (isinstance(wrap_angle, numbers.Integral) and wrap_angle in [180, 360]): + raise ValueError("'wrap_phi_at' must be an integer number: 180 or 360") self._wrap_phi_at = wrap_angle @property def theta_def(self): + """ Definition used for the ``theta`` angle, i.e., latitude or colatitude. + + When ``theta_def`` is either 'elevation', or 'latitude', or 'altitude', + angle ``theta_def`` is measured from the reference plane. + When ``theta_def`` is either 'inclination', or 'polar', or 'colatitude', + angle ``theta_def`` is measured from the pole (zenith direction). + + """ return self._theta_def @theta_def.setter def theta_def(self, theta_def): - if not isinstance(theta_def, str): - raise TypeError("'theta_def' must be a string.") - if theta_def not in ['elevation', 'inclination', 'polar']: - raise ValueError("Allowed 'theta_def' values are: 'elevation', " - "'inclination', or 'polar'.") + if theta_def not in _THETA_NAMES: + raise ValueError( + "'theta_def' must be a string with one of the following " + "values: {:s}".format(','.join(map(repr, _THETA_NAMES))) + ) self._theta_def = theta_def + self._is_theta_latitude = theta_def in _LAT_LIKE def evaluate(self, x, y, z): nquant = [isinstance(i, u.Quantity) for i in (x, y, z)].count(True) @@ -236,7 +281,7 @@ def evaluate(self, x, y, z): phi = np.mod(phi, 360.0 * u.deg if nquant else 360.0) theta = np.rad2deg(np.arctan2(z, h)) - if self._theta_def != 'elevation': + if not self._is_theta_latitude: theta = (90.0 * u.deg if nquant else 90.0) - theta return phi, theta 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 index 5d60edc2..1894d885 100644 --- a/gwcs/schemas/stsci.edu/gwcs/spherical_cartesian-1.0.0.yaml +++ b/gwcs/schemas/stsci.edu/gwcs/spherical_cartesian-1.0.0.yaml @@ -35,13 +35,13 @@ allOf: description: Angle to which to wrap the azimuthal angle. type: integer enum: [180, 360] - default: 180 + default: 360 theta_def: description: | Definition for the theta angle: polar/inclination or elevation. type: string - enum: [elevation, inclination, polar] - default: elevation + enum: [latitude, colatitude, altitude, elevation, inclination, polar] + default: latitude transform_type: description: The type of transform/class to initialize. type: string diff --git a/gwcs/tests/test_geometry.py b/gwcs/tests/test_geometry.py index 88398323..c9e5dba6 100644 --- a/gwcs/tests/test_geometry.py +++ b/gwcs/tests/test_geometry.py @@ -22,7 +22,7 @@ def test_spherical_cartesian_inverse(): @pytest.mark.parametrize( - 'testval,unit,wrap_at,theta_def', + 'testval, unit, wrap_at, theta_def', product( [ (45.0, -90.0, (0.0, 0.0, -1.0)), @@ -48,7 +48,7 @@ def test_spherical_cartesian_inverse(): ], [1, 1 * u.deg, 3600.0 * u.arcsec, np.pi / 180.0 * u.rad], [180, 360], - ['elevation', 'polar'], + geometry._THETA_NAMES, ) ) def test_spherical_to_cartesian(testval, unit, wrap_at, theta_def): @@ -59,7 +59,7 @@ def test_spherical_to_cartesian(testval, unit, wrap_at, theta_def): if wrap_at == 180: phi = np.mod(phi - 180.0, 360.0) - 180.0 - if theta_def == 'polar': + if theta_def not in geometry._LAT_LIKE: theta = 90 - theta xyz = s2c(phi * unit, theta * unit) @@ -69,7 +69,7 @@ def test_spherical_to_cartesian(testval, unit, wrap_at, theta_def): @pytest.mark.parametrize( - 'phi,theta,unit,wrap_at,theta_def', + 'phi, theta, unit, wrap_at, theta_def', list(product( 45.0 * np.arange(8), [-90, -89, -55, 0, 25, 89, 90], @@ -88,7 +88,7 @@ def test_spher2cart_roundrip(phi, theta, unit, wrap_at, theta_def): sl = slice(int(theta == 90 or theta == -90), 2) - if theta_def == 'polar': + if theta_def not in geometry._LAT_LIKE: theta = 90 - theta assert s2c.wrap_phi_at == wrap_at @@ -124,88 +124,48 @@ def test_cartesian_to_spherical_mixed_Q(cart_to_spher, x, y, z): "(i.e., quantity or not).") -@pytest.mark.parametrize('wrap_at', ['1', 180., True, 180j, [180]]) +@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): - with pytest.raises(TypeError) as arg_err: - geometry.SphericalToCartesian(wrap_phi_at=wrap_at) - assert (arg_err.value.args[0] == "'wrap_phi_at' must be an integer number: " - "180 or 360") - - with pytest.raises(TypeError) as arg_err: - spher_to_cart.wrap_phi_at = wrap_at - assert (arg_err.value.args[0] == "'wrap_phi_at' must be an integer number: " - "180 or 360") - - with pytest.raises(TypeError) as arg_err: - geometry.SphericalToCartesian(wrap_phi_at=wrap_at) - assert (arg_err.value.args[0] == "'wrap_phi_at' must be an integer number: " - "180 or 360") - - with pytest.raises(TypeError) as arg_err: - cart_to_spher.wrap_phi_at = wrap_at - assert (arg_err.value.args[0] == "'wrap_phi_at' must be an integer number: " - "180 or 360") - - -@pytest.mark.parametrize('wrap_at', [-180, 0, 1]) -def test_c2s2c_wrong_wrap_val(spher_to_cart, cart_to_spher, wrap_at): + err_msg = "'wrap_phi_at' must be an integer number: 180 or 360" with pytest.raises(ValueError) as arg_err: geometry.SphericalToCartesian(wrap_phi_at=wrap_at) - assert arg_err.value.args[0] == "Allowed 'wrap_phi_at' values are 180 and 360" + assert arg_err.value.args[0] == err_msg with pytest.raises(ValueError) as arg_err: spher_to_cart.wrap_phi_at = wrap_at - assert arg_err.value.args[0] == "Allowed 'wrap_phi_at' values are 180 and 360" + assert arg_err.value.args[0] == err_msg with pytest.raises(ValueError) as arg_err: - geometry.SphericalToCartesian(wrap_phi_at=wrap_at) - assert arg_err.value.args[0] == "Allowed 'wrap_phi_at' values are 180 and 360" + geometry.CartesianToSpherical(wrap_phi_at=wrap_at) + assert arg_err.value.args[0] == err_msg with pytest.raises(ValueError) as arg_err: cart_to_spher.wrap_phi_at = wrap_at - assert arg_err.value.args[0] == "Allowed 'wrap_phi_at' values are 180 and 360" + assert arg_err.value.args[0] == err_msg -@pytest.mark.parametrize('theta_def', [180., True, 180j, [180], np.int(1)]) -def test_c2s2c_wrong_theta_type(spher_to_cart, cart_to_spher, theta_def): - with pytest.raises(TypeError) as arg_err: - geometry.SphericalToCartesian(theta_def=theta_def) - assert arg_err.value.args[0] == "'theta_def' must be a string." - - with pytest.raises(TypeError) as arg_err: - spher_to_cart.theta_def = theta_def - assert arg_err.value.args[0] == "'theta_def' must be a string." - - with pytest.raises(TypeError) as arg_err: - geometry.SphericalToCartesian(theta_def=theta_def) - assert arg_err.value.args[0] == "'theta_def' must be a string." - - with pytest.raises(TypeError) as arg_err: - cart_to_spher.theta_def = theta_def - assert arg_err.value.args[0] == "'theta_def' must be a string." - - -@pytest.mark.parametrize('theta_def', ['ele', '']) +@pytest.mark.parametrize( + 'theta_def', + [180., True, 180j, 'ele', ''] +) def test_c2s2c_wrong_theta_value(spher_to_cart, cart_to_spher, theta_def): + err_msg = ("'theta_def' must be a string with one of the following " + "values: {:s}".format(','.join(map(repr, geometry._THETA_NAMES)))) with pytest.raises(ValueError) as arg_err: geometry.SphericalToCartesian(theta_def=theta_def) - assert (arg_err.value.args[0] == "Allowed 'theta_def' values are: " - "'elevation', 'inclination', or 'polar'.") + assert arg_err.value.args[0] == err_msg with pytest.raises(ValueError) as arg_err: spher_to_cart.theta_def = theta_def - assert (arg_err.value.args[0] == "Allowed 'theta_def' values are: " - "'elevation', 'inclination', or 'polar'.") + assert arg_err.value.args[0] == err_msg with pytest.raises(ValueError) as arg_err: - geometry.SphericalToCartesian(theta_def=theta_def) - assert (arg_err.value.args[0] == "Allowed 'theta_def' values are: " - "'elevation', 'inclination', or 'polar'.") + geometry.CartesianToSpherical(theta_def=theta_def) + assert arg_err.value.args[0] == err_msg with pytest.raises(ValueError) as arg_err: cart_to_spher.theta_def = theta_def - assert (arg_err.value.args[0] == "Allowed 'theta_def' values are: " - "'elevation', 'inclination', or 'polar'.") + assert arg_err.value.args[0] == err_msg def test_cartesian_spherical_asdf(spher_to_cart, cart_to_spher): From 251784e04a419e6b05c1f9cd584e2c582ab6a126 Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Mon, 3 Feb 2020 15:27:48 -0500 Subject: [PATCH 07/38] Remove support for inclination angle --- gwcs/geometry.py | 210 ++++++------------ .../gwcs/spherical_cartesian-1.0.0.yaml | 12 +- gwcs/tags/geometry_models.py | 10 +- gwcs/tests/test_geometry.py | 119 +++++----- 4 files changed, 123 insertions(+), 228 deletions(-) diff --git a/gwcs/geometry.py b/gwcs/geometry.py index 9d053597..f0439116 100644 --- a/gwcs/geometry.py +++ b/gwcs/geometry.py @@ -11,14 +11,6 @@ __all__ = ['ToDirectionCosines', 'FromDirectionCosines', 'SphericalToCartesian', 'CartesianToSpherical'] -# allowable values for angle theta definition in spherical<->Cartesian -# transformations: -_THETA_NAMES = ['latitude', 'colatitude', 'altitude', 'elevation', - 'inclination', 'polar'] - -# theta definitions for which angle is measured from the "reference" plane: -_LAT_LIKE = ['latitude', 'altitude', 'elevation'] - class ToDirectionCosines(Model): """ @@ -70,11 +62,10 @@ def inverse(self): 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 ``phi`` being the *azimuthal angle* (or *longitude*) - ``[0, 360)`` (or ``[-180, 180)``) and ``theta`` being the *elevation angle* - (or *latitude*) ``[-90, 90]`` or the *inclination (polar, colatitude) angle* - in range ``[0, 180]`` depending on the used definition. + 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 @@ -85,102 +76,63 @@ class SphericalToCartesian(Model): n_inputs = 2 n_outputs = 3 - def __init__(self, wrap_phi_at=360, theta_def='latitude', **kwargs): + def __init__(self, wrap_lon_at=360, **kwargs): """ Parameters ---------- - wrap_phi_at : {360, 180}, optional - An **integer number** that specifies the range of the azimuthal - (longitude) angle. When ``wrap_phi_at`` is 180, azimuthal angle - will have a range of ``[-180, 180)`` and when ``wrap_phi_at`` - is 360 (default), the azimuthal angle will have a range of + 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)``. - theta_def : {'latitude', 'colatitude', 'altitude', 'elevation', \ - 'inclination', 'polar'}, optional - Specifies the definition used for the angle ``theta``: - either 'elevation' angle (synonyms: 'latitude', 'altitude') from - the reference plane or 'inclination' angle (synonyms: 'polar', - 'colatitude') measured from the pole (zenith direction). - """ super().__init__(**kwargs) - self.inputs = ('phi', 'theta') + self.inputs = ('lon', 'lat') self.outputs = ('x', 'y', 'z') - self.wrap_phi_at = wrap_phi_at - self.theta_def = theta_def + self.wrap_lon_at = wrap_lon_at @property - def wrap_phi_at(self): - """ An **integer number** that specifies the range of the azimuthal - (longitude) angle. + 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_phi_at`` - is 180, azimuthal angle will have a range of ``[-180, 180)`` and when - ``wrap_phi_at`` is 360 (default), the azimuthal angle will have a + 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_phi_at + return self._wrap_lon_at - @wrap_phi_at.setter - def wrap_phi_at(self, wrap_angle): + @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_phi_at' must be an integer number: 180 or 360") - self._wrap_phi_at = wrap_angle + raise ValueError("'wrap_lon_at' must be an integer number: 180 or 360") + self._wrap_lon_at = wrap_angle - @property - def theta_def(self): - """ Definition used for the ``theta`` angle, i.e., latitude or colatitude. - - When ``theta_def`` is either 'elevation', or 'latitude', or 'altitude', - angle ``theta_def`` is measured from the reference plane. - When ``theta_def`` is either 'inclination', or 'polar', or 'colatitude', - angle ``theta_def`` is measured from the pole (zenith direction). - - """ - return self._theta_def - - @theta_def.setter - def theta_def(self, theta_def): - if theta_def not in _THETA_NAMES: - raise ValueError( - "'theta_def' must be a string with one of the following " - "values: {:s}".format(','.join(map(repr, _THETA_NAMES))) - ) - self._theta_def = theta_def - self._is_theta_latitude = theta_def in _LAT_LIKE - - def evaluate(self, phi, theta): - if isinstance(phi, u.Quantity) != isinstance(theta, u.Quantity): + 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).") - phi = np.deg2rad(phi) - theta = np.deg2rad(theta) + lon = np.deg2rad(lon) + lat = np.deg2rad(lat) - if self._is_theta_latitude: - cs = np.cos(theta) - si = np.sin(theta) - else: - cs = np.sin(theta) - si = np.cos(theta) - - x = cs * np.cos(phi) - y = cs * np.sin(phi) - z = si + 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_phi_at=self._wrap_phi_at, - theta_def=self._theta_def - ) + return CartesianToSpherical(wrap_lon_at=self._wrap_lon_at) @property def input_units(self): - return {'phi': u.deg, 'theta': u.deg} + return {'lon': u.deg, 'lat': u.deg} class CartesianToSpherical(Model): @@ -188,11 +140,10 @@ 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 ``phi`` is the *azimuthal angle* - (or *longitude*) ``[0, 360)`` (or ``[-180, 180)``) and angle ``theta`` - is the *elevation angle* (or *latitude*) ``[-90, 90]`` or the *inclination - (polar, colatitude) angle* in range ``[0, 180]`` depending on the used - definition (see documentation for ``theta_def``). + 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 @@ -200,71 +151,41 @@ class CartesianToSpherical(Model): n_inputs = 3 n_outputs = 2 - def __init__(self, wrap_phi_at=360, theta_def='latitude', **kwargs): + def __init__(self, wrap_lon_at=360, **kwargs): """ Parameters ---------- - wrap_phi_at : {360, 180}, optional - An **integer number** that specifies the range of the azimuthal - (longitude) angle. When ``wrap_phi_at`` is 180, azimuthal angle - will have a range of ``[-180, 180)`` and when ``wrap_phi_at`` - is 360 (default), the azimuthal angle will have a range of + 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)``. - theta_def : {'latitude', 'colatitude', 'altitude', 'elevation', \ - 'inclination', 'polar'}, optional - Specifies the definition used for the angle ``theta``: - either 'elevation' angle (synonyms: 'latitude', 'altitude') from - the reference plane or 'inclination' angle (synonyms: 'polar', - 'colatitude') measured from the pole (zenith direction). - """ super().__init__(**kwargs) self.inputs = ('x', 'y', 'z') - self.outputs = ('phi', 'theta') - self.wrap_phi_at = wrap_phi_at - self.theta_def = theta_def + self.outputs = ('lon', 'lat') + self.wrap_lon_at = wrap_lon_at @property - def wrap_phi_at(self): - """ An **integer number** that specifies the range of the azimuthal - (longitude) angle. + 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_phi_at`` - is 180, azimuthal angle will have a range of ``[-180, 180)`` and when - ``wrap_phi_at`` is 360 (default), the azimuthal angle will have a + 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_phi_at + return self._wrap_lon_at - @wrap_phi_at.setter - def wrap_phi_at(self, wrap_angle): + @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_phi_at' must be an integer number: 180 or 360") - self._wrap_phi_at = wrap_angle - - @property - def theta_def(self): - """ Definition used for the ``theta`` angle, i.e., latitude or colatitude. - - When ``theta_def`` is either 'elevation', or 'latitude', or 'altitude', - angle ``theta_def`` is measured from the reference plane. - When ``theta_def`` is either 'inclination', or 'polar', or 'colatitude', - angle ``theta_def`` is measured from the pole (zenith direction). - - """ - return self._theta_def - - @theta_def.setter - def theta_def(self, theta_def): - if theta_def not in _THETA_NAMES: - raise ValueError( - "'theta_def' must be a string with one of the following " - "values: {:s}".format(','.join(map(repr, _THETA_NAMES))) - ) - self._theta_def = theta_def - self._is_theta_latitude = theta_def in _LAT_LIKE + 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) @@ -273,21 +194,16 @@ def evaluate(self, x, y, z): "(i.e., quantity or not).") h = np.hypot(x, y) - phi = np.rad2deg(np.arctan2(y, x)) + lon = np.rad2deg(np.arctan2(y, x)) if h == 0.0: - phi *= 0.0 + lon *= 0.0 - if self._wrap_phi_at != 180: - phi = np.mod(phi, 360.0 * u.deg if nquant else 360.0) + if self._wrap_lon_at != 180: + lon = np.mod(lon, 360.0 * u.deg if nquant else 360.0) - theta = np.rad2deg(np.arctan2(z, h)) - if not self._is_theta_latitude: - theta = (90.0 * u.deg if nquant else 90.0) - theta + lat = np.rad2deg(np.arctan2(z, h)) - return phi, theta + return lon, lat def inverse(self): - return SphericalToCartesian( - wrap_phi_at=self._wrap_phi_at, - theta_def=self._theta_def - ) + 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 index 1894d885..61f959eb 100644 --- a/gwcs/schemas/stsci.edu/gwcs/spherical_cartesian-1.0.0.yaml +++ b/gwcs/schemas/stsci.edu/gwcs/spherical_cartesian-1.0.0.yaml @@ -31,20 +31,14 @@ allOf: - $ref: "tag:stsci.edu:asdf/transform/transform-1.1.0" - object: properties: - wrap_phi_at: - description: Angle to which to wrap the azimuthal angle. + wrap_lon_at: + description: Angle at which to wrap the longitude angle. type: integer enum: [180, 360] default: 360 - theta_def: - description: | - Definition for the theta angle: polar/inclination or elevation. - type: string - enum: [latitude, colatitude, altitude, elevation, inclination, polar] - default: latitude transform_type: description: The type of transform/class to initialize. type: string enum: [spherical_to_cartesian, cartesian_to_spherical] - required: [transform_type, wrap_phi_at, theta_def] + required: [transform_type, wrap_lon_at] diff --git a/gwcs/tags/geometry_models.py b/gwcs/tags/geometry_models.py index eed6c314..7227826a 100644 --- a/gwcs/tags/geometry_models.py +++ b/gwcs/tags/geometry_models.py @@ -45,12 +45,11 @@ class SphericalCartesianType(GWCSTransformType): @classmethod def from_tree_transform(cls, node, ctx): transform_type = node['transform_type'] - wrap_phi_at = node['wrap_phi_at'] - theta_def = node['theta_def'] + wrap_lon_at = node['wrap_lon_at'] if transform_type == 'spherical_to_cartesian': - return SphericalToCartesian(wrap_phi_at=wrap_phi_at, theta_def=theta_def) + return SphericalToCartesian(wrap_lon_at=wrap_lon_at) elif transform_type == 'cartesian_to_spherical': - return CartesianToSpherical(wrap_phi_at=wrap_phi_at, theta_def=theta_def) + return CartesianToSpherical(wrap_lon_at=wrap_lon_at) else: raise TypeError(f"Unknown model_type {transform_type}") @@ -65,7 +64,6 @@ def to_tree_transform(cls, model, ctx): node = { 'transform_type': transform_type, - 'wrap_phi_at': model.wrap_phi_at, - 'theta_def': model.theta_def + 'wrap_lon_at': model.wrap_lon_at } return yamlutil.custom_tree_to_tagged_tree(node, ctx) diff --git a/gwcs/tests/test_geometry.py b/gwcs/tests/test_geometry.py index c9e5dba6..d4516dc7 100644 --- a/gwcs/tests/test_geometry.py +++ b/gwcs/tests/test_geometry.py @@ -1,10 +1,12 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst from itertools import product, permutations import io -import asdf -import numpy as np 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 @@ -22,7 +24,7 @@ def test_spherical_cartesian_inverse(): @pytest.mark.parametrize( - 'testval, unit, wrap_at, theta_def', + 'testval, unit, wrap_at', product( [ (45.0, -90.0, (0.0, 0.0, -1.0)), @@ -48,57 +50,47 @@ def test_spherical_cartesian_inverse(): ], [1, 1 * u.deg, 3600.0 * u.arcsec, np.pi / 180.0 * u.rad], [180, 360], - geometry._THETA_NAMES, ) ) -def test_spherical_to_cartesian(testval, unit, wrap_at, theta_def): - s2c = geometry.SphericalToCartesian(wrap_phi_at=wrap_at, theta_def=theta_def) +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 - phi, theta, expected = testval + lon, lat, expected = testval if wrap_at == 180: - phi = np.mod(phi - 180.0, 360.0) - 180.0 - - if theta_def not in geometry._LAT_LIKE: - theta = 90 - theta + lon = np.mod(lon - 180.0, 360.0) - 180.0 - xyz = s2c(phi * unit, theta * unit) + 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( - 'phi, theta, unit, wrap_at, theta_def', + 'lon, lat, unit, wrap_at', list(product( 45.0 * np.arange(8), [-90, -89, -55, 0, 25, 89, 90], [1, 1 * u.deg, 3600.0 * u.arcsec, np.pi / 180.0 * u.rad], [180, 360], - ['elevation', 'polar'], )) ) -def test_spher2cart_roundrip(phi, theta, unit, wrap_at, theta_def): - s2c = geometry.SphericalToCartesian(wrap_phi_at=wrap_at, theta_def=theta_def) - c2s = geometry.CartesianToSpherical(wrap_phi_at=wrap_at, theta_def=theta_def) +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: - phi = np.mod(phi - 180.0, 360.0) - 180.0 + lon = np.mod(lon - 180.0, 360.0) - 180.0 - sl = slice(int(theta == 90 or theta == -90), 2) + sl = slice(int(lat == 90 or lat == -90), 2) - if theta_def not in geometry._LAT_LIKE: - theta = 90 - theta - - assert s2c.wrap_phi_at == wrap_at - assert s2c.theta_def == theta_def - assert c2s.wrap_phi_at == wrap_at - assert c2s.theta_def == theta_def + assert s2c.wrap_lon_at == wrap_at + assert c2s.wrap_lon_at == wrap_at assert u.allclose( - c2s(*s2c(phi * unit, theta * unit))[sl], - (phi * ounit, theta * ounit)[sl], + c2s(*s2c(lon * unit, lat * unit))[sl], + (lon * ounit, lat * ounit)[sl], atol=1e-15 * ounit ) @@ -126,51 +118,33 @@ def test_cartesian_to_spherical_mixed_Q(cart_to_spher, x, y, z): @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_phi_at' must be an integer number: 180 or 360" - with pytest.raises(ValueError) as arg_err: - geometry.SphericalToCartesian(wrap_phi_at=wrap_at) - assert arg_err.value.args[0] == err_msg - - with pytest.raises(ValueError) as arg_err: - spher_to_cart.wrap_phi_at = wrap_at - assert arg_err.value.args[0] == err_msg - + err_msg = "'wrap_lon_at' must be an integer number: 180 or 360" with pytest.raises(ValueError) as arg_err: - geometry.CartesianToSpherical(wrap_phi_at=wrap_at) + geometry.SphericalToCartesian(wrap_lon_at=wrap_at) assert arg_err.value.args[0] == err_msg with pytest.raises(ValueError) as arg_err: - cart_to_spher.wrap_phi_at = wrap_at + spher_to_cart.wrap_lon_at = wrap_at assert arg_err.value.args[0] == err_msg - -@pytest.mark.parametrize( - 'theta_def', - [180., True, 180j, 'ele', ''] -) -def test_c2s2c_wrong_theta_value(spher_to_cart, cart_to_spher, theta_def): - err_msg = ("'theta_def' must be a string with one of the following " - "values: {:s}".format(','.join(map(repr, geometry._THETA_NAMES)))) with pytest.raises(ValueError) as arg_err: - geometry.SphericalToCartesian(theta_def=theta_def) + geometry.CartesianToSpherical(wrap_lon_at=wrap_at) assert arg_err.value.args[0] == err_msg with pytest.raises(ValueError) as arg_err: - spher_to_cart.theta_def = theta_def + cart_to_spher.wrap_lon_at = wrap_at assert arg_err.value.args[0] == err_msg - with pytest.raises(ValueError) as arg_err: - geometry.CartesianToSpherical(theta_def=theta_def) - assert arg_err.value.args[0] == err_msg - with pytest.raises(ValueError) as arg_err: - cart_to_spher.theta_def = theta_def - 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) -def test_cartesian_spherical_asdf(spher_to_cart, cart_to_spher): # create file object - f = asdf.AsdfFile({'c2s': cart_to_spher, 's2c': spher_to_cart}) + f = asdf.AsdfFile({'c2s': c2s0, 's2c': s2c0}) # write to... buf = io.BytesIO() @@ -180,23 +154,36 @@ def test_cartesian_spherical_asdf(spher_to_cart, cart_to_spher): buf.seek(0) f = asdf.open(buf) - # retreave transformations: + # retrieve transformations: c2s = f['c2s'] s2c = f['s2c'] - coords = [ + 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), (225.0, -90.0), (225.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 phi, th in coords: - xyz = s2c(phi, th) - assert xyz == spher_to_cart(phi, th) - phi2, th2 = c2s(*xyz) - assert phi2, th2 == cart_to_spher(*xyz) - assert np.allclose((phi, th), (phi2, th2)) + 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)) From 6c8e900fd3bf8be7d034468ce3a3801eaf88251d Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Sun, 9 Feb 2020 00:46:09 -0500 Subject: [PATCH 08/38] Cartesian2Spherical crashes for ndarray inputs --- CHANGES.rst | 3 ++- gwcs/geometry.py | 8 ++++-- gwcs/tests/test_geometry.py | 54 +++++++++++++++++++++++++++++++++---- 3 files changed, 57 insertions(+), 8 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 5c395697..5374f607 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -3,7 +3,8 @@ New Features ^^^^^^^^^^^^ -- Added two new transforms - ``SphericalToCartesian`` and ``CartesianToSpherical``. [#275, #284] +- Added two new transforms - ``SphericalToCartesian`` and + ``CartesianToSpherical``. [#275, #284, #285] 0.12.0 (2019-12-24) diff --git a/gwcs/geometry.py b/gwcs/geometry.py index f0439116..6bf3bc0a 100644 --- a/gwcs/geometry.py +++ b/gwcs/geometry.py @@ -148,6 +148,9 @@ class CartesianToSpherical(Model): """ _separable = False + _input_units_strict = True + _input_units_allow_dimensionless = True + n_inputs = 3 n_outputs = 2 @@ -195,8 +198,9 @@ def evaluate(self, x, y, z): h = np.hypot(x, y) lon = np.rad2deg(np.arctan2(y, x)) - if h == 0.0: - lon *= 0.0 + + if isinstance(h, np.ndarray): + lon[h == 0] *= 0 if self._wrap_lon_at != 180: lon = np.mod(lon, 360.0 * u.deg if nquant else 360.0) diff --git a/gwcs/tests/test_geometry.py b/gwcs/tests/test_geometry.py index d4516dc7..f395f3ef 100644 --- a/gwcs/tests/test_geometry.py +++ b/gwcs/tests/test_geometry.py @@ -69,7 +69,7 @@ def test_spherical_to_cartesian(testval, unit, wrap_at): @pytest.mark.parametrize( 'lon, lat, unit, wrap_at', list(product( - 45.0 * np.arange(8), + [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], @@ -83,18 +83,62 @@ def test_spher2cart_roundrip(lon, lat, unit, wrap_at): if wrap_at == 180: lon = np.mod(lon - 180.0, 360.0) - 180.0 - sl = slice(int(lat == 90 or lat == -90), 2) - assert s2c.wrap_lon_at == wrap_at assert c2s.wrap_lon_at == wrap_at assert u.allclose( - c2s(*s2c(lon * unit, lat * unit))[sl], - (lon * ounit, lat * ounit)[sl], + 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: From 25ba0310837216eaf740a000683ce6e46f7028cb Mon Sep 17 00:00:00 2001 From: Mihai Cara Date: Mon, 10 Feb 2020 12:59:20 -0500 Subject: [PATCH 09/38] address reviewers comments --- gwcs/geometry.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/gwcs/geometry.py b/gwcs/geometry.py index 6bf3bc0a..bcb9f79b 100644 --- a/gwcs/geometry.py +++ b/gwcs/geometry.py @@ -70,7 +70,6 @@ class SphericalToCartesian(Model): """ _separable = False - _input_units_strict = True _input_units_allow_dimensionless = True n_inputs = 2 @@ -148,7 +147,6 @@ class CartesianToSpherical(Model): """ _separable = False - _input_units_strict = True _input_units_allow_dimensionless = True n_inputs = 3 @@ -197,16 +195,13 @@ def evaluate(self, x, y, z): "(i.e., quantity or not).") h = np.hypot(x, y) + lat = np.rad2deg(np.arctan2(z, h)) lon = np.rad2deg(np.arctan2(y, x)) - - if isinstance(h, np.ndarray): - lon[h == 0] *= 0 + lon[h == 0] *= 0 if self._wrap_lon_at != 180: lon = np.mod(lon, 360.0 * u.deg if nquant else 360.0) - lat = np.rad2deg(np.arctan2(z, h)) - return lon, lat def inverse(self): From 33039e60ff2a147b20d50ca1fcd28ba7e5626db2 Mon Sep 17 00:00:00 2001 From: perrygreenfield Date: Tue, 25 Feb 2020 14:12:14 -0500 Subject: [PATCH 10/38] first working version of to_fits_sip method --- gwcs/wcs.py | 198 ++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 198 insertions(+) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 0119e9ad..dc2cb31b 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -4,6 +4,12 @@ import numpy as np 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, Mapping, + AffineTransformation2D, Pix2Sky_TAN, + RotateNative2Celestial, Identity) +from astropy.modeling.fitting import LinearLSQFitter +import astropy.io.fits as fits from .api import GWCSAPIMixin from . import coordinate_frames @@ -634,3 +640,195 @@ 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, max_error=None, max_rms=None, degree=None): + """ + Make a SIP-based approximation to the WCS to be used in FITS WCS + + Parameters + ---------- + bounding_box: a pair of tuples, each consisting of two integers + Represents the range of pixel values in both dimensions + ((ymin, ymax), (xmin, xmax)) + max_error: float + Maximum allowed error over the domain of the pixel array. + max_rms_error: float + Maximum allowed rms error over the domain of the pixel array. + degree: integer + Degree of the SIP polynomial. If supplied, max_error and max_rms_error + must be None. + + Returns + ------- + fits_dict, max_error, rms_error + + fits_dict is a dictionary where the keys are FITS keywords and the values + the corresonding FITS keyword values, to make it easy to merge with + an existing header. + + This assumes a tangent projection. + """ + + # Handle the options related to degree and thresholds. + if degree is not None and (max_error is not None or max_rms is not None): + raise ValueError( + "Degree cannot be specified if either max_error or max_rms have values") + if degree is None and max_error is None and max_rms is None: + raise ValueError( + "At least one of max_error, max_rms, or degree must be specified") + fits_dict = {} + transform = self.forward_transform + # Determine reference points. + (ymin, ymax), (xmin, xmax) = bounding_box + crpix1 = int((xmax - xmin) / 2) + crpix2 = int((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 + hdr['CRPIX2'] = crpix2 + hdr['CRVAL1'] = crval1 + hdr['CRVAL2'] = crval2 + hdr['cd1_1'] = 0 # Placeholder + hdr['cd1_2'] = 0 + hdr['cd2_1'] = 0 + hdr['cd2_2'] = 0 + # Now rotate to native system and deproject + ntransform = (transform | + RotateCelestial2Native(crval1, crval2, 180) | + Sky2Pix_TAN() ) # | + # (Shift(-crpix1) & Shift(-crpix2))) + # Evaluate this transform on all pixels in bounding box region + y, x = np.mgrid[ymin:ymax, xmin:xmax] # Check on the endpoints needing to be 1 larger + u = x - crpix1 + v = y - crpix2 + undist_x, undist_y = ntransform(x, y) + # The fitting section. + llsqfitter = LinearLSQFitter() + # The case of one pass with the specified polynomial degree + if degree is not None: + poly_x = Polynomial2D(degree=degree) + poly_y = Polynomial2D(degree=degree) + fit_poly_x = llsqfitter(poly_x, u, v, undist_x) + fit_poly_y = llsqfitter(poly_y, u, v, undist_y) + max_resid, rms = compute_distance_residual( + undist_x, undist_y, + fit_poly_x(u, v), fit_poly_y(u, v)) + else: + cond = True + degree = 2 + while cond: + poly_x = Polynomial2D(degree=degree) + poly_y = Polynomial2D(degree=degree) + fit_poly_x = llsqfitter(poly_x, u, v, undist_x) + fit_poly_y = llsqfitter(poly_y, u, v, undist_y) + max_resid, rms = compute_distance_residual( + undist_x, undist_y, + fit_poly_x(u, v), fit_poly_y(u, v)) + print(f'Degree = {degree}, max_resid = {max_resid}, rms = {rms}') + if ((max_error is not None and max_resid < max_error and max_rms is None) or + (max_error is not None and max_resid < max_error and max_rms is not None and rms < max_rms) or + (max_error is None and max_resid is not None and rms > max_rms)): + cond = False + print('terminating condition met') + break + else: + degree +=1 + if degree == 10: + raise RuntimeError("Fit conditions not met") + hdr['a_order'] = degree + hdr['b_order'] = degree + print('fit_poly_x', fit_poly_x) + print('fit_poly_y', fit_poly_y) + cd, sip_poly_x, sip_poly_y = reform_poly_coefficients(fit_poly_x, fit_poly_y) + store_2D_coefficients(hdr, sip_poly_x, 'A') + store_2D_coefficients(hdr, sip_poly_y, 'B') + hdr['cd1_1'] = cd[0][0] + hdr['cd1_2'] = cd[0][1] + hdr['cd2_1'] = cd[1][0] + hdr['cd2_2'] = cd[1][1] + hdr['sipmxerr'] = max_resid + hdr['siprms'] = rms + # # for test purposes, make sip equivalent in GWCS + # cdmat = np.array([[cd[0][0], cd[0][1]], [cd[1][0], cd[1][1]]]) + # # cdmat = [[cd[0][0], cd[1][0]], [cd[0][1], cd[1][1]]] + # sip_poly_test_x = sip_poly_x.copy() + # sip_poly_test_y = sip_poly_y.copy() + # # sip_poly_test_x.c1_0 = 1 + # # sip_poly_test_y.c0_1 = 1 + # shift = (Shift(-crpix1) & Shift(-crpix2)) + # aff = AffineTransformation2D(cdmat, [0,0]) + + # transform = (shift | Mapping([0,1,0,1]) | + # (Mapping([0], n_inputs=2) + sip_poly_test_x) & + # (Mapping([1], n_inputs=2) + sip_poly_test_y) | + # aff | + # Pix2Sky_TAN() | RotateNative2Celestial(crval1, crval2, 180)) + # transform2 = ((Shift(-crpix1) & Shift(-crpix2)) | Mapping([0,1,0,1]) | + # (fit_poly_x & fit_poly_y) | Pix2Sky_TAN() | + # RotateNative2Celestial(crval1, crval2, 180)) + return hdr # transform, transform2 + +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() + rms = np.sqrt((dist**2).sum()/dist.size) + return max_resid, rms + +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 recombing + 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 + det = c11 * c22 - c12 * c21 + ff = c22 / det + fg = c12 / det + gg = c11 / det + gf = c21 / det + 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 + new_x = ff * old_x - fg * old_y + setattr(sip_poly_x, f'c{i}_{j}', new_x) + new_y = gg * old_y - gf * old_x + setattr(sip_poly_y, f'c{i}_{j}', new_y) + cd = ((c11, c12), (c21, c22)) + return cd, sip_poly_x, sip_poly_y + + +def store_2D_coefficients(hdr, poly_model, coeff_prefix): + """ + Write the polynomial model coefficients to the header. + """ + degree = poly_model.degree + for i in range(0, degree + 1): + for j in range(0, degree + 1): + if (i + j) > 1 and (i + j < degree + 1): + hdr[f'{coeff_prefix}_{i}_{j}'] = getattr(poly_model, f'c{i}_{j}').value + From d7638468fee54a6e80b61fcb5cd6dab2b0d37c55 Mon Sep 17 00:00:00 2001 From: perrygreenfield Date: Wed, 26 Feb 2020 13:54:13 -0500 Subject: [PATCH 11/38] added support for SIP inverse generation --- gwcs/tests/test_wcs.py | 19 +++++- gwcs/wcs.py | 140 ++++++++++++++++++++++++++--------------- 2 files changed, 107 insertions(+), 52 deletions(-) diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index bb726f50..ef7cddec 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) @@ -486,3 +486,20 @@ 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]) + warg = np.stack((xflat, yflat)).transpose() + af = asdf.open('data/miriwcs.asdf') + miriwcs = af.tree['wcs'] + bounding_box = ((0, 1024), (0, 1024)) + mirisip = miriwcs.to_fits_sip(bounding_box, max_error=1.e-10, max_inv_error=0.1) + fitssip = astwcs.WCS(mirisip) + fitsvals = fitssip.all_pix2world(warg, 1) + assert_allclose(miriwcs(xflat, yflat), + fitsvals.transpose(), atol=10**-10) + fits_inverse_vals = fitssip.all_world2pix(fitsvals, 1) + assert_allclose(warg, fits_inverse_vals, atol=0.1) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index dc2cb31b..c1a736cf 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -2,6 +2,7 @@ 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, @@ -641,7 +642,8 @@ def fix_inputs(self, fixed): new_pipeline.extend(self.pipeline[1:]) return self.__class__(new_pipeline) - def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None): + def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, + max_inv_error=None, max_inv_rms=None, verbose=False): """ Make a SIP-based approximation to the WCS to be used in FITS WCS @@ -657,6 +659,12 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None): degree: integer Degree of the SIP polynomial. If supplied, max_error and max_rms_error must be None. + max_inv_error: float + Maximum allowed inverse error over the domain of the pixel array. + max_inv_rms: float + Maximum allowed inverse rms over the domain of the pixel array + verbose: boolean + print progress of fits Returns ------- @@ -693,60 +701,44 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None): hdr['CRPIX2'] = crpix2 hdr['CRVAL1'] = crval1 hdr['CRVAL2'] = crval2 - hdr['cd1_1'] = 0 # Placeholder + hdr['cd1_1'] = 0 # Placeholders hdr['cd1_2'] = 0 hdr['cd2_1'] = 0 hdr['cd2_2'] = 0 # Now rotate to native system and deproject - ntransform = (transform | - RotateCelestial2Native(crval1, crval2, 180) | - Sky2Pix_TAN() ) # | - # (Shift(-crpix1) & Shift(-crpix2))) + ntransform = (transform + | RotateCelestial2Native(crval1, crval2, 180) + | Sky2Pix_TAN()) # Evaluate this transform on all pixels in bounding box region - y, x = np.mgrid[ymin:ymax, xmin:xmax] # Check on the endpoints needing to be 1 larger + y, x = np.mgrid[ymin:ymax, xmin:xmax] # Check on the endpoints needing to be 1 larger u = x - crpix1 v = y - crpix2 undist_x, undist_y = ntransform(x, y) # The fitting section. - llsqfitter = LinearLSQFitter() - # The case of one pass with the specified polynomial degree - if degree is not None: - poly_x = Polynomial2D(degree=degree) - poly_y = Polynomial2D(degree=degree) - fit_poly_x = llsqfitter(poly_x, u, v, undist_x) - fit_poly_y = llsqfitter(poly_y, u, v, undist_y) - max_resid, rms = compute_distance_residual( - undist_x, undist_y, - fit_poly_x(u, v), fit_poly_y(u, v)) - else: - cond = True - degree = 2 - while cond: - poly_x = Polynomial2D(degree=degree) - poly_y = Polynomial2D(degree=degree) - fit_poly_x = llsqfitter(poly_x, u, v, undist_x) - fit_poly_y = llsqfitter(poly_y, u, v, undist_y) - max_resid, rms = compute_distance_residual( - undist_x, undist_y, - fit_poly_x(u, v), fit_poly_y(u, v)) - print(f'Degree = {degree}, max_resid = {max_resid}, rms = {rms}') - if ((max_error is not None and max_resid < max_error and max_rms is None) or - (max_error is not None and max_resid < max_error and max_rms is not None and rms < max_rms) or - (max_error is None and max_resid is not None and rms > max_rms)): - cond = False - print('terminating condition met') - break - else: - degree +=1 - if degree == 10: - raise RuntimeError("Fit conditions not met") - hdr['a_order'] = degree - hdr['b_order'] = degree - print('fit_poly_x', fit_poly_x) - print('fit_poly_y', fit_poly_y) + # Fit the forward case. + fit_poly_x, fit_poly_y, max_resid, rms = fit_2D_poly(degree, max_error, max_rms, + u, v, undist_x, undist_y, + verbose=verbose) + 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 + fit_inv_poly_u, fit_inv_poly_v, max_resid, rms = fit_2D_poly(degree, + max_inv_error, max_inv_rms, + U, V, u - U, v - V, + verbose=verbose) + pdegree = fit_poly_x.degree + hdr['a_order'] = pdegree + hdr['b_order'] = pdegree + ipdegree = fit_inv_poly_u.degree + hdr['ap_order'] = ipdegree + hdr['bp_order'] = ipdegree cd, sip_poly_x, sip_poly_y = reform_poly_coefficients(fit_poly_x, fit_poly_y) store_2D_coefficients(hdr, sip_poly_x, 'A') store_2D_coefficients(hdr, sip_poly_y, 'B') + store_2D_coefficients(hdr, fit_inv_poly_u, 'AP', keeplinear=True) + store_2D_coefficients(hdr, fit_inv_poly_v, 'BP', keeplinear=True) hdr['cd1_1'] = cd[0][0] hdr['cd1_2'] = cd[0][1] hdr['cd2_1'] = cd[1][0] @@ -764,24 +756,67 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None): # aff = AffineTransformation2D(cdmat, [0,0]) # transform = (shift | Mapping([0,1,0,1]) | - # (Mapping([0], n_inputs=2) + sip_poly_test_x) & + # (Mapping([0], n_inputs=2) + sip_poly_test_x) & # (Mapping([1], n_inputs=2) + sip_poly_test_y) | # aff | # Pix2Sky_TAN() | RotateNative2Celestial(crval1, crval2, 180)) # transform2 = ((Shift(-crpix1) & Shift(-crpix2)) | Mapping([0,1,0,1]) | - # (fit_poly_x & fit_poly_y) | Pix2Sky_TAN() | + # (fit_poly_x & fit_poly_y) | Pix2Sky_TAN() | # RotateNative2Celestial(crval1, crval2, 180)) - return hdr # transform, transform2 + return hdr + + +def fit_2D_poly(degree, max_error, max_rms, u, v, x, y, verbose=False): + """ + Fit a pair of ordinary 2D polynomial to the supplied inputs (u, v) and + outputs (x, y) + """ + llsqfitter = LinearLSQFitter() + # The case of one pass with the specified polynomial degree + if degree is not None: + poly_x = Polynomial2D(degree=degree) + poly_y = Polynomial2D(degree=degree) + fit_poly_x = llsqfitter(poly_x, u, v, x) + fit_poly_y = llsqfitter(poly_y, u, v, y) + max_resid, rms = compute_distance_residual(x, y, + fit_poly_x(u, v), fit_poly_y(u, v)) + else: + cond = True + degree = 2 + while cond: + poly_x = Polynomial2D(degree=degree) + poly_y = Polynomial2D(degree=degree) + fit_poly_x = llsqfitter(poly_x, u, v, x) + fit_poly_y = llsqfitter(poly_y, u, v, y) + max_resid, rms = compute_distance_residual(x, y, + fit_poly_x(u, v), fit_poly_y(u, v)) + if verbose: + print(f'Degree = {degree}, max_resid = {max_resid}, rms = {rms}') + if ((max_error is not None and max_resid < max_error and max_rms is None) + or (max_error is not None and max_resid < max_error and max_rms is not None + and rms < max_rms) + or (max_error is None and max_rms is not None and rms > max_rms)): + cond = False + if verbose: + print('terminating condition met') + break + else: + degree += 1 + if degree == 10: + raise RuntimeError("Fit conditions not met") + return fit_poly_x, fit_poly_y, max_resid, rms + 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) + dist = np.sqrt((undist_x - fit_poly_x)**2 + (undist_y - fit_poly_y)**2) max_resid = dist.max() - rms = np.sqrt((dist**2).sum()/dist.size) + rms = np.sqrt((dist**2).sum() / dist.size) return max_resid, rms + def reform_poly_coefficients(fit_poly_x, fit_poly_y): """ The fit polynomials must be recombined to align with the SIP decomposition @@ -822,13 +857,16 @@ def reform_poly_coefficients(fit_poly_x, fit_poly_y): return cd, sip_poly_x, sip_poly_y -def store_2D_coefficients(hdr, poly_model, coeff_prefix): +def store_2D_coefficients(hdr, poly_model, coeff_prefix, keeplinear=False): """ Write the polynomial model coefficients to the header. """ + if keeplinear: + mindeg = 0 + else: + mindeg = 1 degree = poly_model.degree for i in range(0, degree + 1): for j in range(0, degree + 1): - if (i + j) > 1 and (i + j < 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 - From 80d5777a9a681f5803014e99294f62f75a82d8fa Mon Sep 17 00:00:00 2001 From: perrygreenfield Date: Wed, 26 Feb 2020 14:03:14 -0500 Subject: [PATCH 12/38] added check to ensure WCS is 2D --- gwcs/wcs.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index c1a736cf..8a2005a7 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -676,6 +676,9 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, This assumes a tangent projection. """ + if (self.forward_transform.n_inputs !=2 + or self.forward_transform.n_outputs != 2): + raise ValueError("The to_fits_sip method only works with 2 dimensional transforms") # Handle the options related to degree and thresholds. if degree is not None and (max_error is not None or max_rms is not None): From 9071ec0a81bbbdcb3c2f3e583783f20fa7198c38 Mon Sep 17 00:00:00 2001 From: perrygreenfield Date: Wed, 26 Feb 2020 14:06:09 -0500 Subject: [PATCH 13/38] updated CHANGES.rst --- CHANGES.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index 5374f607..b2180948 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -6,6 +6,8 @@ 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] + 0.12.0 (2019-12-24) ------------------- From 0adec8889c383b44f204727c88efeca2e2369b59 Mon Sep 17 00:00:00 2001 From: Perry Greenfield Date: Wed, 26 Feb 2020 17:14:38 -0500 Subject: [PATCH 14/38] address most comments, waiting on two, added missed asdf file --- gwcs/tests/data/miriwcs.asdf | 347 +++++++++++++++++++++++++++++++++++ gwcs/tests/test_wcs.py | 4 +- gwcs/wcs.py | 115 +++++------- 3 files changed, 397 insertions(+), 69 deletions(-) create mode 100644 gwcs/tests/data/miriwcs.asdf diff --git a/gwcs/tests/data/miriwcs.asdf b/gwcs/tests/data/miriwcs.asdf new file mode 100644 index 00000000..b74ddfd5 --- /dev/null +++ b/gwcs/tests/data/miriwcs.asdf @@ -0,0 +1,347 @@ +#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.1} +history: + extensions: + - !core/extension_metadata-1.0.0 + extension_class: jwst.transforms.jwextension.JWSTExtension + software: {name: jwst, version: 0.13.8a0.dev236+g8e11a52b} + - !core/extension_metadata-1.0.0 + extension_class: gwcs.extension.GWCSExtension + software: {name: gwcs, version: 0.8.dev460+g1c1cb3b.d20200224} + - !core/extension_metadata-1.0.0 + extension_class: asdf.extension.BuiltinExtension + software: {name: asdf, version: 2.5.1} + - !core/extension_metadata-1.0.0 + extension_class: astropy.io.misc.asdf.extension.AstropyAsdfExtension + software: {name: astropy, version: 4.1.dev580+g3d3d98ea1.d20200212} + - !core/extension_metadata-1.0.0 + extension_class: astropy.io.misc.asdf.extension.AstropyExtension + software: {name: astropy, version: 4.1.dev580+g3d3d98ea1.d20200212} +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: + - [-0.5, 1023.5] + - [3.5, 1027.5] + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/shift-1.2.0 {offset: 0.15} + - !transform/shift-1.2.0 {offset: -0.59} + - !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 {offset: -4.0} + - !transform/identity-1.2.0 {} + - !transform/concatenate-1.2.0 + forward: + - !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 0 + datatype: float64 + byteorder: little + shape: [2] + domain: [-1, 1] + inverse: !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 1 + datatype: float64 + byteorder: little + shape: [2] + domain: [-1, 1] + window: [-1, 1] + name: M_column_correction + window: [-1, 1] + - !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 2 + datatype: float64 + byteorder: little + shape: [2] + domain: [-1, 1] + inverse: !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 3 + datatype: float64 + byteorder: little + shape: [2] + domain: [-1, 1] + window: [-1, 1] + name: M_row_correction + window: [-1, 1] + - !transform/remap_axes-1.2.0 + inverse: !transform/identity-1.2.0 {n_dims: 2} + mapping: [0, 1, 0, 1] + - !transform/concatenate-1.2.0 + forward: + - !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 4 + datatype: float64 + byteorder: little + shape: [5, 5] + domain: + - [-1, 1] + - [-1, 1] + inverse: !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 5 + datatype: float64 + byteorder: little + shape: [5, 5] + domain: + - [-1, 1] + - [-1, 1] + window: + - [-1, 1] + - [-1, 1] + name: B_correction + window: + - [-1, 1] + - [-1, 1] + - !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 6 + datatype: float64 + byteorder: little + shape: [5, 5] + domain: + - [-1, 1] + - [-1, 1] + inverse: !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 7 + datatype: float64 + byteorder: little + shape: [5, 5] + domain: + - [-1, 1] + - [-1, 1] + window: + - [-1, 1] + - [-1, 1] + name: A_correction + window: + - [-1, 1] + - [-1, 1] + - !transform/remap_axes-1.2.0 + inverse: !transform/remap_axes-1.2.0 + mapping: [0, 1, 0, 1] + mapping: [0, 1, 0, 1] + - !transform/concatenate-1.2.0 + forward: + - !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 8 + datatype: float64 + byteorder: little + shape: [2, 2] + domain: + - [-1, 1] + - [-1, 1] + name: TI_row_correction + window: + - [-1, 1] + - [-1, 1] + - !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 9 + datatype: float64 + byteorder: little + shape: [2, 2] + domain: + - [-1, 1] + - [-1, 1] + name: TI_column_correction + window: + - [-1, 1] + - [-1, 1] + - !transform/identity-1.2.0 + inverse: !transform/remap_axes-1.2.0 + mapping: [0, 1, 0, 1] + n_dims: 2 + - !transform/remap_axes-1.2.0 + mapping: [1, 0] + inverse: !transform/compose-1.2.0 + forward: + - !transform/compose-1.2.0 + forward: + - !transform/remap_axes-1.2.0 + mapping: [1, 0] + - !transform/compose-1.2.0 + forward: + - !transform/remap_axes-1.2.0 + mapping: [0, 1, 0, 1] + - !transform/compose-1.2.0 + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 10 + datatype: float64 + byteorder: little + shape: [2, 2] + domain: + - [-1, 1] + - [-1, 1] + name: T_row_correction + window: + - [-1, 1] + - [-1, 1] + - !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 11 + datatype: float64 + byteorder: little + shape: [2, 2] + domain: + - [-1, 1] + - [-1, 1] + name: T_column_correction + window: + - [-1, 1] + - [-1, 1] + - !transform/compose-1.2.0 + forward: + - !transform/remap_axes-1.2.0 + mapping: [0, 1, 0, 1] + - !transform/compose-1.2.0 + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 12 + datatype: float64 + byteorder: little + shape: [5, 5] + domain: + - [-1, 1] + - [-1, 1] + window: + - [-1, 1] + - [-1, 1] + - !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 13 + datatype: float64 + byteorder: little + shape: [5, 5] + domain: + - [-1, 1] + - [-1, 1] + window: + - [-1, 1] + - [-1, 1] + - !transform/compose-1.2.0 + forward: + - !transform/identity-1.2.0 {n_dims: 2} + - !transform/compose-1.2.0 + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 14 + datatype: float64 + byteorder: little + shape: [2] + domain: [-1, 1] + window: [-1, 1] + - !transform/polynomial-1.2.0 + coefficients: !core/ndarray-1.0.0 + source: 15 + datatype: float64 + byteorder: little + shape: [2] + domain: [-1, 1] + window: [-1, 1] + - !transform/concatenate-1.2.0 + forward: + - !transform/shift-1.2.0 {offset: 4.0} + - !transform/identity-1.2.0 {} + - !transform/concatenate-1.2.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.2.0 + forward: + - !transform/concatenate-1.2.0 + forward: + - !transform/scale-1.2.0 {factor: 0.0002777777777777778} + - !transform/scale-1.2.0 {factor: 0.0002777777777777778} + - ! + angles: [0.12597594444444443, -0.10374517305555556, -0.0, 72.0545718, 5.630568] + axes_order: zyxyz + inputs: [v2, v3] + name: v23tosky + outputs: [ra, dec] + - ! + 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 +--- +- 12376 +- 12446 +- 12516 +- 12586 +- 12656 +- 12910 +- 13164 +- 13418 +- 13672 +- 13758 +- 13844 +- 13930 +- 14016 +- 14270 +- 14524 +- 14594 +... diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index ef7cddec..81205acf 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -500,6 +500,6 @@ def test_to_fits_sip(): fitssip = astwcs.WCS(mirisip) fitsvals = fitssip.all_pix2world(warg, 1) assert_allclose(miriwcs(xflat, yflat), - fitsvals.transpose(), atol=10**-10) + fitsvals.transpose(), atol=1e-10, rtol=0) fits_inverse_vals = fitssip.all_world2pix(fitsvals, 1) - assert_allclose(warg, fits_inverse_vals, atol=0.1) + assert_allclose(warg, fits_inverse_vals, atol=0.1, rtol=0) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 8a2005a7..be5ab83a 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -668,13 +668,24 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, Returns ------- - fits_dict, max_error, rms_error - - fits_dict is a dictionary where the keys are FITS keywords and the values - the corresonding FITS keyword values, to make it easy to merge with - an existing header. + FITS header with all SIP WCS keywords This assumes a tangent projection. + + Error cases + ----------- + + If the WCS is not 2D, an exception will be raised. If the specified accuracy + (both forward and inverse) 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 (self.forward_transform.n_inputs !=2 or self.forward_transform.n_outputs != 2): @@ -719,7 +730,7 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, undist_x, undist_y = ntransform(x, y) # The fitting section. # Fit the forward case. - fit_poly_x, fit_poly_y, max_resid, rms = fit_2D_poly(degree, max_error, max_rms, + fit_poly_x, fit_poly_y, max_resid, rms = _fit_2D_poly(degree, max_error, max_rms, u, v, undist_x, undist_y, verbose=verbose) cdmat = np.array([[fit_poly_x.c1_0.value, fit_poly_x.c0_1.value], @@ -727,7 +738,7 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, 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 - fit_inv_poly_u, fit_inv_poly_v, max_resid, rms = fit_2D_poly(degree, + fit_inv_poly_u, fit_inv_poly_v, max_resid, rms = _fit_2D_poly(degree, max_inv_error, max_inv_rms, U, V, u - U, v - V, verbose=verbose) @@ -737,80 +748,53 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, ipdegree = fit_inv_poly_u.degree hdr['ap_order'] = ipdegree hdr['bp_order'] = ipdegree - cd, sip_poly_x, sip_poly_y = reform_poly_coefficients(fit_poly_x, fit_poly_y) - store_2D_coefficients(hdr, sip_poly_x, 'A') - store_2D_coefficients(hdr, sip_poly_y, 'B') - store_2D_coefficients(hdr, fit_inv_poly_u, 'AP', keeplinear=True) - store_2D_coefficients(hdr, fit_inv_poly_v, 'BP', keeplinear=True) + cd, sip_poly_x, sip_poly_y = _reform_poly_coefficients(fit_poly_x, fit_poly_y) + _store_2D_coefficients(hdr, sip_poly_x, 'A') + _store_2D_coefficients(hdr, sip_poly_y, 'B') + _store_2D_coefficients(hdr, fit_inv_poly_u, 'AP', keeplinear=True) + _store_2D_coefficients(hdr, fit_inv_poly_v, 'BP', keeplinear=True) hdr['cd1_1'] = cd[0][0] hdr['cd1_2'] = cd[0][1] hdr['cd2_1'] = cd[1][0] hdr['cd2_2'] = cd[1][1] hdr['sipmxerr'] = max_resid hdr['siprms'] = rms - # # for test purposes, make sip equivalent in GWCS - # cdmat = np.array([[cd[0][0], cd[0][1]], [cd[1][0], cd[1][1]]]) - # # cdmat = [[cd[0][0], cd[1][0]], [cd[0][1], cd[1][1]]] - # sip_poly_test_x = sip_poly_x.copy() - # sip_poly_test_y = sip_poly_y.copy() - # # sip_poly_test_x.c1_0 = 1 - # # sip_poly_test_y.c0_1 = 1 - # shift = (Shift(-crpix1) & Shift(-crpix2)) - # aff = AffineTransformation2D(cdmat, [0,0]) - - # transform = (shift | Mapping([0,1,0,1]) | - # (Mapping([0], n_inputs=2) + sip_poly_test_x) & - # (Mapping([1], n_inputs=2) + sip_poly_test_y) | - # aff | - # Pix2Sky_TAN() | RotateNative2Celestial(crval1, crval2, 180)) - # transform2 = ((Shift(-crpix1) & Shift(-crpix2)) | Mapping([0,1,0,1]) | - # (fit_poly_x & fit_poly_y) | Pix2Sky_TAN() | - # RotateNative2Celestial(crval1, crval2, 180)) return hdr -def fit_2D_poly(degree, max_error, max_rms, u, v, x, y, verbose=False): +def _fit_2D_poly(degree, max_error, max_rms, u, v, x, y, verbose=False): """ Fit a pair of ordinary 2D polynomial to the supplied inputs (u, v) and outputs (x, y) """ llsqfitter = LinearLSQFitter() # The case of one pass with the specified polynomial degree - if degree is not None: - poly_x = Polynomial2D(degree=degree) - poly_y = Polynomial2D(degree=degree) + if degree is None: + deglist = range(2,10) + else: + deglist = [degree] + for deg in deglist: + poly_x = Polynomial2D(degree=deg) + poly_y = Polynomial2D(degree=deg) fit_poly_x = llsqfitter(poly_x, u, v, x) fit_poly_y = llsqfitter(poly_y, u, v, y) - max_resid, rms = compute_distance_residual(x, y, - fit_poly_x(u, v), fit_poly_y(u, v)) - else: - cond = True - degree = 2 - while cond: - poly_x = Polynomial2D(degree=degree) - poly_y = Polynomial2D(degree=degree) - fit_poly_x = llsqfitter(poly_x, u, v, x) - fit_poly_y = llsqfitter(poly_y, u, v, y) - max_resid, rms = compute_distance_residual(x, y, - fit_poly_x(u, v), fit_poly_y(u, v)) + max_resid, rms = _compute_distance_residual(x, y, + fit_poly_x(u, v), fit_poly_y(u, v)) + if verbose: + print(f'Degree = {degree}, max_resid = {max_resid}, rms = {rms}') + if ((max_error is not None and max_resid < max_error and max_rms is None) + or (max_error is not None and max_resid < max_error and max_rms is not None + and rms < max_rms) + or (max_error is None and max_rms is not None and rms > max_rms)): if verbose: - print(f'Degree = {degree}, max_resid = {max_resid}, rms = {rms}') - if ((max_error is not None and max_resid < max_error and max_rms is None) - or (max_error is not None and max_resid < max_error and max_rms is not None - and rms < max_rms) - or (max_error is None and max_rms is not None and rms > max_rms)): - cond = False - if verbose: - print('terminating condition met') - break - else: - degree += 1 - if degree == 10: - raise RuntimeError("Fit conditions not met") + print('terminating condition met') + break + if deg == 9: + raise RuntimeError("Unable to achieve required fit accuracy with SIP polynomials") return fit_poly_x, fit_poly_y, max_resid, rms -def compute_distance_residual(undist_x, undist_y, fit_poly_x, fit_poly_y): +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. """ @@ -820,13 +804,13 @@ def compute_distance_residual(undist_x, undist_y, fit_poly_x, fit_poly_y): return max_resid, rms -def reform_poly_coefficients(fit_poly_x, fit_poly_y): +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 recombing + # 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 @@ -860,14 +844,11 @@ def reform_poly_coefficients(fit_poly_x, fit_poly_y): return cd, sip_poly_x, sip_poly_y -def store_2D_coefficients(hdr, poly_model, coeff_prefix, keeplinear=False): +def _store_2D_coefficients(hdr, poly_model, coeff_prefix, keeplinear=False): """ Write the polynomial model coefficients to the header. """ - if keeplinear: - mindeg = 0 - else: - mindeg = 1 + mindeg = int(not keeplinear) degree = poly_model.degree for i in range(0, degree + 1): for j in range(0, degree + 1): From d3612d319bf43f4e44751f7f2dee5ea16ac9cfa7 Mon Sep 17 00:00:00 2001 From: Perry Greenfield Date: Thu, 27 Feb 2020 04:29:09 -0500 Subject: [PATCH 15/38] fix origin issues --- gwcs/tests/test_wcs.py | 4 ++-- gwcs/wcs.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 81205acf..a5a4c3c6 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -498,8 +498,8 @@ def test_to_fits_sip(): bounding_box = ((0, 1024), (0, 1024)) mirisip = miriwcs.to_fits_sip(bounding_box, max_error=1.e-10, max_inv_error=0.1) fitssip = astwcs.WCS(mirisip) - fitsvals = fitssip.all_pix2world(warg, 1) + fitsvals = fitssip.all_pix2world(warg + 1, 1) assert_allclose(miriwcs(xflat, yflat), fitsvals.transpose(), atol=1e-10, rtol=0) fits_inverse_vals = fitssip.all_world2pix(fitsvals, 1) - assert_allclose(warg, fits_inverse_vals, atol=0.1, rtol=0) + assert_allclose(warg, fits_inverse_vals - 1, atol=0.1, rtol=0) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index be5ab83a..4eae7778 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -711,8 +711,8 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, hdr['naxis2'] = ymax hdr['ctype1'] = 'RA---TAN-SIP' hdr['ctype2'] = 'DEC--TAN-SIP' - hdr['CRPIX1'] = crpix1 - hdr['CRPIX2'] = crpix2 + hdr['CRPIX1'] = crpix1 + 1 + hdr['CRPIX2'] = crpix2 + 1 hdr['CRVAL1'] = crval1 hdr['CRVAL2'] = crval2 hdr['cd1_1'] = 0 # Placeholders From 599b63a32e68e0e67fc9c975e193dfff14116e48 Mon Sep 17 00:00:00 2001 From: Perry Greenfield Date: Thu, 27 Feb 2020 05:12:50 -0500 Subject: [PATCH 16/38] address further comments, add support for linear case --- gwcs/tests/test_wcs.py | 2 +- gwcs/wcs.py | 39 +++++++++++++++++++++++---------------- 2 files changed, 24 insertions(+), 17 deletions(-) diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index a5a4c3c6..85aca436 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -493,7 +493,7 @@ def test_to_fits_sip(): xflat = np.ravel(x[1:-1, 1:-1]) yflat = np.ravel(y[1:-1, 1:-1]) warg = np.stack((xflat, yflat)).transpose() - af = asdf.open('data/miriwcs.asdf') + 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_error=1.e-10, max_inv_error=0.1) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 4eae7778..7f29645e 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -672,11 +672,12 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, This assumes a tangent projection. - Error cases - ----------- + Raises + ------ If the WCS is not 2D, an exception will be raised. If the specified accuracy - (both forward and inverse) is not achieved an exception will be raised. + (both forward and inverse, both rms and maximum) is not achieved an exception + will be raised. Notes ----- @@ -743,22 +744,28 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, U, V, u - U, v - V, verbose=verbose) pdegree = fit_poly_x.degree - hdr['a_order'] = pdegree - hdr['b_order'] = pdegree - ipdegree = fit_inv_poly_u.degree - hdr['ap_order'] = ipdegree - hdr['bp_order'] = ipdegree - cd, sip_poly_x, sip_poly_y = _reform_poly_coefficients(fit_poly_x, fit_poly_y) - _store_2D_coefficients(hdr, sip_poly_x, 'A') - _store_2D_coefficients(hdr, sip_poly_y, 'B') - _store_2D_coefficients(hdr, fit_inv_poly_u, 'AP', keeplinear=True) - _store_2D_coefficients(hdr, fit_inv_poly_v, 'BP', keeplinear=True) + if pdegree > 1: + hdr['a_order'] = pdegree + hdr['b_order'] = pdegree + ipdegree = fit_inv_poly_u.degree + hdr['ap_order'] = ipdegree + hdr['bp_order'] = ipdegree + cd, sip_poly_x, sip_poly_y = _reform_poly_coefficients(fit_poly_x, fit_poly_y) + _store_2D_coefficients(hdr, sip_poly_x, 'A') + _store_2D_coefficients(hdr, sip_poly_y, 'B') + _store_2D_coefficients(hdr, fit_inv_poly_u, 'AP', keeplinear=True) + _store_2D_coefficients(hdr, fit_inv_poly_v, 'BP', keeplinear=True) + hdr['sipmxerr'] = (max_resid, 'Maximum difference from the GWCS model SIP is fit to.') + hdr['siprms'] = (rms, 'RMS difference from the GWCS model SIP is fit to.') + + else: + hdr['ctype1'] = 'RA---TAN' + hdr['ctype2'] = 'DEC--TAN' + hdr['cd1_1'] = cd[0][0] hdr['cd1_2'] = cd[0][1] hdr['cd2_1'] = cd[1][0] hdr['cd2_2'] = cd[1][1] - hdr['sipmxerr'] = max_resid - hdr['siprms'] = rms return hdr @@ -770,7 +777,7 @@ def _fit_2D_poly(degree, max_error, max_rms, u, v, x, y, verbose=False): llsqfitter = LinearLSQFitter() # The case of one pass with the specified polynomial degree if degree is None: - deglist = range(2,10) + deglist = range(1,10) else: deglist = [degree] for deg in deglist: From 40916005936cbe82b0302950b4a3c75b0dba62c0 Mon Sep 17 00:00:00 2001 From: Perry Greenfield Date: Thu, 27 Feb 2020 16:34:45 -0500 Subject: [PATCH 17/38] address further comments from @nden --- gwcs/wcs.py | 47 ++++++++++++++++++++++------------------------- 1 file changed, 22 insertions(+), 25 deletions(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 7f29645e..d2f73d61 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -645,39 +645,39 @@ def fix_inputs(self, fixed): def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, max_inv_error=None, max_inv_rms=None, verbose=False): """ - Make a SIP-based approximation to the WCS to be used in FITS WCS + Construct a SIP-based approximation to the WCS in the form of a FITS header + + This assumes a tangent projection. Parameters ---------- - bounding_box: a pair of tuples, each consisting of two integers + bounding_box : a pair of tuples, each consisting of two integers Represents the range of pixel values in both dimensions ((ymin, ymax), (xmin, xmax)) - max_error: float + max_error : float Maximum allowed error over the domain of the pixel array. - max_rms_error: float + max_rms_error : float Maximum allowed rms error over the domain of the pixel array. - degree: integer + degree : int Degree of the SIP polynomial. If supplied, max_error and max_rms_error must be None. - max_inv_error: float + max_inv_error : float Maximum allowed inverse error over the domain of the pixel array. - max_inv_rms: float + max_inv_rms : float Maximum allowed inverse rms over the domain of the pixel array - verbose: boolean + verbose : bool print progress of fits Returns ------- FITS header with all SIP WCS keywords - This assumes a tangent projection. - Raises ------ - - 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. + 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 ----- @@ -831,11 +831,9 @@ def _reform_poly_coefficients(fit_poly_x, fit_poly_y): sip_poly_x.c0_1 = 0 sip_poly_y.c1_0 = 0 sip_poly_y.c0_1 = 0 - det = c11 * c22 - c12 * c21 - ff = c22 / det - fg = c12 / det - gg = c11 / det - gf = c21 / det + + 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): @@ -843,12 +841,11 @@ def _reform_poly_coefficients(fit_poly_x, fit_poly_y): 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 - new_x = ff * old_x - fg * old_y - setattr(sip_poly_x, f'c{i}_{j}', new_x) - new_y = gg * old_y - gf * old_x - setattr(sip_poly_y, f'c{i}_{j}', new_y) - cd = ((c11, c12), (c21, c22)) - return cd, sip_poly_x, sip_poly_y + 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): From b6d3221e884581b79eec75be7ce5cadb4b287a2a Mon Sep 17 00:00:00 2001 From: Perry Greenfield Date: Thu, 27 Feb 2020 16:47:08 -0500 Subject: [PATCH 18/38] ripped out rms threshold options --- gwcs/wcs.py | 42 +++++++++++++++++------------------------- 1 file changed, 17 insertions(+), 25 deletions(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index d2f73d61..9ebbba86 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -656,15 +656,10 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, ((ymin, ymax), (xmin, xmax)) max_error : float Maximum allowed error over the domain of the pixel array. - max_rms_error : float - Maximum allowed rms error over the domain of the pixel array. degree : int - Degree of the SIP polynomial. If supplied, max_error and max_rms_error - must be None. + Degree of the SIP polynomial. If supplied, max_error must be None. max_inv_error : float Maximum allowed inverse error over the domain of the pixel array. - max_inv_rms : float - Maximum allowed inverse rms over the domain of the pixel array verbose : bool print progress of fits @@ -690,15 +685,15 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, """ if (self.forward_transform.n_inputs !=2 or self.forward_transform.n_outputs != 2): - raise ValueError("The to_fits_sip method only works with 2 dimensional transforms") + raise ValueError("The to_fits_sip method only works with 2-dimensional transforms") # Handle the options related to degree and thresholds. - if degree is not None and (max_error is not None or max_rms is not None): + if degree is not None and max_error is not None: raise ValueError( - "Degree cannot be specified if either max_error or max_rms have values") + "Degree cannot be specified if max_error has a value") if degree is None and max_error is None and max_rms is None: raise ValueError( - "At least one of max_error, max_rms, or degree must be specified") + "max_error must be specified if degree is None") fits_dict = {} transform = self.forward_transform # Determine reference points. @@ -731,7 +726,7 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, undist_x, undist_y = ntransform(x, y) # The fitting section. # Fit the forward case. - fit_poly_x, fit_poly_y, max_resid, rms = _fit_2D_poly(degree, max_error, max_rms, + fit_poly_x, fit_poly_y, max_resid = _fit_2D_poly(degree, max_error, u, v, undist_x, undist_y, verbose=verbose) cdmat = np.array([[fit_poly_x.c1_0.value, fit_poly_x.c0_1.value], @@ -739,8 +734,8 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, 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 - fit_inv_poly_u, fit_inv_poly_v, max_resid, rms = _fit_2D_poly(degree, - max_inv_error, max_inv_rms, + fit_inv_poly_u, fit_inv_poly_v, max_inv_resid = _fit_2D_poly(degree, + max_inv_error, U, V, u - U, v - V, verbose=verbose) pdegree = fit_poly_x.degree @@ -756,7 +751,8 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, _store_2D_coefficients(hdr, fit_inv_poly_u, 'AP', keeplinear=True) _store_2D_coefficients(hdr, fit_inv_poly_v, 'BP', keeplinear=True) hdr['sipmxerr'] = (max_resid, 'Maximum difference from the GWCS model SIP is fit to.') - hdr['siprms'] = (rms, 'RMS difference from the GWCS model SIP is fit to.') + hdr['sipiverr'] = (max_inv_resid, 'Maximum difference for the inverse transform') + else: hdr['ctype1'] = 'RA---TAN' @@ -769,7 +765,7 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, return hdr -def _fit_2D_poly(degree, max_error, max_rms, u, v, x, y, verbose=False): +def _fit_2D_poly(degree, max_error, u, v, x, y, verbose=False): """ Fit a pair of ordinary 2D polynomial to the supplied inputs (u, v) and outputs (x, y) @@ -785,20 +781,17 @@ def _fit_2D_poly(degree, max_error, max_rms, u, v, x, y, verbose=False): poly_y = Polynomial2D(degree=deg) fit_poly_x = llsqfitter(poly_x, u, v, x) fit_poly_y = llsqfitter(poly_y, u, v, y) - max_resid, rms = _compute_distance_residual(x, y, - fit_poly_x(u, v), fit_poly_y(u, v)) + max_resid = _compute_distance_residual(x, y, + fit_poly_x(u, v), fit_poly_y(u, v)) if verbose: - print(f'Degree = {degree}, max_resid = {max_resid}, rms = {rms}') - if ((max_error is not None and max_resid < max_error and max_rms is None) - or (max_error is not None and max_resid < max_error and max_rms is not None - and rms < max_rms) - or (max_error is None and max_rms is not None and rms > max_rms)): + print(f'Degree = {degree}, max_resid = {max_resid}') + if max_error is not None and max_resid < max_error: if verbose: print('terminating condition met') break if deg == 9: raise RuntimeError("Unable to achieve required fit accuracy with SIP polynomials") - return fit_poly_x, fit_poly_y, max_resid, rms + return fit_poly_x, fit_poly_y, max_resid def _compute_distance_residual(undist_x, undist_y, fit_poly_x, fit_poly_y): @@ -807,8 +800,7 @@ def _compute_distance_residual(undist_x, undist_y, fit_poly_x, fit_poly_y): """ dist = np.sqrt((undist_x - fit_poly_x)**2 + (undist_y - fit_poly_y)**2) max_resid = dist.max() - rms = np.sqrt((dist**2).sum() / dist.size) - return max_resid, rms + return max_resid def _reform_poly_coefficients(fit_poly_x, fit_poly_y): From f8415b2c9a26c9535662270cf019bf757381fa00 Mon Sep 17 00:00:00 2001 From: Perry Greenfield Date: Thu, 27 Feb 2020 16:54:28 -0500 Subject: [PATCH 19/38] reordered bounding box elements to be x, y order --- gwcs/wcs.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 9ebbba86..5e5f2376 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -653,7 +653,7 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, ---------- bounding_box : a pair of tuples, each consisting of two integers Represents the range of pixel values in both dimensions - ((ymin, ymax), (xmin, xmax)) + ((xmin, xmax), (ymin, ymax)) max_error : float Maximum allowed error over the domain of the pixel array. degree : int @@ -697,7 +697,7 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, fits_dict = {} transform = self.forward_transform # Determine reference points. - (ymin, ymax), (xmin, xmax) = bounding_box + (xmin, xmax), (ymin, ymax) = bounding_box crpix1 = int((xmax - xmin) / 2) crpix2 = int((ymax - ymin) / 2) crval1, crval2 = transform(crpix1, crpix2) From d2521abc2d46789778e9ce950ea64f0e007046ff Mon Sep 17 00:00:00 2001 From: Perry Greenfield Date: Thu, 27 Feb 2020 16:57:01 -0500 Subject: [PATCH 20/38] removed extra line --- gwcs/wcs.py | 1 - 1 file changed, 1 deletion(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 5e5f2376..b9796e59 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -753,7 +753,6 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, hdr['sipmxerr'] = (max_resid, 'Maximum difference from the GWCS model SIP is fit to.') hdr['sipiverr'] = (max_inv_resid, 'Maximum difference for the inverse transform') - else: hdr['ctype1'] = 'RA---TAN' hdr['ctype2'] = 'DEC--TAN' From 60e12c6f85792000d15bc85a879121d37faff481 Mon Sep 17 00:00:00 2001 From: Perry Greenfield Date: Thu, 27 Feb 2020 17:12:06 -0500 Subject: [PATCH 21/38] fix bug in fit termination test --- gwcs/wcs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index b9796e59..35262506 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -788,7 +788,7 @@ def _fit_2D_poly(degree, max_error, u, v, x, y, verbose=False): if verbose: print('terminating condition met') break - if deg == 9: + if degree is None and deg == 9: raise RuntimeError("Unable to achieve required fit accuracy with SIP polynomials") return fit_poly_x, fit_poly_y, max_resid From 23adbd19310f16ec23f467f52a34095e1e49cf44 Mon Sep 17 00:00:00 2001 From: Perry Greenfield Date: Wed, 4 Mar 2020 12:03:23 -0500 Subject: [PATCH 22/38] add automatic fitting mode --- gwcs/tests/data/miriwcs.asdf | 231 +++++++++++------------------------ gwcs/tests/test_wcs.py | 2 +- gwcs/wcs.py | 147 ++++++++++++++++------ 3 files changed, 182 insertions(+), 198 deletions(-) diff --git a/gwcs/tests/data/miriwcs.asdf b/gwcs/tests/data/miriwcs.asdf index b74ddfd5..bde9439b 100644 --- a/gwcs/tests/data/miriwcs.asdf +++ b/gwcs/tests/data/miriwcs.asdf @@ -7,21 +7,12 @@ asdf_library: !core/software-1.0.0 {author: Space Telescope Science Institute, h name: asdf, version: 2.5.1} history: extensions: - - !core/extension_metadata-1.0.0 - extension_class: jwst.transforms.jwextension.JWSTExtension - software: {name: jwst, version: 0.13.8a0.dev236+g8e11a52b} - !core/extension_metadata-1.0.0 extension_class: gwcs.extension.GWCSExtension - software: {name: gwcs, version: 0.8.dev460+g1c1cb3b.d20200224} + software: {name: gwcs, version: 0.12.0} - !core/extension_metadata-1.0.0 extension_class: asdf.extension.BuiltinExtension software: {name: asdf, version: 2.5.1} - - !core/extension_metadata-1.0.0 - extension_class: astropy.io.misc.asdf.extension.AstropyAsdfExtension - software: {name: astropy, version: 4.1.dev580+g3d3d98ea1.d20200212} - - !core/extension_metadata-1.0.0 - extension_class: astropy.io.misc.asdf.extension.AstropyExtension - software: {name: astropy, version: 4.1.dev580+g3d3d98ea1.d20200212} wcs: ! name: '' steps: @@ -32,262 +23,190 @@ wcs: ! 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 + transform: !transform/compose-1.1.0 bounding_box: - [-0.5, 1023.5] - [3.5, 1027.5] forward: - - !transform/concatenate-1.2.0 + - !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.2.0 + - !transform/compose-1.1.0 forward: - - !transform/compose-1.2.0 + - !transform/compose-1.1.0 forward: - - !transform/compose-1.2.0 + - !transform/compose-1.1.0 forward: - - !transform/compose-1.2.0 + - !transform/compose-1.1.0 forward: - - !transform/compose-1.2.0 + - !transform/compose-1.1.0 forward: - - !transform/compose-1.2.0 + - !transform/compose-1.1.0 forward: - - !transform/compose-1.2.0 + - !transform/compose-1.1.0 forward: - - !transform/concatenate-1.2.0 + - !transform/concatenate-1.1.0 forward: - !transform/shift-1.2.0 {offset: -4.0} - - !transform/identity-1.2.0 {} - - !transform/concatenate-1.2.0 + - !transform/identity-1.1.0 {} + - !transform/concatenate-1.1.0 forward: - - !transform/polynomial-1.2.0 + - !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 0 datatype: float64 byteorder: little shape: [2] - domain: [-1, 1] - inverse: !transform/polynomial-1.2.0 + inverse: !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 1 datatype: float64 byteorder: little shape: [2] - domain: [-1, 1] - window: [-1, 1] name: M_column_correction - window: [-1, 1] - - !transform/polynomial-1.2.0 + - !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 2 datatype: float64 byteorder: little shape: [2] - domain: [-1, 1] - inverse: !transform/polynomial-1.2.0 + inverse: !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 3 datatype: float64 byteorder: little shape: [2] - domain: [-1, 1] - window: [-1, 1] name: M_row_correction - window: [-1, 1] - - !transform/remap_axes-1.2.0 - inverse: !transform/identity-1.2.0 {n_dims: 2} + - !transform/remap_axes-1.1.0 + inverse: !transform/identity-1.1.0 {n_dims: 2} mapping: [0, 1, 0, 1] - - !transform/concatenate-1.2.0 + - !transform/concatenate-1.1.0 forward: - - !transform/polynomial-1.2.0 + - !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 4 datatype: float64 byteorder: little shape: [5, 5] - domain: - - [-1, 1] - - [-1, 1] - inverse: !transform/polynomial-1.2.0 + inverse: !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 5 datatype: float64 byteorder: little shape: [5, 5] - domain: - - [-1, 1] - - [-1, 1] - window: - - [-1, 1] - - [-1, 1] name: B_correction - window: - - [-1, 1] - - [-1, 1] - - !transform/polynomial-1.2.0 + - !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 6 datatype: float64 byteorder: little shape: [5, 5] - domain: - - [-1, 1] - - [-1, 1] - inverse: !transform/polynomial-1.2.0 + inverse: !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 7 datatype: float64 byteorder: little shape: [5, 5] - domain: - - [-1, 1] - - [-1, 1] - window: - - [-1, 1] - - [-1, 1] name: A_correction - window: - - [-1, 1] - - [-1, 1] - - !transform/remap_axes-1.2.0 - inverse: !transform/remap_axes-1.2.0 + - !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.2.0 + - !transform/concatenate-1.1.0 forward: - - !transform/polynomial-1.2.0 + - !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 8 datatype: float64 byteorder: little shape: [2, 2] - domain: - - [-1, 1] - - [-1, 1] name: TI_row_correction - window: - - [-1, 1] - - [-1, 1] - - !transform/polynomial-1.2.0 + - !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 9 datatype: float64 byteorder: little shape: [2, 2] - domain: - - [-1, 1] - - [-1, 1] name: TI_column_correction - window: - - [-1, 1] - - [-1, 1] - - !transform/identity-1.2.0 - inverse: !transform/remap_axes-1.2.0 + - !transform/identity-1.1.0 + inverse: !transform/remap_axes-1.1.0 mapping: [0, 1, 0, 1] n_dims: 2 - - !transform/remap_axes-1.2.0 + - !transform/remap_axes-1.1.0 mapping: [1, 0] - inverse: !transform/compose-1.2.0 + inverse: !transform/compose-1.1.0 forward: - - !transform/compose-1.2.0 + - !transform/compose-1.1.0 forward: - - !transform/remap_axes-1.2.0 + - !transform/remap_axes-1.1.0 mapping: [1, 0] - - !transform/compose-1.2.0 + - !transform/compose-1.1.0 forward: - - !transform/remap_axes-1.2.0 + - !transform/remap_axes-1.1.0 mapping: [0, 1, 0, 1] - - !transform/compose-1.2.0 + - !transform/compose-1.1.0 forward: - - !transform/concatenate-1.2.0 + - !transform/concatenate-1.1.0 forward: - - !transform/polynomial-1.2.0 + - !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 10 datatype: float64 byteorder: little shape: [2, 2] - domain: - - [-1, 1] - - [-1, 1] name: T_row_correction - window: - - [-1, 1] - - [-1, 1] - - !transform/polynomial-1.2.0 + - !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 11 datatype: float64 byteorder: little shape: [2, 2] - domain: - - [-1, 1] - - [-1, 1] name: T_column_correction - window: - - [-1, 1] - - [-1, 1] - - !transform/compose-1.2.0 + - !transform/compose-1.1.0 forward: - - !transform/remap_axes-1.2.0 + - !transform/remap_axes-1.1.0 mapping: [0, 1, 0, 1] - - !transform/compose-1.2.0 + - !transform/compose-1.1.0 forward: - - !transform/concatenate-1.2.0 + - !transform/concatenate-1.1.0 forward: - - !transform/polynomial-1.2.0 + - !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 12 datatype: float64 byteorder: little shape: [5, 5] - domain: - - [-1, 1] - - [-1, 1] - window: - - [-1, 1] - - [-1, 1] - - !transform/polynomial-1.2.0 + - !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 13 datatype: float64 byteorder: little shape: [5, 5] - domain: - - [-1, 1] - - [-1, 1] - window: - - [-1, 1] - - [-1, 1] - - !transform/compose-1.2.0 + - !transform/compose-1.1.0 forward: - - !transform/identity-1.2.0 {n_dims: 2} - - !transform/compose-1.2.0 + - !transform/identity-1.1.0 {n_dims: 2} + - !transform/compose-1.1.0 forward: - - !transform/concatenate-1.2.0 + - !transform/concatenate-1.1.0 forward: - - !transform/polynomial-1.2.0 + - !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 14 datatype: float64 byteorder: little shape: [2] - domain: [-1, 1] - window: [-1, 1] - - !transform/polynomial-1.2.0 + - !transform/polynomial-1.1.0 coefficients: !core/ndarray-1.0.0 source: 15 datatype: float64 byteorder: little shape: [2] - domain: [-1, 1] - window: [-1, 1] - - !transform/concatenate-1.2.0 + - !transform/concatenate-1.1.0 forward: - !transform/shift-1.2.0 {offset: 4.0} - - !transform/identity-1.2.0 {} - - !transform/concatenate-1.2.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} @@ -298,18 +217,16 @@ wcs: ! 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.2.0 + transform: !transform/compose-1.1.0 forward: - - !transform/concatenate-1.2.0 + - !transform/concatenate-1.1.0 forward: - !transform/scale-1.2.0 {factor: 0.0002777777777777778} - !transform/scale-1.2.0 {factor: 0.0002777777777777778} - ! angles: [0.12597594444444443, -0.10374517305555556, -0.0, 72.0545718, 5.630568] axes_order: zyxyz - inputs: [v2, v3] name: v23tosky - outputs: [ra, dec] - ! frame: ! axes_names: [lon, lat] @@ -328,20 +245,20 @@ X,,R ؐ=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 --- -- 12376 -- 12446 -- 12516 -- 12586 -- 12656 -- 12910 -- 13164 -- 13418 -- 13672 -- 13758 -- 13844 -- 13930 -- 14016 -- 14270 -- 14524 -- 14594 +- 9552 +- 9622 +- 9692 +- 9762 +- 9832 +- 10086 +- 10340 +- 10594 +- 10848 +- 10934 +- 11020 +- 11106 +- 11192 +- 11446 +- 11700 +- 11770 ... diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 85aca436..2b34443d 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -496,7 +496,7 @@ def test_to_fits_sip(): 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_error=1.e-10, max_inv_error=0.1) + mirisip = miriwcs.to_fits_sip(bounding_box, max_pix_error=1.e-1, max_inv_pix_error=0.1, verbose=True) fitssip = astwcs.WCS(mirisip) fitsvals = fitssip.all_pix2world(warg + 1, 1) assert_allclose(miriwcs(xflat, yflat), diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 35262506..36496bb0 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -642,24 +642,34 @@ def fix_inputs(self, fixed): new_pipeline.extend(self.pipeline[1:]) return self.__class__(new_pipeline) - def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, - max_inv_error=None, max_inv_rms=None, verbose=False): + def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, + max_inv_pix_error=0.25, inv_degree=None, 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 : a pair of tuples, each consisting of two integers Represents the range of pixel values in both dimensions ((xmin, xmax), (ymin, ymax)) - max_error : float - Maximum allowed error over the domain of the pixel array. + max_pix_error : float + 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 - Degree of the SIP polynomial. If supplied, max_error must be None. + Degree of the SIP polynomial. If supplied, max_pixel_error is ignored. max_inv_error : float - Maximum allowed inverse error over the domain of the pixel array. + Maximum allowed inverse error over the domain of the pixel array + in pixel units. If None, no inverse is generated. + inv_degree : int + Degree of the inverse SIP polynomial. If supplied max_inv_pixel_error + is ignored. verbose : bool print progress of fits @@ -683,18 +693,11 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, """ - if (self.forward_transform.n_inputs !=2 - or self.forward_transform.n_outputs != 2): - raise ValueError("The to_fits_sip method only works with 2-dimensional transforms") - - # Handle the options related to degree and thresholds. - if degree is not None and max_error is not None: - raise ValueError( - "Degree cannot be specified if max_error has a value") - if degree is None and max_error is None and max_rms is None: + if (self.forward_transform.n_inputs != 2 + or self.forward_transform.n_outputs != 2): raise ValueError( - "max_error must be specified if degree is None") - fits_dict = {} + "The to_fits_sip method only works with 2-dimensional transforms") + transform = self.forward_transform # Determine reference points. (xmin, xmax), (ymin, ymax) = bounding_box @@ -724,19 +727,31 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, u = x - crpix1 v = y - crpix2 undist_x, undist_y = ntransform(x, y) + # 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) + dist1 = np.sqrt((xx - x0)**2 + (xy-y0)**2) + dist2 = np.sqrt((yx - x0)**2 + (yy-y0)**2) + # Average the distances in the two different directions. + plate_scale = 2. / (dist1 + dist2) + max_error = max_pix_error / plate_scale # The fitting section. # Fit the forward case. - fit_poly_x, fit_poly_y, max_resid = _fit_2D_poly(degree, max_error, - u, v, undist_x, undist_y, + fit_poly_x, fit_poly_y, max_resid = _fit_2D_poly(ntransform, degree, max_error, + bounding_box, verbose=verbose) 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 - fit_inv_poly_u, fit_inv_poly_v, max_inv_resid = _fit_2D_poly(degree, - max_inv_error, - U, V, u - U, v - V, + fit_inv_poly_u, fit_inv_poly_v, max_inv_resid = _fit_2D_poly(ntransform, None, + max_inv_pix_error, + bounding_box, + inverse_fit=True, + #U, V, u - U, v - V, verbose=verbose) pdegree = fit_poly_x.degree if pdegree > 1: @@ -750,8 +765,8 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, _store_2D_coefficients(hdr, sip_poly_y, 'B') _store_2D_coefficients(hdr, fit_inv_poly_u, 'AP', keeplinear=True) _store_2D_coefficients(hdr, fit_inv_poly_v, 'BP', keeplinear=True) - hdr['sipmxerr'] = (max_resid, 'Maximum difference from the GWCS model SIP is fit to.') - hdr['sipiverr'] = (max_inv_resid, 'Maximum difference for the inverse transform') + hdr['sipmxerr'] = (max_resid, 'Max diff from GWCS model.') + hdr['sipiverr'] = (max_inv_resid, 'Max diff for inverse transform') else: hdr['ctype1'] = 'RA---TAN' @@ -764,35 +779,87 @@ def to_fits_sip(self, bounding_box, max_error=None, max_rms=None, degree=None, return hdr -def _fit_2D_poly(degree, max_error, u, v, x, y, verbose=False): +def _fit_2D_poly(ntransform, degree, max_error, bounding_box, + inverse_fit=False, verbose=False): """ - Fit a pair of ordinary 2D polynomial to the supplied inputs (u, v) and - outputs (x, y) + Fit a pair of ordinary 2D polynomial to the supplied transform. + + Use adaptive sampling based on the current degree. """ llsqfitter = LinearLSQFitter() + # The case of one pass with the specified polynomial degree - if degree is None: - deglist = range(1,10) - else: + if degree: deglist = [degree] + else: + deglist = range(10) + prev_max_error = -1 + # Use 32 points in each direction since numpy is not much faster for fewer + # than 1000 values + npoints = 32 + u, v = _make_sampling_grid(npoints, bounding_box) + undist_x, undist_y = ntransform(u, v) + 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, u, v, x) - fit_poly_y = llsqfitter(poly_y, u, v, y) - max_resid = _compute_distance_residual(x, y, - fit_poly_x(u, v), fit_poly_y(u, v)) + if not inverse_fit: + fit_poly_x = llsqfitter(poly_x, u, v, undist_x) + fit_poly_y = llsqfitter(poly_y, u, v, undist_y) + max_resid = _compute_distance_residual(undist_x, undist_y, + fit_poly_x(u, v), + fit_poly_y(u, v)) + else: + fit_poly_x = llsqfitter(poly_x, undist_x, undist_y, u) + fit_poly_y = llsqfitter(poly_y, undist_x, undist_y, v) + max_resid = _compute_distance_residual(u, v, + fit_poly_x(undist_x, undist_y), + fit_poly_y(undist_x, undist_y )) + if prev_max_error < 0: + prev_max_error = max_resid + else: + if max_resid > prev_max_error: + raise RuntimeError('Failed to achieve required error tolerance') if verbose: - print(f'Degree = {degree}, max_resid = {max_resid}') - if max_error is not None and max_resid < max_error: + print(f'Degree = {deg}, max_resid = {max_resid}') + if max_resid < max_error: + # Check to see if double sampling meets error requirement. + ud, vd = _make_sampling_grid(2 * npoints, bounding_box) + undist_xd, undist_yd = ntransform(ud, vd) + if not inverse_fit: + max_resid = _compute_distance_residual(undist_xd, undist_yd, + fit_poly_x(ud, vd), + fit_poly_y(ud, vd)) + else: + max_resid = _compute_distance_residual(ud, vd, + fit_poly_x(undist_xd, undist_yd), + fit_poly_y(undist_xd, undist_yd)) if verbose: - print('terminating condition met') - break - if degree is None and deg == 9: - raise RuntimeError("Unable to achieve required fit accuracy with SIP polynomials") + 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_x, : xmax + 1: stepsize_y] + 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. From 15ca53eff7df4327938fc785cfdaee16cea0ed70 Mon Sep 17 00:00:00 2001 From: James Davies Date: Fri, 6 Mar 2020 10:02:31 -0500 Subject: [PATCH 23/38] Remove pytest runtime dependency --- gwcs/tags/spectroscopy_models.py | 16 +++++++++++++--- 1 file changed, 13 insertions(+), 3 deletions(-) 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) From 0502b46ac28df3af83bc20c3c223cade6df046ba Mon Sep 17 00:00:00 2001 From: Perry Greenfield Date: Sun, 8 Mar 2020 13:57:20 -0400 Subject: [PATCH 24/38] fixed automatic mode --- gwcs/tests/test_wcs.py | 12 +++--- gwcs/wcs.py | 94 ++++++++++++++++++++++-------------------- 2 files changed, 55 insertions(+), 51 deletions(-) diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 2b34443d..bdf74dde 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -496,10 +496,10 @@ def test_to_fits_sip(): 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_pix_error=1.e-1, max_inv_pix_error=0.1, verbose=True) + mirisip = miriwcs.to_fits_sip(bounding_box, max_inv_pix_error=0.1, verbose=True) #, max_pix_error=1.e-6, max_inv_pix_error=0.1) fitssip = astwcs.WCS(mirisip) - fitsvals = fitssip.all_pix2world(warg + 1, 1) - assert_allclose(miriwcs(xflat, yflat), - fitsvals.transpose(), atol=1e-10, rtol=0) - fits_inverse_vals = fitssip.all_world2pix(fitsvals, 1) - assert_allclose(warg, fits_inverse_vals - 1, atol=0.1, rtol=0) + fitsvalx, fitsvaly = fitssip.all_pix2world(xflat+1, yflat+1, 1) + gwcsvalx, gwcsvaly = miriwcs(xflat, yflat) + assert_allclose(gwcsvalx, fitsvalx, atol=1e-10, rtol=0) + #fits_inverse_vals = fitssip.all_world2pix(fitsvals, 1) + #assert_allclose(warg, fits_inverse_vals - 1, atol=0.1, rtol=0) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 36496bb0..ab45b9f6 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -643,7 +643,8 @@ def fix_inputs(self, fixed): return self.__class__(new_pipeline) def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, - max_inv_pix_error=0.25, inv_degree=None, verbose=False): + 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 @@ -718,69 +719,73 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, hdr['cd1_2'] = 0 hdr['cd2_1'] = 0 hdr['cd2_2'] = 0 - # Now rotate to native system and deproject - ntransform = (transform + # 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()) - # Evaluate this transform on all pixels in bounding box region - y, x = np.mgrid[ymin:ymax, xmin:xmax] # Check on the endpoints needing to be 1 larger - u = x - crpix1 - v = y - crpix2 - undist_x, undist_y = ntransform(x, y) + u, v = _make_sampling_grid(npoints, bounding_box) + undist_x, undist_y = ntransform(u, v) + # undist_x, undist_y = ntransform(x, y) # 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) - dist1 = np.sqrt((xx - x0)**2 + (xy-y0)**2) - dist2 = np.sqrt((yx - x0)**2 + (yy-y0)**2) + dist1 = np.sqrt((xx - x0)**2 + (xy - y0)**2) + dist2 = np.sqrt((yx - x0)**2 + (yy - y0)**2) # Average the distances in the two different directions. - plate_scale = 2. / (dist1 + dist2) - max_error = max_pix_error / plate_scale + plate_scale = (dist1 + dist2) / 2. + max_error = max_pix_error * plate_scale # The fitting section. - # Fit the forward case. - fit_poly_x, fit_poly_y, max_resid = _fit_2D_poly(ntransform, degree, max_error, + fit_poly_x, fit_poly_y, max_resid = _fit_2D_poly(ntransform, npoints, + degree, max_error, bounding_box, verbose=verbose) - 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 - fit_inv_poly_u, fit_inv_poly_v, max_inv_resid = _fit_2D_poly(ntransform, None, + # 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 + 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, bounding_box, inverse_fit=True, - #U, V, u - U, v - V, + UV=(U, V), verbose=verbose) pdegree = fit_poly_x.degree if pdegree > 1: hdr['a_order'] = pdegree hdr['b_order'] = pdegree - ipdegree = fit_inv_poly_u.degree - hdr['ap_order'] = ipdegree - hdr['bp_order'] = ipdegree - cd, sip_poly_x, sip_poly_y = _reform_poly_coefficients(fit_poly_x, fit_poly_y) _store_2D_coefficients(hdr, sip_poly_x, 'A') _store_2D_coefficients(hdr, sip_poly_y, 'B') - _store_2D_coefficients(hdr, fit_inv_poly_u, 'AP', keeplinear=True) - _store_2D_coefficients(hdr, fit_inv_poly_v, 'BP', keeplinear=True) - hdr['sipmxerr'] = (max_resid, 'Max diff from GWCS model.') - hdr['sipiverr'] = (max_inv_resid, 'Max diff for inverse transform') - + if max_inv_pix_error: + _store_2D_coefficients(hdr, fit_inv_poly_u, 'AP', keeplinear=True) + _store_2D_coefficients(hdr, fit_inv_poly_v, 'BP', keeplinear=True) + hdr['sipmxerr'] = (max_resid * plate_scale, 'Max diff from GWCS (equiv pix).') + if max_inv_pix_error: + ipdegree = fit_inv_poly_u.degree + hdr['ap_order'] = ipdegree + hdr['bp_order'] = ipdegree + hdr['sipiverr'] = (max_inv_resid, 'Max diff for inverse (pixels)') else: hdr['ctype1'] = 'RA---TAN' hdr['ctype2'] = 'DEC--TAN' - hdr['cd1_1'] = cd[0][0] - hdr['cd1_2'] = cd[0][1] - hdr['cd2_1'] = cd[1][0] - hdr['cd2_2'] = cd[1][1] + 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, degree, max_error, bounding_box, - inverse_fit=False, verbose=False): +def _fit_2D_poly(ntransform, npoints, degree, max_error, bounding_box, + inverse_fit=False, UV=None, verbose=False): """ Fit a pair of ordinary 2D polynomial to the supplied transform. @@ -794,9 +799,6 @@ def _fit_2D_poly(ntransform, degree, max_error, bounding_box, else: deglist = range(10) prev_max_error = -1 - # Use 32 points in each direction since numpy is not much faster for fewer - # than 1000 values - npoints = 32 u, v = _make_sampling_grid(npoints, bounding_box) undist_x, undist_y = ntransform(u, v) if verbose: @@ -811,11 +813,13 @@ def _fit_2D_poly(ntransform, degree, max_error, bounding_box, fit_poly_x(u, v), fit_poly_y(u, v)) else: - fit_poly_x = llsqfitter(poly_x, undist_x, undist_y, u) - fit_poly_y = llsqfitter(poly_y, undist_x, undist_y, v) - max_resid = _compute_distance_residual(u, v, - fit_poly_x(undist_x, undist_y), - fit_poly_y(undist_x, undist_y )) + U, V = UV + fit_poly_x = llsqfitter(poly_x, U, V, u - U) + fit_poly_y = llsqfitter(poly_y, U, V, v - V) + max_resid = _compute_distance_residual(u-U, v-V, + fit_poly_x(U, V), + fit_poly_y(U, V)) + #print(u, v, U, ) if prev_max_error < 0: prev_max_error = max_resid else: @@ -853,7 +857,7 @@ def _make_sampling_grid(npoints, bounding_box): 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_x, : xmax + 1: stepsize_y] + y, x = np.mgrid[: ymax + 1: stepsize_y, : xmax + 1: stepsize_x] x[:, -1] = xsize - 1 y[-1, :] = ysize - 1 u = x - crpix1 From 93187bf17cdbeb965f465896aa344beed8ca11e9 Mon Sep 17 00:00:00 2001 From: perrygreenfield Date: Sun, 8 Mar 2020 14:38:27 -0400 Subject: [PATCH 25/38] refactored automatic version --- gwcs/tests/test_wcs.py | 8 +++--- gwcs/wcs.py | 55 ++++++++++++++++++------------------------ 2 files changed, 28 insertions(+), 35 deletions(-) diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index bdf74dde..8595eff1 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -496,10 +496,12 @@ def test_to_fits_sip(): 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, verbose=True) #, max_pix_error=1.e-6, max_inv_pix_error=0.1) + 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) - #fits_inverse_vals = fitssip.all_world2pix(fitsvals, 1) - #assert_allclose(warg, fits_inverse_vals - 1, atol=0.1, 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) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index ab45b9f6..8479262f 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -727,6 +727,9 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, | 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) # undist_x, undist_y = ntransform(x, y) # Determine approximate pixel scale in order to compute error threshold # from the specified pixel error. Computed at the center of the array. @@ -741,7 +744,8 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, # The fitting section. fit_poly_x, fit_poly_y, max_resid = _fit_2D_poly(ntransform, npoints, degree, max_error, - bounding_box, + 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) @@ -750,13 +754,16 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, 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, - bounding_box, - inverse_fit=True, - UV=(U, V), + U, V, u-U, v-V, + Ud, Vd, ud-Ud, vd-Vd, verbose=verbose) pdegree = fit_poly_x.degree if pdegree > 1: @@ -784,8 +791,10 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, return hdr -def _fit_2D_poly(ntransform, npoints, degree, max_error, bounding_box, - inverse_fit=False, UV=None, verbose=False): +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 polynomial to the supplied transform. @@ -799,27 +808,16 @@ def _fit_2D_poly(ntransform, npoints, degree, max_error, bounding_box, else: deglist = range(10) prev_max_error = -1 - u, v = _make_sampling_grid(npoints, bounding_box) - undist_x, undist_y = ntransform(u, v) if verbose: print(f'maximum_specified_error: {max_error}') for deg in deglist: poly_x = Polynomial2D(degree=deg) poly_y = Polynomial2D(degree=deg) - if not inverse_fit: - fit_poly_x = llsqfitter(poly_x, u, v, undist_x) - fit_poly_y = llsqfitter(poly_y, u, v, undist_y) - max_resid = _compute_distance_residual(undist_x, undist_y, - fit_poly_x(u, v), - fit_poly_y(u, v)) - else: - U, V = UV - fit_poly_x = llsqfitter(poly_x, U, V, u - U) - fit_poly_y = llsqfitter(poly_y, U, V, v - V) - max_resid = _compute_distance_residual(u-U, v-V, - fit_poly_x(U, V), - fit_poly_y(U, V)) - #print(u, v, U, ) + 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 prev_max_error < 0: prev_max_error = max_resid else: @@ -829,16 +827,9 @@ def _fit_2D_poly(ntransform, npoints, degree, max_error, bounding_box, print(f'Degree = {deg}, max_resid = {max_resid}') if max_resid < max_error: # Check to see if double sampling meets error requirement. - ud, vd = _make_sampling_grid(2 * npoints, bounding_box) - undist_xd, undist_yd = ntransform(ud, vd) - if not inverse_fit: - max_resid = _compute_distance_residual(undist_xd, undist_yd, - fit_poly_x(ud, vd), - fit_poly_y(ud, vd)) - else: - max_resid = _compute_distance_residual(ud, vd, - fit_poly_x(undist_xd, undist_yd), - fit_poly_y(undist_xd, undist_yd)) + 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: From a3c66947e81f4b4bb0b783ff2cb2992e7559ddf2 Mon Sep 17 00:00:00 2001 From: perrygreenfield Date: Mon, 16 Mar 2020 09:23:53 -0400 Subject: [PATCH 26/38] address latest review comments --- gwcs/wcs.py | 40 ++++++++++++++++++---------------------- 1 file changed, 18 insertions(+), 22 deletions(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 8479262f..8a48b67d 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -658,20 +658,23 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, bounding_box : a pair of tuples, each consisting of two integers Represents the range of pixel values in both dimensions ((xmin, xmax), (ymin, ymax)) - max_pix_error : float + 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 + degree : int, optional Degree of the SIP polynomial. If supplied, max_pixel_error is ignored. - max_inv_error : float + 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 + inv_degree : int, optional Degree of the inverse SIP polynomial. If supplied max_inv_pixel_error is ignored. - verbose : bool + 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 @@ -727,19 +730,16 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, | 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 + # Double sampling to check if sampling is sufficient. ud, vd = _make_sampling_grid(2 * npoints, bounding_box) undist_xd, undist_yd = ntransform(ud, vd) - # undist_x, undist_y = ntransform(x, y) # 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) - dist1 = np.sqrt((xx - x0)**2 + (xy - y0)**2) - dist2 = np.sqrt((yx - x0)**2 + (yy - y0)**2) - # Average the distances in the two different directions. - plate_scale = (dist1 + dist2) / 2. + 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, @@ -747,7 +747,7 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, 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. + # 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]]) @@ -771,15 +771,15 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, 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) - hdr['sipmxerr'] = (max_resid * plate_scale, 'Max diff from GWCS (equiv pix).') if max_inv_pix_error: ipdegree = fit_inv_poly_u.degree hdr['ap_order'] = ipdegree hdr['bp_order'] = ipdegree - hdr['sipiverr'] = (max_inv_resid, 'Max diff for inverse (pixels)') else: hdr['ctype1'] = 'RA---TAN' hdr['ctype2'] = 'DEC--TAN' @@ -796,9 +796,8 @@ def _fit_2D_poly(ntransform, npoints, degree, max_error, xind, yind, xoutd, youtd, verbose=False): """ - Fit a pair of ordinary 2D polynomial to the supplied transform. + Fit a pair of ordinary 2D polynomials to the supplied transform. - Use adaptive sampling based on the current degree. """ llsqfitter = LinearLSQFitter() @@ -807,7 +806,7 @@ def _fit_2D_poly(ntransform, npoints, degree, max_error, deglist = [degree] else: deglist = range(10) - prev_max_error = -1 + prev_max_error = float(np.inf) if verbose: print(f'maximum_specified_error: {max_error}') for deg in deglist: @@ -818,11 +817,8 @@ def _fit_2D_poly(ntransform, npoints, degree, max_error, max_resid = _compute_distance_residual(xout, yout, fit_poly_x(xin, yin), fit_poly_y(xin, yin)) - if prev_max_error < 0: - prev_max_error = max_resid - else: - if max_resid > prev_max_error: - raise RuntimeError('Failed to achieve required error tolerance') + 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: From 4a0cf76147a68d87a8ce3e36626d3d01536add54 Mon Sep 17 00:00:00 2001 From: perrygreenfield Date: Mon, 16 Mar 2020 15:04:11 -0400 Subject: [PATCH 27/38] add changes from review --- gwcs/wcs.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 8a48b67d..bd246b39 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -17,6 +17,7 @@ from .utils import CoordinateFrameError from .utils import _toindex from . import utils +from gwcs import coordinate_frames as cf HAS_FIX_INPUTS = True @@ -655,7 +656,7 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, Parameters ---------- - bounding_box : a pair of tuples, each consisting of two integers + bounding_box : 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 @@ -697,16 +698,15 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, """ - if (self.forward_transform.n_inputs != 2 - or self.forward_transform.n_outputs != 2): + if not isinstance(self.output_frame, cf.CelestialFrame): raise ValueError( "The to_fits_sip method only works with 2-dimensional transforms") transform = self.forward_transform # Determine reference points. (xmin, xmax), (ymin, ymax) = bounding_box - crpix1 = int((xmax - xmin) / 2) - crpix2 = int((ymax - ymin) / 2) + crpix1 = (xmax - xmin) // 2 + crpix2 = (ymax - ymin) // 2 crval1, crval2 = transform(crpix1, crpix2) hdr = fits.Header() hdr['naxis'] = 2 @@ -718,10 +718,10 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, hdr['CRPIX2'] = crpix2 + 1 hdr['CRVAL1'] = crval1 hdr['CRVAL2'] = crval2 - hdr['cd1_1'] = 0 # Placeholders + hdr['cd1_1'] = 1 # Placeholders for FITS card order, all will change. hdr['cd1_2'] = 0 hdr['cd2_1'] = 0 - hdr['cd2_2'] = 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. From 5482b358a4ccc74f7ee0834295123d84aa6ad1c1 Mon Sep 17 00:00:00 2001 From: perrygreenfield Date: Mon, 16 Mar 2020 15:11:47 -0400 Subject: [PATCH 28/38] updated exception message for not getting celestial frame transform --- gwcs/wcs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index bd246b39..1be57114 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -700,7 +700,7 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, """ if not isinstance(self.output_frame, cf.CelestialFrame): raise ValueError( - "The to_fits_sip method only works with 2-dimensional transforms") + "The to_fits_sip method only works with celestial frame transforms") transform = self.forward_transform # Determine reference points. From a385f7256aae755cadd4cb8bc9daada45cb85263 Mon Sep 17 00:00:00 2001 From: perrygreenfield Date: Mon, 16 Mar 2020 16:57:31 -0400 Subject: [PATCH 29/38] using fixed asdf file for testing to_fits_sip --- gwcs/tests/data/miriwcs.asdf | 44 ++++++++++++++++++++---------------- 1 file changed, 24 insertions(+), 20 deletions(-) diff --git a/gwcs/tests/data/miriwcs.asdf b/gwcs/tests/data/miriwcs.asdf index bde9439b..e70bc620 100644 --- a/gwcs/tests/data/miriwcs.asdf +++ b/gwcs/tests/data/miriwcs.asdf @@ -4,15 +4,18 @@ %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.1} + 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.1} + software: {name: asdf, version: 2.5.2} wcs: ! name: '' steps: @@ -223,10 +226,11 @@ wcs: ! forward: - !transform/scale-1.2.0 {factor: 0.0002777777777777778} - !transform/scale-1.2.0 {factor: 0.0002777777777777778} - - ! - angles: [0.12597594444444443, -0.10374517305555556, -0.0, 72.0545718, 5.630568] + - !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] @@ -245,20 +249,20 @@ X,,R ؐ=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 --- -- 9552 -- 9622 -- 9692 -- 9762 -- 9832 -- 10086 -- 10340 -- 10594 -- 10848 -- 10934 -- 11020 -- 11106 -- 11192 -- 11446 -- 11700 -- 11770 +- 9730 +- 9800 +- 9870 +- 9940 +- 10010 +- 10264 +- 10518 +- 10772 +- 11026 +- 11112 +- 11198 +- 11284 +- 11370 +- 11624 +- 11878 +- 11948 ... From 2265530f6810ea5285eb7fc97b0b37f504e86fda Mon Sep 17 00:00:00 2001 From: perrygreenfield Date: Tue, 17 Mar 2020 10:28:01 -0400 Subject: [PATCH 30/38] fix flake issues --- gwcs/tests/test_wcs.py | 1 - gwcs/wcs.py | 4 +--- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index 8595eff1..a8ef665c 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -492,7 +492,6 @@ 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]) - warg = np.stack((xflat, yflat)).transpose() af = asdf.open(get_pkg_data_filename('data/miriwcs.asdf')) miriwcs = af.tree['wcs'] bounding_box = ((0, 1024), (0, 1024)) diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 1be57114..1cc1e6d3 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -6,9 +6,7 @@ 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, Mapping, - AffineTransformation2D, Pix2Sky_TAN, - RotateNative2Celestial, Identity) + RotateCelestial2Native) from astropy.modeling.fitting import LinearLSQFitter import astropy.io.fits as fits From 5ec523d20a35502afd53be0f41090bdf625e5125 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Mon, 23 Mar 2020 12:07:14 -0400 Subject: [PATCH 31/38] add ucd-to-ctype mapping --- gwcs/coordinate_frames.py | 34 +++++++++++++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/gwcs/coordinate_frames.py b/gwcs/coordinate_frames.py index 20dc95fb..5bb3bffc 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,11 @@ "GALACTIC_CENTER", "LOCAL_GROUP_CENTER"] +def get_ctype_from_ucd(ucd): + """ Return CTYPE string corresponding to a UCD1 value.""" + return UCD1_TO_CTYPE.get(ucd, "") + + class CoordinateFrame: """ Base class for Coordinate Frames. @@ -475,7 +508,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 From eaa795fba4bdf3c2f0570f1b51728a6cbdafb720 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Mon, 23 Mar 2020 14:27:25 -0400 Subject: [PATCH 32/38] add change log entry --- CHANGES.rst | 1 + gwcs/coordinate_frames.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index b2180948..82582816 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -8,6 +8,7 @@ New Features - Added ``to_fits_sip`` method to generate FITS header with SIP keywords [#286] +- Added ``get_ctype_from_ucd`` function. [#288] 0.12.0 (2019-12-24) ------------------- diff --git a/gwcs/coordinate_frames.py b/gwcs/coordinate_frames.py index 5bb3bffc..bcd235ad 100644 --- a/gwcs/coordinate_frames.py +++ b/gwcs/coordinate_frames.py @@ -54,7 +54,7 @@ def get_ctype_from_ucd(ucd): - """ Return CTYPE string corresponding to a UCD1 value.""" + """ Return the FITS ``CTYPE`` corresponding to a UCD1 value.""" return UCD1_TO_CTYPE.get(ucd, "") From 02d3d67a53c7866dac497b88b57d63c46f3e058b Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Mon, 23 Mar 2020 18:15:54 -0400 Subject: [PATCH 33/38] make bounding-box optional --- gwcs/coordinate_frames.py | 14 +++++++++++++- gwcs/wcs.py | 21 ++++++++++++++------- 2 files changed, 27 insertions(+), 8 deletions(-) diff --git a/gwcs/coordinate_frames.py b/gwcs/coordinate_frames.py index bcd235ad..6a955024 100644 --- a/gwcs/coordinate_frames.py +++ b/gwcs/coordinate_frames.py @@ -54,7 +54,19 @@ def get_ctype_from_ucd(ucd): - """ Return the FITS ``CTYPE`` corresponding to a UCD1 value.""" + """ + 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, "") diff --git a/gwcs/wcs.py b/gwcs/wcs.py index 1cc1e6d3..3aee00fd 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -5,7 +5,7 @@ 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, +from astropy.modeling.models import (Shift, Polynomial2D, Sky2Pix_TAN, RotateCelestial2Native) from astropy.modeling.fitting import LinearLSQFitter import astropy.io.fits as fits @@ -641,8 +641,8 @@ def fix_inputs(self, fixed): new_pipeline.extend(self.pipeline[1:]) return self.__class__(new_pipeline) - def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, - max_inv_pix_error=0.25, inv_degree=None, + 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 @@ -654,7 +654,8 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, Parameters ---------- - bounding_box : a pair of tuples, each consisting of two numbers + 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 @@ -702,7 +703,13 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, transform = self.forward_transform # Determine reference points. - (xmin, xmax), (ymin, ymax) = bounding_box + if bounding_box is not None: + (xmin, xmax), (ymin, ymax) = bounding_box + else: + try: + (xmin, xmax), (ymin, ymax) = self.bounding_box + except KeyError: + raise TypeError("A bounding_box is needed to proceed.") crpix1 = (xmax - xmin) // 2 crpix2 = (ymax - ymin) // 2 crval1, crval2 = transform(crpix1, crpix2) @@ -755,7 +762,7 @@ def to_fits_sip(self, bounding_box, max_pix_error=0.25, degree=None, 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, @@ -848,7 +855,7 @@ def _make_sampling_grid(npoints, bounding_box): 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. From 0773ba54ead191a3cc654e8031ab7659afd468da Mon Sep 17 00:00:00 2001 From: csimpson Date: Mon, 23 Mar 2020 13:52:59 -1000 Subject: [PATCH 34/38] Use CRPIXi-1 to shift pixel coordinates to reference point --- gwcs/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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) From 2ced523782105ee89b9b9951df53ea9baaaa8769 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Tue, 24 Mar 2020 09:24:12 -0400 Subject: [PATCH 35/38] fix and add test --- gwcs/tests/test_wcs.py | 12 +++++++++++- gwcs/wcs.py | 13 ++++++------- 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/gwcs/tests/test_wcs.py b/gwcs/tests/test_wcs.py index a8ef665c..a8702ab8 100644 --- a/gwcs/tests/test_wcs.py +++ b/gwcs/tests/test_wcs.py @@ -500,7 +500,17 @@ def test_to_fits_sip(): 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) + 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/wcs.py b/gwcs/wcs.py index 3aee00fd..5bd66b88 100644 --- a/gwcs/wcs.py +++ b/gwcs/wcs.py @@ -703,13 +703,12 @@ def to_fits_sip(self, bounding_box=None, max_pix_error=0.25, degree=None, transform = self.forward_transform # Determine reference points. - if bounding_box is not None: - (xmin, xmax), (ymin, ymax) = bounding_box - else: - try: - (xmin, xmax), (ymin, ymax) = self.bounding_box - except KeyError: - raise TypeError("A bounding_box is needed to proceed.") + 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) From c4009c5a6b384f93b718ad6fcccc17553ae0bd30 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Thu, 26 Mar 2020 14:59:35 -0400 Subject: [PATCH 36/38] fix test --- gwcs/tests/test_utils.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/gwcs/tests/test_utils.py b/gwcs/tests/test_utils.py index b500b954..729bded0 100644 --- a/gwcs/tests/test_utils.py +++ b/gwcs/tests/test_utils.py @@ -30,7 +30,7 @@ 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(): @@ -61,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) From bd39ab26ff2ef0c94b6fc472e45aaed8659d0c84 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Thu, 26 Mar 2020 15:01:15 -0400 Subject: [PATCH 37/38] add chanage log --- CHANGES.rst | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index 82582816..483ddd4a 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -10,6 +10,11 @@ New Features - 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 From 2d089716f70a23e5d0cb9392ec5493cb7f1cf9e2 Mon Sep 17 00:00:00 2001 From: Nadia Dencheva Date: Thu, 26 Mar 2020 15:28:33 -0400 Subject: [PATCH 38/38] prepare release --- CHANGES.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index 483ddd4a..6970afc2 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,4 +1,4 @@ -0.12.1 (Unreleased) +0.13.0 (2020-03-26) ------------------- New Features ^^^^^^^^^^^^