Ellipsoids

.github/workflows/testsuite.yaml

@@ -15,10 +15,10 @@ jobs:
       - uses: actions/checkout@v3
         with:
          fetch-depth: 1
-      - name: Set up Python 3.8
+      - name: Set up Python 3.9
         uses: actions/setup-python@v4
         with:
-          python-version: 3.8
+          python-version: 3.9
       - name: Setup Environment
         run: |
           pip install pre-commit
@@ -36,7 +36,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v4
         with:
-          python-version: 3.8
+          python-version: 3.9
       - name: Get current date
         id: date
         run: echo "::set-output name=date::$(date +'%Y-%m-%d')"
@@ -63,7 +63,7 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest]
-        python-version: [3.7, 3.8, 3.9]
+        python-version: ["3.9", "3.10", "3.11"]
     steps:
       - uses: actions/checkout@v3
         with:

.pre-commit-config.yaml

@@ -26,6 +26,6 @@ repos:
       rev: 22.3.0
       hooks:
           - id: black
-            language_version: python3.8
+            language_version: python3.9
             args:
               - --line-length=90

CHANGELOG.md

@@ -1,5 +1,13 @@
 # Changelog
 
+## [0.2.2] -- 2022-03-19
+
+## Added
+- Support for ellipsoid models for selenodetic coordinates (non-spherical)
+
+## Deprecated
+- Removed support for Python < 3.9. Astropy >= 6.0 is required for ellipsoid support
+
 ## [0.2.0] -- 2022-10-12
 
 ## Fixed

lunarsky/moon.py

@@ -4,28 +4,96 @@
 from astropy.units.quantity import QuantityInfoBase
 from astropy.coordinates.angles import Longitude, Latitude
 from astropy.coordinates.earth import GeodeticLocation
+from astropy.coordinates.representation.geodetic import BaseGeodeticRepresentation
 from astropy.coordinates.representation import (
     CartesianRepresentation,
-    SphericalRepresentation,
 )
 from astropy.coordinates.attributes import Attribute
 
-from .spice_utils import remove_topo, LUNAR_RADIUS
+from .spice_utils import remove_topo
+
+LUNAR_RADIUS = 1737.1e3  # m
 
 __all__ = ["MoonLocation", "MoonLocationAttribute"]
 
 
+class SPHERESelenodeticRepresentation(BaseGeodeticRepresentation):
+    """Lunar ellipsoid as a sphere
+
+    Radius defined by lunarsky.spice_utils.LUNAR_RADIUS
+    """
+
+    _equatorial_radius = LUNAR_RADIUS * u.m
+    _flattening = 0.0
+
+
+class GSFCSelenodeticRepresentation(BaseGeodeticRepresentation):
+    """Lunar ellipsoid from NASA/GSFC "Planetary Fact Sheet"
+
+    https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html
+    """
+
+    _equatorial_radius = 1738.1e3 * u.m
+    _flattening = 0.0012
+
+
+class GRAIL23SelenodeticRepresentation(BaseGeodeticRepresentation):
+    """Lunar ellipsoid defined by gravimetry of GRAIL data.
+
+    https://doi.org/10.1007/s40328-023-00415-w
+    """
+
+    _equatorial_radius = 1737576.6 * u.m
+    _flattening = 0.000305
+
+
+class CE1LAM10SelenodeticRepresentation(BaseGeodeticRepresentation):
+    """Lunar ellipsoid from Chang'e 1 laser altimetry.
+
+    Rotation ellipsoid = CE-1-LAM-GEO
+
+    https://doi.org/10.1007/s11430-010-4060-6
+    """
+
+    _equatorial_radius = 1737.632 * u.km
+    _flattening = 1 / 973.463
+
+
+# Define reference ellipsoids
+SELENOIDS = {
+    "SPHERE": SPHERESelenodeticRepresentation,
+    "GSFC": GSFCSelenodeticRepresentation,
+    "GRAIL23": GRAIL23SelenodeticRepresentation,
+    "CE-1-LAM-GEO": CE1LAM10SelenodeticRepresentation,
+}
+
+
+class SelenodeticLocation(GeodeticLocation):
+    """Rename GeodeticLocation class for clarity"""
+
+
+def _check_ellipsoid(ellipsoid=None, default="SPHERE"):
+    """Defaulting lunar ellipsoid"""
+    if ellipsoid is None:
+        ellipsoid = default
+    if ellipsoid not in SELENOIDS:
+        raise ValueError(f"Ellipsoid {ellipsoid} not among known ones ({SELENOIDS})")
+    return ellipsoid
+
+
 class MoonLocationInfo(QuantityInfoBase):
     """
     Container for meta information like name, description, format.  This is
     required when the object is used as a mixin column within a table, but can
     be used as a general way to store meta information.
     """
 
-    _represent_as_dict_attrs = ("x", "y", "z")
+    _represent_as_dict_attrs = ("x", "y", "z", "ellipsoid")
 
     def _construct_from_dict(self, map):
+        ellipsoid = map.pop("ellipsoid")
         out = self._parent_cls(**map)
+        out.ellipsoid = ellipsoid
         return out
 
     def new_like(self, cols, length, metadata_conflicts="warn", name=None):
@@ -53,7 +121,7 @@
             Empty instance of this class consistent with ``cols``
         """
         # Very similar to QuantityInfo.new_like, but the creation of the
-        # map is different enough that this needs its own rouinte.
+        # map is different enough that this needs its own routine.
         # Get merged info attributes shape, dtype, format, description.
         attrs = self.merge_cols_attributes(
             cols, metadata_conflicts, name, ("meta", "format", "description")
@@ -93,8 +161,10 @@
     This class uses the ME frame.
 
     Positions may be defined in Cartesian (x, y, z) coordinates with respect to the
-    center of mass of the Moon, or in ``selenodetic'' coordinates (longitude, latitude).
-    In selenodetic coordinates, positions are on the surface exactly.
+    center of mass of the Moon, or in ``selenodetic'' coordinates (longitude, latitude, height).
+
+    Selenodetic coordinates are defined with respect to a reference ellipsoid. The default is a
+    sphere of radius 1731.1e3 km, but other ellipsoids are available. See lunarsky.moon.SELENOIDS.
 
     See:
         "A Standardized Lunar Coordinate System for the Lunar Reconnaissance
@@ -110,11 +180,10 @@
     property.
     """
 
+    _ellipsoid = "SPHERE"
     _location_dtype = np.dtype({"names": ["x", "y", "z"], "formats": [np.float64] * 3})
     _array_dtype = np.dtype((np.float64, (3,)))
 
-    _lunar_radius = LUNAR_RADIUS
-
     # Manage the set of defined ephemerides.
     # Class attributes only
     _inuse_stat_ids = []
@@ -170,7 +239,9 @@
             raise ValueError("Too many unique MoonLocation objects open at once.")
 
         for llh in llh_arr:
-            lonlatheight = "_".join(["{:.4f}".format(ll) for ll in llh])
+            lonlatheight = "_".join(
+                ["{:.4f}".format(ll) for ll in llh] + [inst._ellipsoid]
+            )
             if lonlatheight not in cls._existing_locs:
                 new_stat_id = cls._avail_stat_ids.pop()
                 cls._existing_locs.append(lonlatheight)
@@ -248,7 +319,7 @@
         return inst
 
     @classmethod
-    def from_selenodetic(cls, lon, lat, height=0.0):
+    def from_selenodetic(cls, lon, lat, height=0.0, ellipsoid=None):
         """
         Location on the Moon, from latitude and longitude.
 
@@ -264,6 +335,10 @@
             Height above reference sphere (if float, in meters; default: 0).
             The reference sphere is a sphere of radius 1737.1 kilometers,
             from the center of mass of the Moon.
+        ellipsoid : str, optional
+            Name of the reference ellipsoid to use (default: 'SPHERE').
+            Available ellipsoids are:  'SPHERE', 'GRAIL23', 'CE-1-LAM-GEO'.
+            See docstrings for classes in ELLIPSOIDS dictionary for references.
 
         Raises
         ------
@@ -281,7 +356,11 @@
         relative to a prime meridian, which is itself given by
         the mean position of the "sub-Earth" point on the lunar surface.
 
+        For the conversion to selenocentric coordinates, the ERFA routine
+        ``gd2gce`` is used.  See https://github.com/liberfa/erfa
+
         """
+        ellipsoid = _check_ellipsoid(ellipsoid, default=cls._ellipsoid)
         lon = Longitude(lon, u.degree, copy=False).wrap_at(180 * u.degree)
         lat = Latitude(lat, u.degree, copy=False)
         # don't convert to m by default, so we can use the height unit below.
@@ -294,23 +373,14 @@
                     str(lon.shape), str(lat.shape)
                 )
             )
-        # get selenocentric coordinates. Have to give one-dimensional array.
-
-        lunar_radius = u.Quantity(cls._lunar_radius, u.m, copy=False)
-
-        Npts = lon.size
-        xyz = np.zeros((Npts, 3))
-        xyz[:, 0] = (lunar_radius + height) * np.cos(lat) * np.cos(lon)
-        xyz[:, 1] = (lunar_radius + height) * np.cos(lat) * np.sin(lon)
-        xyz[:, 2] = (lunar_radius + height) * np.sin(lat)
 
-        xyz = np.squeeze(xyz)
-
-        self = xyz.ravel().view(cls._location_dtype, cls).reshape(xyz.shape[:-1])
-        self._unit = u.meter
-        inst = self.to(height.unit)
+        # get selenocentric coordinates. Have to give one-dimensional array.
 
-        return inst
+        selenodetic = SELENOIDS[ellipsoid](lon, lat, height, copy=False)
+        xyz = selenodetic.to_cartesian().get_xyz(xyz_axis=-1) << height.unit
+        self = xyz.view(cls._location_dtype, cls).reshape(selenodetic.shape)
+        self.ellipsoid = ellipsoid
+        return self
 
     def __str__(self):
         return self.__repr__()
@@ -351,11 +421,19 @@
         """Convert to selenodetic coordinates."""
         return self.to_selenodetic()
 
-    def to_selenodetic(self):
+    @property
+    def ellipsoid(self):
+        """The default ellipsoid used to convert to selenodetic coordinates."""
+        return self._ellipsoid
+
+    @ellipsoid.setter
+    def ellipsoid(self, ellipsoid):
+        self._ellipsoid = _check_ellipsoid(ellipsoid)
+
+    def to_selenodetic(self, ellipsoid=None):
         """Convert to selenodetic coordinates (lat, lon, height).
 
-        Height is in reference to a sphere with radius `_lunar_radius`,
-        centered at the center of mass.
+        Height is in reference to the ellipsoid.
 
         Returns
         -------
@@ -364,16 +442,15 @@
             `~astropy.coordinates.Latitude`, and `~astropy.units.Quantity`
 
         """
+        ellipsoid = _check_ellipsoid(ellipsoid, default=self.ellipsoid)
         xyz = self.view(self._array_dtype, u.Quantity)
-        lld = CartesianRepresentation(xyz, xyz_axis=-1, copy=False).represent_as(
-            SphericalRepresentation
+        llh = CartesianRepresentation(xyz, xyz_axis=-1, copy=False).represent_as(
+            SELENOIDS[ellipsoid],
         )
-        return GeodeticLocation(
-            Longitude(lld.lon, u.degree, wrap_angle=180.0 * u.degree, copy=False),
-            Latitude(lld.lat, u.degree, copy=False),
-            u.Quantity(
-                lld.distance - (self._lunar_radius * u.meter), self.unit, copy=False
-            ),
+        return SelenodeticLocation(
+            Longitude(llh.lon, u.degree, wrap_angle=180.0 * u.degree, copy=False),
+            Latitude(llh.lat, u.degree, copy=False),
+            u.Quantity(llh.height, self.unit, copy=False),
         )
 
     @property
@@ -476,6 +553,8 @@
 
     def __array_finalize__(self, obj):
         super().__array_finalize__(obj)
+        if hasattr(obj, "_ellipsoid"):
+            self._ellipsoid = obj._ellipsoid
 
     def __len__(self):
         if self.shape == ():

lunarsky/spice_utils.py

@@ -12,13 +12,9 @@
 from .data import DATA_PATH
 
 _J2000 = Time("J2000")
-LUNAR_RADIUS = 1737.1e3  # m
 
 TEMPORARY_KERNEL_DIR = tempfile.TemporaryDirectory()
 
-# TODO --> Look for updated kernels and coordinate definitions.
-#          Use https://naif.jpl.nasa.gov/pub/naif/MER/misc/pck/pck00008.tpc for lunar radii
-
 
 def check_is_loaded(search):
     """
@@ -73,18 +69,16 @@ def furnish_kernels():
     return kernel_paths
 
 
-def lunar_surface_ephem(latitude, longitude, station_num=98):
+def lunar_surface_ephem(pos_x, pos_y, pos_z, station_num=98):
     """
     Make an SPK for the point on the lunar surface
 
     Creates a temporary file and furnishes from that.
 
     Parameters
     ----------
-    latitude: float
-        Mean-Earth frame selenodetic latitude in degrees.
-    longitude: float
-        Mean-Earth frame selenodetic longitude in degrees.
+    pos_x, pos_y, pos_z: float
+        MCMF frame cartesian position in km
     station_num: int
         Station number
 
@@ -95,15 +89,10 @@ def lunar_surface_ephem(latitude, longitude, station_num=98):
     """
     point_id = 301000 + station_num
 
-    lat = np.radians(latitude)
-    lon = np.radians(longitude)
     ets = np.array([spice.str2et("1950-01-01"), spice.str2et("2150-01-01")])
-    pos_mcmf = spice.latrec(
-        LUNAR_RADIUS / 1e3, lon, lat
-    )  # TODO Use MoonLocation instead?
 
     states = np.zeros((len(ets), 6))
-    states[:, :3] = np.repeat(pos_mcmf[None, :], len(ets), axis=0)
+    states[:, :3] = np.repeat([[pos_x], [pos_y], [pos_z]], len(ets), axis=1).T
 
     center = 301
     frame = "MOON_ME"