WIP: Proof of concept for new frames

Objectives:

* Abolish frames with data
  Frames only contain the parameters for the transformations
* No metaprogramming or global variables to define new frames

Now instead of a graph, there's a central frame ICRS
(the only one without obstime as parameter)
and all the other frames define to/from transformation methods.

Missing:

* Unit tests
* Regression tests
* Short-circuit transforms (for example GCRS<->ITRS or GCRS<->HCRS)
* Factories for custom frames or specific conventions
* Transformation graph (necessary?)

astropy/coordinates/new_frames/__init__.py

@@ -0,0 +1,29 @@
+"""
+Alternative implementation of coordinate frames.
+
+Objectives:
+
+* Abolish frames with data
+  Frames only contain the parameters for the transformations
+* No metaprogramming or global variables to define new frames
+
+"""
+# TODO: Short-circuit transforms (for example GCRS<->ITRS or GCRS<->HCRS)
+# TODO: Evaluate need for transformation graph
+# TODO: Factories for custom frames or specific conventions
+
+from .earth import CIRS, GCRS
+from .enums import NutationModels
+from .utils import transform
+from .ecliptic import BarycentricEcliptic, HeliocentricEcliptic
+from .equatorial import HCRS, ICRS
+
+__all__ = [
+    "CIRS",
+    "GCRS",
+    "BarycentricEcliptic",
+    "HeliocentricEcliptic",
+    "HCRS",
+    "ICRS",
+    "NutationModels",
+]

astropy/coordinates/new_frames/base.py

@@ -0,0 +1,34 @@
+class BaseCoordinateFrame:
+    def to_icrs(self, coords_self):
+        raise NotImplementedError
+
+    def from_icrs(self, coords_icrs):
+        raise NotImplementedError
+
+
+class HasObstime:
+    @property
+    def obstime(self):
+        return self._obstime
+
+
+class HasEquinox:
+    @property
+    def equinox(self):
+        return self._equinox
+
+
+class HasNutation:
+    @property
+    def nutation(self):
+        return self._nutation
+
+
+class HasObserver:
+    @property
+    def obsgeoloc(self):
+        return self._obsgeoloc
+
+    @property
+    def obsgeovel(self):
+        return self._obsgeovel

astropy/coordinates/new_frames/earth.py

@@ -0,0 +1,160 @@
+from astropy import _erfa as erfa
+from astropy import units as u
+from astropy.coordinates import CartesianRepresentation, SphericalRepresentation
+from astropy.coordinates.builtin_frames.utils import (
+    DEFAULT_OBSTIME,
+    get_cip,
+    get_jd12,
+    prepare_earth_position_vel,
+)
+
+from .base import HasObstime, HasObserver, BaseCoordinateFrame
+
+
+class CIRS(BaseCoordinateFrame, HasObstime):
+    def __init__(self, obstime=DEFAULT_OBSTIME):
+        self._obstime = obstime
+
+    def to_icrs(self, coords_self):
+        srepr = coords_self.represent_as(SphericalRepresentation)
+        cirs_ra = srepr.lon.to_value(u.radian)
+        cirs_dec = srepr.lat.to_value(u.radian)
+
+        # set up the astrometry context for ICRS<->cirs and then convert to
+        # astrometric coordinate direction
+        jd1, jd2 = get_jd12(self.obstime, "tt")
+        x, y, s = get_cip(jd1, jd2)
+        earth_pv, earth_heliocentric = prepare_earth_position_vel(self.obstime)
+        astrom = erfa.apci(jd1, jd2, earth_pv, earth_heliocentric, x, y, s)
+        i_ra, i_dec = erfa.aticq(cirs_ra, cirs_dec, astrom)
+
+        # When there is a distance, apply the parallax/offset to the SSB as the
+        # last step - ensures round-tripping with the icrs_to_cirs transform
+
+        # the distance in intermedrep is *not* a real distance as it does not
+        # include the offset back to the SSB
+        intermedrep = SphericalRepresentation(
+            lat=u.Quantity(i_dec, u.radian, copy=False),
+            lon=u.Quantity(i_ra, u.radian, copy=False),
+            distance=srepr.distance,
+            copy=False,
+        )
+
+        astrom_eb = CartesianRepresentation(
+            astrom["eb"], unit=u.au, xyz_axis=-1, copy=False
+        )
+
+        return intermedrep + astrom_eb
+
+    def from_icrs(self, coords_icrs):
+        jd1, jd2 = get_jd12(self.obstime, "tt")
+        x, y, s = get_cip(jd1, jd2)
+        earth_pv, earth_heliocentric = prepare_earth_position_vel(self.obstime)
+        astrom = erfa.apci(jd1, jd2, earth_pv, earth_heliocentric, x, y, s)
+
+        # When there is a distance, we first offset for parallax to get the
+        # astrometric coordinate direction and *then* run the ERFA transform for
+        # no parallax/PM. This ensures reversibility and is more sensible for
+        # inside solar system objects
+        astrom_eb = CartesianRepresentation(
+            astrom["eb"], unit=u.au, xyz_axis=-1, copy=False
+        )
+        newcart = coords_icrs.represent_as(CartesianRepresentation) - astrom_eb
+
+        srepr = newcart.represent_as(SphericalRepresentation)
+        i_ra = srepr.lon.to_value(u.radian)
+        i_dec = srepr.lat.to_value(u.radian)
+        cirs_ra, cirs_dec = erfa.atciqz(i_ra, i_dec, astrom)
+
+        newrep = SphericalRepresentation(
+            lat=u.Quantity(cirs_dec, u.radian, copy=False),
+            lon=u.Quantity(cirs_ra, u.radian, copy=False),
+            distance=srepr.distance,
+            copy=False,
+        )
+        return newrep
+
+
+class GCRS(BaseCoordinateFrame, HasObstime, HasObserver):
+    def __init__(
+        self,
+        obstime=DEFAULT_OBSTIME,
+        obsgeoloc=CartesianRepresentation([0, 0, 0], unit=u.m),
+        obsgeovel=CartesianRepresentation([0, 0, 0], unit=u.m / u.s),
+    ):
+        self._obstime = obstime
+        self._obsgeoloc = obsgeoloc
+        self._obsgeovel = obsgeovel
+
+    def to_icrs(self, coords_self):
+        srepr = coords_self.represent_as(SphericalRepresentation)
+        gcrs_ra = srepr.lon.to_value(u.radian)
+        gcrs_dec = srepr.lat.to_value(u.radian)
+
+        # set up the astrometry context for ICRS<->GCRS and then convert to BCRS
+        # coordinate direction
+        obs_pv = erfa.pav2pv(
+            self.obsgeoloc.get_xyz(xyz_axis=-1).to_value(u.m),
+            self.obsgeovel.get_xyz(xyz_axis=-1).to_value(u.m / u.s),
+        )
+
+        jd1, jd2 = get_jd12(self.obstime, "tt")
+        earth_pv, earth_heliocentric = prepare_earth_position_vel(self.obstime)
+        astrom = erfa.apcs(jd1, jd2, obs_pv, earth_pv, earth_heliocentric)
+        i_ra, i_dec = erfa.aticq(gcrs_ra, gcrs_dec, astrom)
+
+        # When there is a distance, apply the parallax/offset to the SSB as the
+        # last step - ensures round-tripping with the icrs_to_gcrs transform
+
+        # the distance in intermedrep is *not* a real distance as it does not
+        # include the offset back to the SSB
+        intermedrep = SphericalRepresentation(
+            lat=u.Quantity(i_dec, u.radian, copy=False),
+            lon=u.Quantity(i_ra, u.radian, copy=False),
+            distance=srepr.distance,
+            copy=False,
+        )
+
+        astrom_eb = CartesianRepresentation(
+            astrom["eb"], unit=u.au, xyz_axis=-1, copy=False
+        )
+        newrep = intermedrep + astrom_eb
+
+        return newrep
+
+    def from_icrs(self, coords_icrs):
+        # first set up the astrometry context for ICRS<->GCRS. There are a few steps...
+        # get the position and velocity arrays for the observatory.  Need to
+        # have xyz in last dimension, and pos/vel in one-but-last.
+        # (Note could use np.stack once our minimum numpy version is >=1.10.)
+        obs_pv = erfa.pav2pv(
+            self.obsgeoloc.get_xyz(xyz_axis=-1).to_value(u.m),
+            self.obsgeovel.get_xyz(xyz_axis=-1).to_value(u.m / u.s),
+        )
+
+        # find the position and velocity of earth
+        jd1, jd2 = get_jd12(self.obstime, "tt")
+        earth_pv, earth_heliocentric = prepare_earth_position_vel(self.obstime)
+        astrom = erfa.apcs(jd1, jd2, obs_pv, earth_pv, earth_heliocentric)
+
+        # When there is a distance,  we first offset for parallax to get the
+        # BCRS coordinate direction and *then* run the ERFA transform for no
+        # parallax/PM. This ensures reversibility and is more sensible for
+        # inside solar system objects
+        astrom_eb = CartesianRepresentation(
+            astrom["eb"], unit=u.au, xyz_axis=-1, copy=False
+        )
+        newcart = coords_icrs.represent_as(CartesianRepresentation) - astrom_eb
+
+        srepr = newcart.represent_as(SphericalRepresentation)
+        i_ra = srepr.lon.to_value(u.radian)
+        i_dec = srepr.lat.to_value(u.radian)
+        gcrs_ra, gcrs_dec = erfa.atciqz(i_ra, i_dec, astrom)
+
+        newrep = SphericalRepresentation(
+            lat=u.Quantity(gcrs_dec, u.radian, copy=False),
+            lon=u.Quantity(gcrs_ra, u.radian, copy=False),
+            distance=srepr.distance,
+            copy=False,
+        )
+        return newrep

astropy/coordinates/new_frames/ecliptic.py

@@ -0,0 +1,110 @@
+from astropy import _erfa as erfa
+from astropy import units as u
+from astropy.coordinates.matrix_utilities import matrix_product, rotation_matrix
+from astropy.coordinates.builtin_frames.utils import (
+    EQUINOX_J2000,
+    DEFAULT_OBSTIME,
+    get_jd12,
+)
+
+from .base import HasEquinox, HasObstime, HasNutation, BaseCoordinateFrame
+from .enums import NutationModels
+from .helpers import apply_affine
+from .equatorial import HCRS
+
+
+def _mean_ecliptic_rotation_matrix(equinox):
+    jd1, jd2 = get_jd12(equinox, "tt")
+    rmat = erfa.ecm06(jd1, jd2)
+    return rmat
+
+
+def _true_ecliptic_rotation_matrix(equinox):
+    # This code calls pnm06a from ERFA, which retrieves the precession
+    # matrix (including frame bias) according to the IAU 2006 model, and
+    # including the nutation. This family of systems is less popular
+    # than the "mean" ecliptic ones, obtained using ecm06
+    # (see https://github.com/astropy/astropy/pull/6508).
+    jd1, jd2 = get_jd12(equinox, "tt")
+    rnpb = erfa.pnm06a(jd1, jd2)
+    obl = erfa.obl06(jd1, jd2) * u.radian
+    return matrix_product(rotation_matrix(obl, "x"), rnpb)
+
+
+class BarycentricEcliptic(BaseCoordinateFrame, HasEquinox, HasNutation):
+    def __init__(self, equinox=EQUINOX_J2000, nutation=NutationModels.MEAN):
+        self._equinox = equinox
+        self._nutation = NutationModels(nutation)
+
+    def to_icrs(self, coords_self):
+        if self.nutation is NutationModels.MEAN:
+            rotation = _mean_ecliptic_rotation_matrix(self.equinox).T
+        elif self.nutation is NutationModels.TRUE:
+            rotation = _true_ecliptic_rotation_matrix(self.equinox).T
+
+        return apply_affine(coords_self, rotation)
+
+    def from_icrs(self, coords_icrs):
+        if self.nutation is NutationModels.MEAN:
+            rotation = _mean_ecliptic_rotation_matrix(self.equinox)
+        elif self.nutation is NutationModels.TRUE:
+            rotation = _true_ecliptic_rotation_matrix(self.equinox)
+
+        return apply_affine(coords_icrs, rotation)
+
+
+class BarycentricMeanEcliptic(BarycentricEcliptic):
+    def __init__(self, equinox=EQUINOX_J2000):
+        super().__init__(equinox=equinox, nutation=NutationModels.MEAN)
+
+
+class BarycentricTrueEcliptic(BarycentricEcliptic):
+    def __init__(self, equinox=EQUINOX_J2000):
+        super().__init__(equinox=equinox, nutation=NutationModels.TRUE)
+
+
+class HeliocentricEcliptic(BaseCoordinateFrame, HasObstime, HasEquinox, HasNutation):
+    def __init__(
+        self,
+        obstime=DEFAULT_OBSTIME,
+        equinox=EQUINOX_J2000,
+        nutation=NutationModels.MEAN,
+    ):
+        self._obstime = obstime
+        self._equinox = equinox
+        self._nutation = NutationModels(nutation)
+
+    def to_icrs(self, coords_self):
+        if self.nutation is NutationModels.MEAN:
+            coords_heliocentric_equatorial = BarycentricMeanEcliptic(
+                self.equinox
+            ).to_icrs(coords_self)
+        elif self.nutation is NutationModels.TRUE:
+            coords_heliocentric_equatorial = BarycentricTrueEcliptic(
+                self.equinox
+            ).to_icrs(coords_self)
+
+        return HCRS(self.obstime).to_icrs(coords_heliocentric_equatorial)
+
+    def from_icrs(self, coords_icrs):
+        coords_heliocentric_equatorial = HCRS(self.obstime).from_icrs(coords_icrs)
+        if self.nutation is NutationModels.MEAN:
+            coords_self = BarycentricMeanEcliptic(self.equinox).from_icrs(
+                coords_heliocentric_equatorial
+            )
+        elif self.nutation is NutationModels.TRUE:
+            coords_self = BarycentricTrueEcliptic(self.equinox).from_icrs(
+                coords_heliocentric_equatorial
+            )
+
+        return coords_self
+
+
+class HeliocentricMeanEcliptic(HeliocentricEcliptic):
+    def __init__(self, obstime=DEFAULT_OBSTIME, equinox=EQUINOX_J2000):
+        super().__init__(obstime=obstime, equinox=equinox, nutation=NutationModels.MEAN)
+
+
+class HeliocentricTrueEcliptic(HeliocentricEcliptic):
+    def __init__(self, obstime=DEFAULT_OBSTIME, equinox=EQUINOX_J2000):
+        super().__init__(obstime=obstime, equinox=equinox, nutation=NutationModels.TRUE)

astropy/coordinates/new_frames/enums.py

@@ -0,0 +1,6 @@
+from enum import Enum, auto
+
+
+class NutationModels(Enum):
+    MEAN = auto()
+    TRUE = auto()

astropy/coordinates/new_frames/equatorial.py

@@ -0,0 +1,39 @@
+from astropy.coordinates.builtin_frames.utils import DEFAULT_OBSTIME
+
+from .base import HasObstime, BaseCoordinateFrame
+from .helpers import apply_affine
+
+
+class ICRS(BaseCoordinateFrame):
+    def to_icrs(self, coords_self):
+        return coords_self
+
+    def from_icrs(self, coords_icrs):
+        return coords_icrs
+
+
+class HCRS(BaseCoordinateFrame, HasObstime):
+    def __init__(self, obstime=DEFAULT_OBSTIME):
+        self._obstime = obstime
+
+    def to_icrs(self, coords_self):
+        if coords_self.differentials:
+            raise NotImplementedError
+        else:
+            # TODO: This should be cached! Most of the times it will be DEFAULT_OBSTIME
+            from astropy.coordinates.solar_system import get_body_barycentric
+
+            bary_sun_pos = get_body_barycentric("sun", self.obstime)
+
+        return apply_affine(coords_self, None, bary_sun_pos)
+
+    def from_icrs(self, coords_icrs):
+        if coords_icrs.differentials:
+            raise NotImplementedError
+        else:
+            # TODO: This should be cached! Most of the times it will be DEFAULT_OBSTIME
+            from astropy.coordinates.solar_system import get_body_barycentric
+
+            bary_sun_pos = -get_body_barycentric("sun", self.obstime)
+
+        return apply_affine(coords_icrs, None, bary_sun_pos)

astropy/coordinates/new_frames/helpers.py

@@ -0,0 +1,14 @@
+from astropy.coordinates import CartesianRepresentation
+
+
+def apply_affine(coords, left_rotation=None, offset=None, right_rotation=None):
+    if left_rotation is not None:
+        coords = coords.represent_as(CartesianRepresentation).transform(left_rotation)
+
+    if offset is not None:
+        coords = coords + offset
+
+    if right_rotation is not None:
+        coords = coords.represent_as(CartesianRepresentation).transform(right_rotation)
+
+    return coords