# 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
from astropy.coordinates.baseframe import base_doc
from astropy.utils.decorators import format_doc
import numpy as np
from .helpers import StringValidatedAttribute
__all__ = ["GreatCircleICRSFrame", "make_greatcircle_cls", "pole_from_endpoints"]
def get_xhat(zhat, ra0, tol=1e-10):
"""
Helper to get the x-hat vector along a great circle defined by the input zhat that
intersects with the specified longitude (ra0).
"""
ra0 = 90 * u.deg - ra0
z1, z2, z3 = zhat
if np.isclose(z3, 0, atol=tol): # pole in x-y - can't satisfy ra0
raise ValueError(
"Pole is in the x-y plane, so can't satisfy the ra0 requirement"
)
denom = (
z2**2
+ z3**2
+ 2 * z1 * z2 * np.tan(ra0)
+ (z1**2 + z3**2) * np.tan(ra0) ** 2
)
x1 = -np.tan(ra0) * np.sqrt(z3**2 / denom)
x2 = x1 / np.tan(ra0)
if np.isclose(z3, 1, atol=tol):
x3 = 0.0
else:
x3 = (z2 + z1 * np.tan(ra0)) * np.abs(x2) / z3
return np.array([x1, x2, x3])
def get_origin_from_pole_ra0(pole, ra0, origin_disambiguate=None):
"""
Figure out the coordinate system origin (i.e. the x-axis, expressed in the old
coordinate frame). Given just a pole and ra0, there is an ambiguity to the direction
of the x-axis because the two great circles (defined by pole and ra0) intersect at
two points. To resolve this ambiguity, you can specify ``origin_disambiguate``,
which is a coordinate in the old system (ICRS) used to pick the x-axis closest to
that location. If this is not specified, it uses (RA, Dec)=(0, 0).
"""
if origin_disambiguate is None:
origin_disambiguate = coord.SkyCoord(0, 0, unit=u.deg, frame=pole)
# figure out origin from ra0
zhat = np.squeeze((pole.cartesian / pole.cartesian.norm()).xyz)
xhat1 = coord.CartesianRepresentation(get_xhat(zhat, ra0))
xhat2 = -xhat1
origin1 = coord.SkyCoord(xhat1, frame=pole, representation_type="unitspherical")
origin2 = coord.SkyCoord(xhat2, frame=pole, representation_type="unitspherical")
sep1 = origin_disambiguate.separation(origin1).to_value(u.deg)
sep2 = origin_disambiguate.separation(origin2).to_value(u.deg)
# Convention:
if sep1 <= sep2:
return origin1
else:
return origin2
[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.
"""
cart1 = coord1.cartesian.without_differentials()
cart2 = coord2.cartesian.without_differentials()
if isinstance(coord1, coord.SkyCoord):
frame1 = coord1.frame
elif isinstance(coord1, coord.BaseCoordinateFrame):
frame1 = coord1
else:
raise TypeError(
"Input coordinate must be a SkyCoord or coordinate frame instance."
)
c1 = cart1 / cart1.norm()
coord2 = coord2.transform_to(frame1)
c2 = cart2 / cart2.norm()
pole = c1.cross(c2)
pole = pole / pole.norm()
return frame1.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.
"""
cart1 = coord1.cartesian.without_differentials()
cart2 = coord2.cartesian.without_differentials()
if isinstance(coord1, coord.SkyCoord):
frame1 = coord1.frame
elif isinstance(coord1, coord.BaseCoordinateFrame):
frame1 = coord1
else:
raise TypeError(
"Input coordinate must be a SkyCoord or coordinate frame instance."
)
c1 = cart1 / cart1.norm()
coord2 = coord2.transform_to(frame1)
c2 = cart2 / cart2.norm()
midpt = 0.5 * (c1 + c2)
usph = midpt.represent_as(coord.UnitSphericalRepresentation)
return frame1.realize_frame(usph)
def ensure_orthogonal(pole, origin, priority="origin", tol=1e-10):
"""
Makes sure the pole and origin are unit vectors, and are orthogonal. Adjusts either
the pole or origin to make orthogonal if not.
Parameters
----------
x : array_like
Must be a unit vector.
z : array_like
Must be a unit vector.
"""
origin = origin.realize_frame(
origin.represent_as("unitspherical").without_differentials()
).squeeze()
pole = pole.realize_frame(
pole.represent_as("unitspherical").without_differentials()
).squeeze()
x = np.squeeze(origin.cartesian.xyz)
z = np.squeeze(pole.cartesian.xyz)
if np.abs(np.dot(x, z)) > tol:
if priority == "origin":
msg = "Keeping the origin fixed and adjusting the pole to be orthogonal."
z = z - (z @ x) * x
pole = pole.realize_frame(
coord.CartesianRepresentation(z), representation_type="unitspherical"
)
else: # validated by class attribute, so assume "pole"
msg = "Keeping the pole fixed and adjusting the origin to be orthogonal."
x = x - (x @ z) * z
origin = origin.realize_frame(
coord.CartesianRepresentation(x), representation_type="unitspherical"
)
warn(
f"Input origin and pole are not orthogonal. {msg} Use "
"warnings.simplefilter('ignore') to ignore this warning.",
RuntimeWarning,
)
return pole, origin
def pole_origin_to_R(pole, origin):
"""
Compute the Cartesian rotation matrix from the given pole and origin.
This functiona assumes that ``pole`` and ``origin`` are orthogonal.
"""
if not pole.is_equivalent_frame(origin):
raise ValueError("The coordinate frame of the input pole and origin must match")
xaxis = np.squeeze((origin.cartesian / origin.cartesian.norm()).xyz)
zaxis = np.squeeze((pole.cartesian / pole.cartesian.norm()).xyz)
yaxis = np.cross(zaxis, xaxis)
R = np.stack((xaxis, yaxis, zaxis))
return R
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."""
return greatcircle_frame._R
def greatcircle_to_reference(greatcircle_coord, reference_frame):
"""Convert a 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 R.T
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 pole of the new coordinate frame, defined in the old frame (ICRS).
origin : `~astropy.coordinates.SkyCoord`, `~astropy.coordinates.ICRS`
The x-axis (spherical origin) of the new coordinate frame, defined in the old
frame (ICRS).
"""
[docs]
@format_doc(dedent(base_doc), components=_components, footer=_footer)
@greatcircle_transforms(self_transform=True)
class GreatCircleICRSFrame(coord.BaseCoordinateFrame):
"""
A coordinate frame defined by a pole and origin.
``GreatCircleICRSFrame``s always have component names for spherical coordinates of
``phi1`` and ``phi2`` (so, proper motion components are ``pm_phi1_cosphi2``, etc.).
"""
pole = CoordinateAttribute(default=None, frame=coord.ICRS)
origin = CoordinateAttribute(default=None, frame=coord.ICRS)
priority = StringValidatedAttribute(
default="origin", valid_values=["origin", "pole"]
)
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):
if "ra0" in kwargs:
raise ValueError(
"Initializing a GreatCircleICRSFrame with a pole and ra0 is no longer "
"supported because this does not uniquely determine a coordinate frame."
" To initialize a frame with a pole and ra0 and ignore the ambiguity, "
"use the .from_pole_ra0() classmethod."
)
if "rotation" in kwargs:
raise ValueError(
"Initializing a GreatCircleICRSFrame with a `rotation` is no longer "
"supported."
)
wrap = kwargs.pop("wrap_longitude", True)
super().__init__(*args, **kwargs)
if self.pole is None or self.origin is None:
raise ValueError("You must specify both a pole and an origin")
pole, origin = ensure_orthogonal(self.pole, self.origin, priority=self.priority)
self._pole = pole
self._origin = origin
self._R = pole_origin_to_R(self.pole, self.origin)
if wrap and isinstance(
self._data,
(coord.UnitSphericalRepresentation, coord.SphericalRepresentation),
):
self._data.lon.wrap_angle = self._default_wrap_angle
# 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)
if hasattr(r, "lon"):
r.lon.wrap_angle = self._default_wrap_angle
return r
represent_as.__doc__ = coord.BaseCoordinateFrame.represent_as.__doc__
[docs]
@classmethod
def from_pole_ra0(cls, pole, ra0, origin_disambiguate=None):
f"""
Compute the great circle frame from a pole and RA of longitude=0.
{get_origin_from_pole_ra0.__doc__}
Parameters
----------
pole : `~astropy.coordinates.SkyCoord`
The pole of the new coordinate frame, defined in the old frame (ICRS).
ra0 : `~astropy.units.Quantity`, `~astropy.coordinates.Angle` (optional)
Right Ascension of longitude zero.
origin_disambiguate : `~astropy.coordinates.SkyCoord` (optional)
A sky coordinate in the old frame (ICRS) used to disambiguate the coordinate
system origin. The x-axis closest to this coordinate is chosen as the new
system origin / x-axis.
"""
origin = get_origin_from_pole_ra0(
pole, ra0, origin_disambiguate=origin_disambiguate
)
return cls(pole=pole, origin=origin)
[docs]
@classmethod
def from_endpoints(cls, coord1, coord2, origin=None, ra0=None, priority=None):
"""
Compute the great circle frame from two endpoints of an arc on the unit sphere.
If you specify an ``origin``, it should be orthogonal to the pole of the great
circle defined by ``coord1`` and ``coord2``. If it is not orthogonal to the
pole, by default, the pole will be adjusted along the great circle connecting
the pole to the input ``origin``. If you would instead like to keep the pole
fixed and orthogonalize the ``origin``, pass in ``priority='pole'``.
Parameters
----------
coord1 : `~astropy.coordinates.SkyCoord`
One endpoint of the great circle arc.
coord2 : `~astropy.coordinates.SkyCoord`
The other endpoint of the great circle arc.
origin : `~astropy.coordinates.SkyCoord` (optional)
The x-axis (spherical origin) of the new coordinate frame, defined in the
old frame (ICRS). This defines the (phi1,phi2)=(0,0)º coordinate.
ra0 : `~astropy.units.Quantity`, `~astropy.coordinates.Angle` (optional)
Right Ascension of longitude zero. You can only specify one of ``origin`` or
``ra0``.
priority : str (optional)
Defines the priority of keeping either the pole or origin fixed when they
are not orthogonal based on the input.
"""
if ra0 is not None and origin is not None:
raise ValueError("You can only pass one of `ra0` or `origin`, not both")
pole = pole_from_endpoints(coord1.squeeze(), coord2.squeeze())
if ra0 is not None:
midpt = sph_midpoint(coord1.squeeze(), coord2.squeeze())
origin = get_origin_from_pole_ra0(pole, ra0, midpt)
elif ra0 is None and origin is None:
origin = sph_midpoint(coord1.squeeze(), coord2.squeeze())
return cls(pole=pole, origin=origin, priority=priority)
[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")
origin = coord.SkyCoord(xnew, frame="icrs")
return cls(pole=pole, origin=origin)
[docs]
@classmethod
def from_R(cls, R):
"""
Compute the great circle frame from a rotation matrix that specifies the
transformation from ICRS to the new frame.
Parameters
----------
R : array_like
The transformation matrix.
"""
pole = coord.SkyCoord(
coord.CartesianRepresentation([0.0, 0.0, 1.0]).transform(R.T),
frame="icrs",
representation_type="unitspherical",
)
origin = coord.SkyCoord(
coord.CartesianRepresentation([1.0, 0.0, 0.0]).transform(R.T),
frame="icrs",
representation_type="unitspherical",
)
return cls(pole=pole, origin=origin)
[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)
origin = kwargs.get("origin", None)
GCFrame.__name__ = cls_name
if docstring_header:
GCFrame.__doc__ = "{0}\n{1}".format(docstring_header, GCFrame.__doc__)
return GCFrame
