Add BaseBodycentricRepresentation

astropy/coordinates/representation/__init__.py

@@ -8,6 +8,7 @@
 from .cylindrical import CylindricalRepresentation, CylindricalDifferential
 from .geodetic import (
     BaseGeodeticRepresentation,
+    BaseBodycentricRepresentation,
     WGS84GeodeticRepresentation,
     WGS72GeodeticRepresentation,
     GRS80GeodeticRepresentation,
@@ -56,6 +57,7 @@
     "CylindricalDifferential",
     "PhysicsSphericalDifferential",
     "BaseGeodeticRepresentation",
+    "BaseBodycentricRepresentation",
     "WGS84GeodeticRepresentation",
     "WGS72GeodeticRepresentation",
     "GRS80GeodeticRepresentation",

astropy/coordinates/representation/geodetic.py

@@ -1,5 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 
+import numpy as np
 import erfa
 
 from astropy import units as u
@@ -45,6 +46,7 @@
     to quantities holding correct values (with units of length and dimensionless,
     respectively), or alternatively an ``_ellipsoid`` attribute to the relevant ERFA
     index (as passed in to `erfa.eform`).
+    Longitudes are east positive and span from 0 to 360 degrees by default.
     """
 
     attr_classes = {"lon": Longitude, "lat": Latitude, "height": u.Quantity}
@@ -94,10 +96,93 @@
         Converts 3D rectangular cartesian coordinates (assumed geocentric) to
         geodetic coordinates.
         """
+        # Compute geodetic/planetodetic angles
         lon, lat, height = erfa.gc2gde(
             cls._equatorial_radius, cls._flattening, cart.get_xyz(xyz_axis=-1)
         )
-        return cls(lon, lat, height, copy=False)
+        return cls(
+            lon,
+            lat,
+            height,
+            copy=False,
+        )
+
+
+@format_doc(geodetic_base_doc)
+class BaseBodycentricRepresentation(BaseRepresentation):
+    """Representation of points in bodycentric 3D coordinates.
+
+    Subclasses need to set attributes ``_equatorial_radius`` and ``_flattening``
+    to quantities holding correct values (with units of length and dimensionless,
+    respectively).
+    Longitudes are east positive and span from 0 to 360 degrees by default.
+    """
+
+    attr_classes = {"lon": Longitude, "lat": Latitude, "height": u.Quantity}
+
+    def __init_subclass__(cls, **kwargs):
+        if (
+            "_equatorial_radius" not in cls.__dict__
+            or "_flattening" not in cls.__dict__
+        ):
+            raise AttributeError(
+                f"{cls.__name__} requires '_equatorial_radius' and '_flattening'."
+            )
+        super().__init_subclass__(**kwargs)
+
+    def __init__(self, lon, lat=None, height=None, copy=True):
+        if height is None and not isinstance(lon, self.__class__):
+            height = 0 << u.m
+
+        super().__init__(lon, lat, height, copy=copy)
+        if not self.height.unit.is_equivalent(u.m):
+            raise u.UnitTypeError(
+                f"{self.__class__.__name__} requires height with units of length."
+            )
+
+    def to_cartesian(self):
+        """
+        Converts bodycentric coordinates to 3D rectangular (geocentric)
+        cartesian coordinates.
+        """
+        coslat = np.cos(self.lat)
+        sinlat = np.sin(self.lat)
+        coslon = np.cos(self.lon)
+        sinlon = np.sin(self.lon)
+        x_spheroid = self._equatorial_radius * coslat * coslon
+        y_spheroid = self._equatorial_radius * coslat * sinlon
+        z_spheroid = self._equatorial_radius * (1 - self._flattening) * sinlat
+        r = (
+            self._equatorial_radius
+            * np.sqrt(coslat**2 + ((1 - self._flattening) * sinlat) ** 2)
+            + self.height
+        )
+        x = r * coslon * coslat
+        y = r * sinlon * coslat
+        z = r * sinlat
+        return CartesianRepresentation(x=x, y=y, z=z, copy=False)
+
+    @classmethod
+    def from_cartesian(cls, cart):
+        """
+        Converts 3D rectangular cartesian coordinates (assumed geocentric) to
+        bodycentric coordinates.
+        """
+        # Compute bodycentric latitude
+        p = np.hypot(cart.x, cart.y)
+        d = np.hypot(p, cart.z)
+        lat = np.arctan2(cart.z, p)
+        p_spheroid = cls._equatorial_radius * np.cos(lat)
+        z_spheroid = (cls._equatorial_radius * (1 - cls._flattening)) * np.sin(lat)
+        r_spheroid = np.hypot(p_spheroid, z_spheroid)
+        height = d - r_spheroid
+        lon = np.arctan2(cart.y, cart.x)
+        return cls(
+            lon,
+            lat,
+            height,
+            copy=False,
+        )
 
 
 @format_doc(geodetic_base_doc)

astropy/coordinates/tests/test_geodetic_representations.py

@@ -8,6 +8,7 @@
     CartesianRepresentation,
     REPRESENTATION_CLASSES,
     BaseGeodeticRepresentation,
+    BaseBodycentricRepresentation,
     GRS80GeodeticRepresentation,
     WGS72GeodeticRepresentation,
     WGS84GeodeticRepresentation,
@@ -29,8 +30,46 @@ class CustomGeodetic(BaseGeodeticRepresentation):
     _equatorial_radius = 4000000.0 * u.m
 
 
+class CustomSphericGeodetic(BaseGeodeticRepresentation):
+    _flattening = 0.0
+    _equatorial_radius = 4000000.0 * u.m
+
+
+class CustomSphericBodycentric(BaseBodycentricRepresentation):
+    _flattening = 0.0
+    _equatorial_radius = 4000000.0 * u.m
+
+
+class IAUMARS2000GeodeticRepresentation(BaseGeodeticRepresentation):
+    _equatorial_radius = 3396190.0 * u.m
+    _flattening = 0.5886007555512007 * u.percent
+
+
+class IAUMARS2000BodycentricRepresentation(BaseBodycentricRepresentation):
+    _equatorial_radius = 3396190.0 * u.m
+    _flattening = 0.5886007555512007 * u.percent
+
+
+def test_geodetic_bodycentric_equivalence_spherical_bodies():
+    initial_cartesian = CartesianRepresentation(
+        x=[1, 3000.0] * u.km, y=[7000.0, 4.0] * u.km, z=[5.0, 6000.0] * u.km
+    )
+
+    gd_transformed = CustomSphericGeodetic.from_representation(initial_cartesian)
+    bc_transformed = CustomSphericBodycentric.from_representation(initial_cartesian)
+    assert_quantity_allclose(gd_transformed.lon, bc_transformed.lon)
+    assert_quantity_allclose(gd_transformed.lat, bc_transformed.lat)
+    assert_quantity_allclose(gd_transformed.height, bc_transformed.height)
+
+
 @pytest.mark.parametrize(
-    "geodeticrepresentation", [CustomGeodetic, WGS84GeodeticRepresentation]
+    "geodeticrepresentation",
+    [
+        CustomGeodetic,
+        WGS84GeodeticRepresentation,
+        IAUMARS2000GeodeticRepresentation,
+        IAUMARS2000BodycentricRepresentation,
+    ],
 )
 def test_cartesian_geodetic_roundtrip(geodeticrepresentation):
     # Test array-valued input in the process.
@@ -48,7 +87,13 @@ def test_cartesian_geodetic_roundtrip(geodeticrepresentation):
 
 
 @pytest.mark.parametrize(
-    "geodeticrepresentation", [CustomGeodetic, WGS84GeodeticRepresentation]
+    "geodeticrepresentation",
+    [
+        CustomGeodetic,
+        WGS84GeodeticRepresentation,
+        IAUMARS2000GeodeticRepresentation,
+        IAUMARS2000BodycentricRepresentation,
+    ],
 )
 def test_geodetic_cartesian_roundtrip(geodeticrepresentation):
     initial_geodetic = geodeticrepresentation(
@@ -122,41 +167,69 @@ def test_geodetic_to_geocentric():
     vvd(xyz[2], -3040908.6861467111, 1e-7, "eraGd2gc", "2/3", status)
 
 
-def test_default_height_is_zero():
-    gd = WGS84GeodeticRepresentation(10 * u.deg, 20 * u.deg)
+@pytest.mark.parametrize(
+    "representation",
+    [
+        WGS84GeodeticRepresentation,
+        IAUMARS2000BodycentricRepresentation,
+    ],
+)
+def test_default_height_is_zero(representation):
+    gd = representation(10 * u.deg, 20 * u.deg)
     assert gd.lon == 10 * u.deg
     assert gd.lat == 20 * u.deg
     assert gd.height == 0 * u.m
 
 
-def test_non_angle_error():
-    with pytest.raises(u.UnitTypeError):
-        WGS84GeodeticRepresentation(20 * u.m, 20 * u.deg, 20 * u.m)
+@pytest.mark.parametrize(
+    "representation",
+    [
+        WGS84GeodeticRepresentation,
+        IAUMARS2000BodycentricRepresentation,
+    ],
+)
+def test_non_angle_error(representation):
+    with pytest.raises(u.UnitTypeError, match="require units equivalent to 'rad'"):
+        representation(20 * u.m, 20 * u.deg, 20 * u.m)
 
 
-def test_non_length_error():
+@pytest.mark.parametrize(
+    "representation",
+    [
+        WGS84GeodeticRepresentation,
+        IAUMARS2000BodycentricRepresentation,
+    ],
+)
+def test_non_length_error(representation):
     with pytest.raises(u.UnitTypeError, match="units of length"):
-        WGS84GeodeticRepresentation(10 * u.deg, 20 * u.deg, 30)
+        representation(10 * u.deg, 20 * u.deg, 30)
 
 
-def test_geodetic_subclass_bad_ellipsoid():
+def test_subclass_bad_ellipsoid():
     # Test incomplete initialization.
 
     msg = "module 'erfa' has no attribute 'foo'"
     with pytest.raises(AttributeError, match=msg):
 
-        class InvalidCustomGeodeticEllipsoid(BaseGeodeticRepresentation):
+        class InvalidCustomEllipsoid(BaseGeodeticRepresentation):
             _ellipsoid = "foo"
 
     assert "foo" not in ELLIPSOIDS
-    assert "customgeodeticellipsoiderror" not in REPRESENTATION_CLASSES
+    assert "invalidcustomellipsoid" not in REPRESENTATION_CLASSES
 
 
-def test_geodetic_subclass_missing_equatorial_radius():
-    msg = "MissingCustomGeodeticAttribute requires '_ellipsoid' or '_equatorial_radius' and '_flattening'."
+@pytest.mark.parametrize(
+    "baserepresentation",
+    [
+        BaseGeodeticRepresentation,
+        BaseBodycentricRepresentation,
+    ],
+)
+def test_geodetic_subclass_missing_equatorial_radius(baserepresentation):
+    msg = "'_equatorial_radius' and '_flattening'."
     with pytest.raises(AttributeError, match=msg):
 
-        class MissingCustomGeodeticAttribute(BaseGeodeticRepresentation):
+        class MissingCustomAttribute(baserepresentation):
             _flattening = 0.075 * u.dimensionless_unscaled
 
-    assert "customgeodeticerror" not in REPRESENTATION_CLASSES
+    assert "missingcustomattribute" not in REPRESENTATION_CLASSES

astropy/coordinates/tests/test_representation.py

@@ -1870,8 +1870,9 @@ def test_represent_as(self):
                 # TODO: Converting a CartesianDifferential to a
                 #       RadialDifferential fails, even on `main`
                 continue
-            elif name.endswith("geodetic"):
-                # TODO: Geodetic representations do not have differentials yet
+            elif "geodetic" in name or "bodycentric" in name:
+                # TODO: spheroidal representations (geodetic or bodycentric)
+                # do not have differentials yet
                 continue
             new_rep = rep1.represent_as(
                 REPRESENTATION_CLASSES[name], DIFFERENTIAL_CLASSES[name]

docs/changes/coordinates/14851.feature.rst

@@ -0,0 +1,2 @@
+Add ``BaseBodycentricRepresentation``, a new spheroidal representation for bodycentric
+latitudes and longitudes spanning from 0 to 360 degrees.

docs/coordinates/representations.rst

@@ -32,9 +32,11 @@ The built-in representation classes are:
   (``rho``), azimuthal angle (``phi``), and height (``z``).
 
 
-Astropy also offers a `~astropy.coordinates.BaseGeodeticRepresentation` useful to
+Astropy also offers a `~astropy.coordinates.BaseGeodeticRepresentation` and
+a `~astropy.coordinates.BaseBodycentricRepresentation` useful to
 create specific representations on spheroidal bodies.
-This is used internally for the standard Earth ellipsoids used in
+`~astropy.coordinates.BaseGeodeticRepresentation` is used internally for the standard
+Earth ellipsoids used in
 `~astropy.coordinates.EarthLocation`
 (`~astropy.coordinates.WGS84GeodeticRepresentation`,
 `~astropy.coordinates.WGS72GeodeticRepresentation`, and
@@ -48,6 +50,13 @@ between the vertical to the surface at a specific point of the spheroid and its
 projection onto the equatorial plane.
 The latitude is a value ranging from -90 to 90 degrees, the longitude from 0 to 360
 degrees.
+`~astropy.coordinates.BaseBodycentricRepresentation` is the coordinate representation
+recommended by the Cartographic Coordinates & Rotational Elements Working Group
+(see for example its `2019 report <https://rdcu.be/b32WL>`_): the bodicentric latitude
+and longitude are spherical latitude and longitude originated at the barycenter of the
+body, the height is the distance from the spheroid surface on the line of sight.
+The latitude is a value ranging from -90 to 90 degrees, the longitude from 0 to 360
+degrees.
 
 .. Note::
    For information about using and changing the representation of
@@ -706,11 +715,13 @@ In pseudo-code, this means that a class will look like::
 
 .. _astropy-coordinates-create-geodetic:
 
-Creating Your Own Geodetic Representation
------------------------------------------
+Creating Your Own Geodetic and Bodycentric Representation
+---------------------------------------------------------
 
 If you would like to use geodetic coordinates on planetary bodies other than the Earth,
-you can define a new class that inherits from  `~astropy.coordinates.BaseGeodeticRepresentation`.
+you can define a new class that inherits from
+`~astropy.coordinates.BaseGeodeticRepresentation` or
+`~astropy.coordinates.BaseBodycentricRepresentation`.
 The equatorial radius and flattening must be both assigned via the attributes
 `_equatorial_radius` and `_flattening`.
 
@@ -721,3 +732,11 @@ For example the spheroid describing Mars as in the
 
         _equatorial_radius = 3393400.0 * u.m
         _flattening = 0.518650 * u.percent
+
+The bodycentric coordinate system representing Mars as in the
+`2000 IAU standard <https://doi.org/10.1023/A:1013939327465>`_ could be defined as::
+
+    class IAUMARS2000BodycentricRepresentation(BaseBodycentricRepresentation):
+
+        _equatorial_radius = 3396190.0 * u.m
+        _flattening = 0.5886008 * u.percent