Source code for gala.coordinates.greatcircle

# Built-in
from textwrap import dedent
from warnings import warn

# Third-party
import astropy.units as u
import astropy.coordinates as coord
from astropy.coordinates.transformations import (DynamicMatrixTransform,
                                                 FunctionTransform)
from astropy.coordinates.attributes import (CoordinateAttribute,
                                            QuantityAttribute)
from astropy.coordinates.matrix_utilities import (rotation_matrix,
                                                  matrix_product,
                                                  matrix_transpose)
from astropy.coordinates.baseframe import base_doc
from astropy.utils.decorators import format_doc
import numpy as np

__all__ = ['GreatCircleICRSFrame', 'make_greatcircle_cls',
           'pole_from_endpoints']


def greatcircle_to_greatcircle(from_greatcircle_coord,
                               to_greatcircle_frame):
    """Transform between two greatcircle frames."""

    # This transform goes through the parent frames on each side.
    # from_frame -> from_frame.origin -> to_frame.origin -> to_frame
    intermediate_from = from_greatcircle_coord.transform_to(
        from_greatcircle_coord.pole)
    intermediate_to = intermediate_from.transform_to(
        to_greatcircle_frame.pole)
    return intermediate_to.transform_to(to_greatcircle_frame)


def reference_to_greatcircle(reference_frame, greatcircle_frame):
    """Convert a reference coordinate to a great circle frame."""

    # Define rotation matrices along the position angle vector, and
    # relative to the origin.
    pole = greatcircle_frame.pole.transform_to(coord.ICRS)
    ra0 = greatcircle_frame.ra0
    center = greatcircle_frame.center
    R_rot = rotation_matrix(greatcircle_frame.rotation, 'z')

    if not np.isnan(ra0) and np.abs(pole.dec.value) > 1e-15:
        zaxis = pole.cartesian.xyz.value
        xaxis = np.array([np.cos(ra0), np.sin(ra0), 0.])
        if np.abs(zaxis[2]) >= 1e-15:
            xaxis[2] = -(zaxis[0]*xaxis[0] + zaxis[1]*xaxis[1]) / zaxis[2]
        else:
            xaxis[2] = 0.
        xaxis = xaxis / np.sqrt(np.sum(xaxis**2))

        yaxis = np.cross(zaxis, xaxis)
        yaxis = yaxis / np.sqrt(np.sum(yaxis**2))

        R = np.stack((xaxis, yaxis, zaxis))

    elif center is not None:
        R1 = rotation_matrix(pole.ra, 'z')
        R2 = rotation_matrix(90*u.deg - pole.dec, 'y')
        Rtmp = matrix_product(R2, R1)

        rot = center.cartesian.transform(Rtmp)
        rot_lon = rot.represent_as(coord.UnitSphericalRepresentation).lon
        R3 = rotation_matrix(rot_lon, 'z')
        R = matrix_product(R3, R2, R1)

    else:
        if not np.isnan(ra0) and np.abs(pole.dec.value) < 1e-15:
            warn("Ignoring input ra0 because the pole is along dec=0",
                 RuntimeWarning)

        R1 = rotation_matrix(pole.ra, 'z')
        R2 = rotation_matrix(90*u.deg - pole.dec, 'y')
        R = matrix_product(R2, R1)

    return matrix_product(R_rot, R)


def greatcircle_to_reference(greatcircle_coord, reference_frame):
    """Convert an great circle frame coordinate to the reference frame"""

    # use the forward transform, but just invert it
    R = reference_to_greatcircle(reference_frame, greatcircle_coord)
    # transpose is the inverse because R is a rotation matrix
    return matrix_transpose(R)


def greatcircle_transforms(self_transform=False):
    def set_greatcircle_transforms(cls):
        DynamicMatrixTransform(reference_to_greatcircle,
                               coord.ICRS, cls,
                               register_graph=coord.frame_transform_graph)

        DynamicMatrixTransform(greatcircle_to_reference,
                               cls, coord.ICRS,
                               register_graph=coord.frame_transform_graph)

        if self_transform:
            FunctionTransform(greatcircle_to_greatcircle,
                              cls, cls,
                              register_graph=coord.frame_transform_graph)
        return cls

    return set_greatcircle_transforms


_components = """
phi1 : `~astropy.units.Quantity`
    Longitude component.
phi2 : `~astropy.units.Quantity`
    Latitude component.
distance : `~astropy.units.Quantity`
    Distance.

pm_phi1_cosphi2 : `~astropy.units.Quantity`
    Proper motion in longitude.
pm_phi2 : `~astropy.units.Quantity`
    Proper motion in latitude.
radial_velocity : `~astropy.units.Quantity`
    Line-of-sight or radial velocity.
"""

_footer = """
Frame attributes
----------------
pole : `~astropy.coordinates.SkyCoord`, `~astropy.coordinates.ICRS`
    The coordinate specifying the pole of this frame.
ra0 : `~astropy.units.Quantity`, `~astropy.coordinates.Angle` (optional)
    If specified, an additional transformation will be applied to make
    this right ascension the longitude zero-point of the resulting
    coordinate frame.
rotation : `~astropy.units.Quantity`, `~astropy.coordinates.Angle` (optional)
    If specified, a final rotation about the pole (i.e. the resulting z
    axis) applied.
"""


[docs]@format_doc(dedent(base_doc), components=_components, footer=_footer) @greatcircle_transforms(self_transform=True) class GreatCircleICRSFrame(coord.BaseCoordinateFrame): """A frame rotated into great circle coordinates with the pole and longitude specified as frame attributes. ``GreatCircleICRSFrame``s always have component names for spherical coordinates of ``phi1``/``phi2``. """ pole = CoordinateAttribute(default=None, frame=coord.ICRS) center = CoordinateAttribute(default=None, frame=coord.ICRS) ra0 = QuantityAttribute(default=np.nan*u.deg, unit=u.deg) rotation = QuantityAttribute(default=0, unit=u.deg) frame_specific_representation_info = { coord.SphericalRepresentation: [ coord.RepresentationMapping('lon', 'phi1'), coord.RepresentationMapping('lat', 'phi2'), coord.RepresentationMapping('distance', 'distance')] } default_representation = coord.SphericalRepresentation default_differential = coord.SphericalCosLatDifferential _default_wrap_angle = 180*u.deg def __init__(self, *args, **kwargs): wrap = kwargs.pop('wrap_longitude', True) super().__init__(*args, **kwargs) if wrap and isinstance(self._data, (coord.UnitSphericalRepresentation, coord.SphericalRepresentation)): self._data.lon.wrap_angle = self._default_wrap_angle if self.center is not None and np.isfinite(self.ra0): raise ValueError("Both `center` and `ra0` were specified for this " "{} object: you can only specify one or the other." .format(self.__class__.__name__)) # TODO: remove this. This is a hack required as of astropy v3.1 in order # to have the longitude components wrap at the desired angle
[docs] def represent_as(self, base, s='base', in_frame_units=False): r = super().represent_as(base, s=s, in_frame_units=in_frame_units) r.lon.wrap_angle = self._default_wrap_angle return r
represent_as.__doc__ = coord.BaseCoordinateFrame.represent_as.__doc__
[docs] @classmethod def from_endpoints(cls, coord1, coord2, ra0=None, rotation=None): """Compute the great circle frame from two endpoints of an arc on the unit sphere. Parameters ---------- coord1 : `~astropy.coordinates.SkyCoord` One endpoint of the great circle arc. coord2 : `~astropy.coordinates.SkyCoord` The other endpoint of the great circle arc. ra0 : `~astropy.units.Quantity`, `~astropy.coordinates.Angle` (optional) If specified, an additional transformation will be applied to make this right ascension the longitude zero-point of the resulting coordinate frame. rotation : `~astropy.units.Quantity`, `~astropy.coordinates.Angle` (optional) If specified, a final rotation about the pole (i.e. the resulting z axis) applied. """ pole = pole_from_endpoints(coord1, coord2) kw = dict(pole=pole) if ra0 is not None: kw['ra0'] = ra0 if rotation is not None: kw['rotation'] = rotation if ra0 is None and rotation is None: midpt = sph_midpoint(coord1, coord2) kw['ra0'] = midpt.ra return cls(**kw)
[docs] @classmethod def from_xyz(cls, xnew=None, ynew=None, znew=None): """Compute the great circle frame from a specification of the coordinate axes in the new system. Parameters ---------- xnew : astropy ``Representation`` object The x-axis in the new system. ynew : astropy ``Representation`` object The y-axis in the new system. znew : astropy ``Representation`` object The z-axis in the new system. """ is_none = [xnew is None, ynew is None, znew is None] if np.sum(is_none) > 1: raise ValueError("At least 2 axes must be specified.") if xnew is not None: xnew = xnew.to_cartesian() if ynew is not None: ynew = ynew.to_cartesian() if znew is not None: znew = znew.to_cartesian() if znew is None: znew = xnew.cross(ynew) if ynew is None: ynew = -xnew.cross(znew) if xnew is None: xnew = ynew.cross(znew) pole = coord.SkyCoord(znew, frame='icrs') center = coord.SkyCoord(xnew, frame='icrs') return cls(pole=pole, center=center)
[docs]def make_greatcircle_cls(cls_name, docstring_header=None, **kwargs): @format_doc(base_doc, components=_components, footer=_footer) @greatcircle_transforms(self_transform=False) class GCFrame(GreatCircleICRSFrame): pole = kwargs.get('pole', None) ra0 = kwargs.get('ra0', np.nan*u.deg) center = kwargs.get('center', None) rotation = kwargs.get('rotation', 0*u.deg) GCFrame.__name__ = cls_name if docstring_header: GCFrame.__doc__ = "{0}\n{1}".format(docstring_header, GCFrame.__doc__) return GCFrame
[docs]def pole_from_endpoints(coord1, coord2): """Compute the pole from a great circle that connects the two specified coordinates. This assumes a right-handed rule from coord1 to coord2: the pole is the north pole under that assumption. Parameters ---------- coord1 : `~astropy.coordinates.SkyCoord` Coordinate of one point on a great circle. coord2 : `~astropy.coordinates.SkyCoord` Coordinate of the other point on a great circle. Returns ------- pole : `~astropy.coordinates.SkyCoord` The coordinates of the pole. """ c1 = coord1.cartesian / coord1.cartesian.norm() coord2 = coord2.transform_to(coord1.frame) c2 = coord2.cartesian / coord2.cartesian.norm() pole = c1.cross(c2) pole = pole / pole.norm() return coord1.frame.realize_frame(pole)
def sph_midpoint(coord1, coord2): """Compute the midpoint between two points on the sphere. Parameters ---------- coord1 : `~astropy.coordinates.SkyCoord` Coordinate of one point on a great circle. coord2 : `~astropy.coordinates.SkyCoord` Coordinate of the other point on a great circle. Returns ------- midpt : `~astropy.coordinates.SkyCoord` The coordinates of the spherical midpoint. """ c1 = coord1.cartesian / coord1.cartesian.norm() coord2 = coord2.transform_to(coord1.frame) c2 = coord2.cartesian / coord2.cartesian.norm() midpt = 0.5 * (c1 + c2) usph = midpt.represent_as(coord.UnitSphericalRepresentation) return coord1.frame.realize_frame(usph)