Merge pull request #291 from spacetelescope/0.13.0x

0.13.0x

CHANGES.rst

@@ -1,3 +1,20 @@
+0.13.0 (2020-03-26)
+-------------------
+New Features
+^^^^^^^^^^^^
+
+- Added two new transforms - ``SphericalToCartesian`` and
+  ``CartesianToSpherical``. [#275, #284, #285]
+
+- Added ``to_fits_sip`` method to generate FITS header with SIP keywords [#286]
+
+- Added ``get_ctype_from_ucd`` function. [#288]
+
+Bug Fixes
+^^^^^^^^^
+
+- Fixed an off by one issue in ``utils.make_fitswcs_transform``. [#290]
+
 0.12.0 (2019-12-24)
 -------------------
 New Features

gwcs/coordinate_frames.py

@@ -18,13 +18,58 @@
            '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",
                                "TOPOCENTER", "LSR", "LSRK", "LSRD",
                                "GALACTIC_CENTER", "LOCAL_GROUP_CENTER"]
 
 
+def get_ctype_from_ucd(ucd):
+    """
+    Return the FITS ``CTYPE`` corresponding to a UCD1 value.
+
+    Parameters
+    ----------
+    ucd : str
+        UCD string, for example one of ```WCS.world_axis_physical_types``.
+
+    Returns
+    -------
+    CTYPE : str
+        The corresponding FITS ``CTYPE`` value or an empty string.
+    """
+    return UCD1_TO_CTYPE.get(ucd, "")
+
+
 class CoordinateFrame:
     """
     Base class for Coordinate Frames.
@@ -475,7 +520,6 @@ def _convert_to_time(self, dt, *, unit, **kwargs):
         else:
             return time.Time(dt, **kwargs)
 
-
     def coordinate_to_quantity(self, *coords):
         if isinstance(coords[0], time.Time):
             ref_value = self.reference_frame.value
@@ -629,11 +673,11 @@ def from_index(cls, indexes):
         indexes = np.asanyarray(indexes, dtype=int)
         out = np.empty_like(indexes, dtype=object)
 
-        out[nans] = np.nan
-
         for profile, index in cls.profiles.items():
             out[indexes == index] = profile
 
+        out[nans] = np.nan
+
         if out.size == 1 and not nans:
             return StokesProfile(out.item())
         elif nans.all():

gwcs/geometry.py

@@ -3,11 +3,13 @@
 Spectroscopy related models.
 """
 
+import numbers
 import numpy as np
 from astropy.modeling.core import Model
+from astropy import units as u
 
-
-__all__ = ['ToDirectionCosines', 'FromDirectionCosines']
+__all__ = ['ToDirectionCosines', 'FromDirectionCosines',
+           'SphericalToCartesian', 'CartesianToSpherical']
 
 
 class ToDirectionCosines(Model):
@@ -55,3 +57,152 @@ def evaluate(self, cosa, cosb, cosc, length):
 
     def inverse(self):
         return ToDirectionCosines()
+
+
+class SphericalToCartesian(Model):
+    """
+    Convert spherical coordinates on a unit sphere to cartesian coordinates.
+    Spherical coordinates when not provided as ``Quantity`` are assumed
+    to be in degrees with ``lon`` being the *longitude (or azimuthal) angle*
+    ``[0, 360)`` (or ``[-180, 180)``) and angle ``lat`` is the *latitude*
+    (or *elevation angle*) in range ``[-90, 90]``.
+
+    """
+    _separable = False
+
+    _input_units_allow_dimensionless = True
+
+    n_inputs = 2
+    n_outputs = 3
+
+    def __init__(self, wrap_lon_at=360, **kwargs):
+        """
+        Parameters
+        ----------
+        wrap_lon_at : {360, 180}, optional
+            An **integer number** that specifies the range of the longitude
+            (azimuthal) angle. When ``wrap_lon_at`` is 180, the longitude angle
+            will have a range of ``[-180, 180)`` and when ``wrap_lon_at``
+            is 360 (default), the longitude angle will have a range of
+            ``[0, 360)``.
+
+        """
+        super().__init__(**kwargs)
+        self.inputs = ('lon', 'lat')
+        self.outputs = ('x', 'y', 'z')
+        self.wrap_lon_at = wrap_lon_at
+
+    @property
+    def wrap_lon_at(self):
+        """ An **integer number** that specifies the range of the longitude
+        (azimuthal) angle.
+
+        Allowed values are 180 and 360. When ``wrap_lon_at``
+        is 180, the longitude angle will have a range of ``[-180, 180)`` and
+        when ``wrap_lon_at`` is 360 (default), the longitude angle will have a
+        range of ``[0, 360)``.
+
+        """
+        return self._wrap_lon_at
+
+    @wrap_lon_at.setter
+    def wrap_lon_at(self, wrap_angle):
+        if not (isinstance(wrap_angle, numbers.Integral) and wrap_angle in [180, 360]):
+            raise ValueError("'wrap_lon_at' must be an integer number: 180 or 360")
+        self._wrap_lon_at = wrap_angle
+
+    def evaluate(self, lon, lat):
+        if isinstance(lon, u.Quantity) != isinstance(lat, u.Quantity):
+            raise TypeError("All arguments must be of the same type "
+                            "(i.e., quantity or not).")
+
+        lon = np.deg2rad(lon)
+        lat = np.deg2rad(lat)
+
+        cs = np.cos(lat)
+        x = cs * np.cos(lon)
+        y = cs * np.sin(lon)
+        z = np.sin(lat)
+
+        return x, y, z
+
+    def inverse(self):
+        return CartesianToSpherical(wrap_lon_at=self._wrap_lon_at)
+
+    @property
+    def input_units(self):
+        return {'lon': u.deg, 'lat': u.deg}
+
+
+class CartesianToSpherical(Model):
+    """
+    Convert cartesian coordinates to spherical coordinates on a unit sphere.
+    Output spherical coordinates are in degrees. When input cartesian
+    coordinates are quantities (``Quantity`` objects), output angles
+    will also be quantities in degrees. Angle ``lon`` is the *longitude*
+    (or *azimuthal angle*) in range ``[0, 360)`` (or ``[-180, 180)``) and
+    angle ``lat`` is the *latitude* (or *elevation angle*) in the
+    range ``[-90, 90]``.
+
+    """
+    _separable = False
+
+    _input_units_allow_dimensionless = True
+
+    n_inputs = 3
+    n_outputs = 2
+
+    def __init__(self, wrap_lon_at=360, **kwargs):
+        """
+        Parameters
+        ----------
+        wrap_lon_at : {360, 180}, optional
+            An **integer number** that specifies the range of the longitude
+            (azimuthal) angle. When ``wrap_lon_at`` is 180, the longitude angle
+            will have a range of ``[-180, 180)`` and when ``wrap_lon_at``
+            is 360 (default), the longitude angle will have a range of
+            ``[0, 360)``.
+
+        """
+        super().__init__(**kwargs)
+        self.inputs = ('x', 'y', 'z')
+        self.outputs = ('lon', 'lat')
+        self.wrap_lon_at = wrap_lon_at
+
+    @property
+    def wrap_lon_at(self):
+        """ An **integer number** that specifies the range of the longitude
+        (azimuthal) angle.
+
+        Allowed values are 180 and 360. When ``wrap_lon_at``
+        is 180, the longitude angle will have a range of ``[-180, 180)`` and
+        when ``wrap_lon_at`` is 360 (default), the longitude angle will have a
+        range of ``[0, 360)``.
+
+        """
+        return self._wrap_lon_at
+
+    @wrap_lon_at.setter
+    def wrap_lon_at(self, wrap_angle):
+        if not (isinstance(wrap_angle, numbers.Integral) and wrap_angle in [180, 360]):
+            raise ValueError("'wrap_lon_at' must be an integer number: 180 or 360")
+        self._wrap_lon_at = wrap_angle
+
+    def evaluate(self, x, y, z):
+        nquant = [isinstance(i, u.Quantity) for i in (x, y, z)].count(True)
+        if nquant in [1, 2]:
+            raise TypeError("All arguments must be of the same type "
+                            "(i.e., quantity or not).")
+
+        h = np.hypot(x, y)
+        lat = np.rad2deg(np.arctan2(z, h))
+        lon = np.rad2deg(np.arctan2(y, x))
+        lon[h == 0] *= 0
+
+        if self._wrap_lon_at != 180:
+            lon = np.mod(lon, 360.0 * u.deg if nquant else 360.0)
+
+        return lon, lat
+
+    def inverse(self):
+        return SphericalToCartesian(wrap_lon_at=self._wrap_lon_at)

gwcs/schemas/stsci.edu/gwcs/spherical_cartesian-1.0.0.yaml

@@ -0,0 +1,44 @@
+%YAML 1.1
+---
+$schema: "http://stsci.edu/schemas/yaml-schema/draft-01"
+id: "http://stsci.edu/schemas/gwcs/spherical_cartesian-1.0.0"
+tag: "tag:stsci.edu:gwcs/spherical_cartesian-1.0.0"
+
+title: >
+  Convert coordinates between spherical and Cartesian coordinates.
+
+description: |
+  This schema is for transforms which convert between spherical coordinates
+  (on the unit sphere) and Cartesian coordinates.
+
+examples:
+  -
+    - Convert spherical coordinates to Cartesian coordinates.
+
+    - |
+
+        !<tag:stsci.edu:gwcs/spherical_cartesian-1.0.0>
+          transform_type: spherical_to_cartesian
+
+  -
+    - Convert Cartesian coordinates to spherical coordinates.
+
+    - |
+        !<tag:stsci.edu:gwcs/spherical_cartesian-1.0.0>
+          transform_type: cartesian_to_spherical
+
+allOf:
+  - $ref:  "tag:stsci.edu:asdf/transform/transform-1.1.0"
+  - object:
+    properties:
+      wrap_lon_at:
+        description: Angle at which to wrap the longitude angle.
+        type: integer
+        enum: [180, 360]
+        default: 360
+      transform_type:
+        description: The type of transform/class to initialize.
+        type: string
+        enum: [spherical_to_cartesian, cartesian_to_spherical]
+
+    required: [transform_type, wrap_lon_at]

gwcs/tags/geometry_models.py

@@ -3,10 +3,11 @@
 """
 from asdf import yamlutil
 from ..gwcs_types import GWCSTransformType
-from .. geometry import ToDirectionCosines, FromDirectionCosines
+from .. geometry import (ToDirectionCosines, FromDirectionCosines,
+                         SphericalToCartesian, CartesianToSpherical)
 
 
-__all__ = ['DirectionCosinesType']
+__all__ = ['DirectionCosinesType', 'SphericalCartesianType']
 
 
 class DirectionCosinesType(GWCSTransformType):
@@ -34,3 +35,35 @@ def to_tree_transform(cls, model, ctx):
             raise TypeError(f"Model of type {model.__class__} is not supported.")
         node = {'transform_type': transform_type}
         return yamlutil.custom_tree_to_tagged_tree(node, ctx)
+
+
+class SphericalCartesianType(GWCSTransformType):
+    name = "spherical_cartesian"
+    types = [SphericalToCartesian, CartesianToSpherical]
+    version = "1.0.0"
+
+    @classmethod
+    def from_tree_transform(cls, node, ctx):
+        transform_type = node['transform_type']
+        wrap_lon_at = node['wrap_lon_at']
+        if transform_type == 'spherical_to_cartesian':
+            return SphericalToCartesian(wrap_lon_at=wrap_lon_at)
+        elif transform_type == 'cartesian_to_spherical':
+            return CartesianToSpherical(wrap_lon_at=wrap_lon_at)
+        else:
+            raise TypeError(f"Unknown model_type {transform_type}")
+
+    @classmethod
+    def to_tree_transform(cls, model, ctx):
+        if isinstance(model, SphericalToCartesian):
+            transform_type = 'spherical_to_cartesian'
+        elif isinstance(model, CartesianToSpherical):
+            transform_type = 'cartesian_to_spherical'
+        else:
+            raise TypeError(f"Model of type {model.__class__} is not supported.")
+
+        node = {
+            'transform_type': transform_type,
+            'wrap_lon_at': model.wrap_lon_at
+        }
+        return yamlutil.custom_tree_to_tagged_tree(node, ctx)

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)