use jplephem and python instead of SPICE toolkit

.github/workflows/testsuite.yaml

@@ -27,8 +27,12 @@ jobs:
           pre-commit install
           pre-commit run -a
   datacache:
-    name: Cache Data
-    runs-on: ubuntu-latest
+    name: Cache Data (${{ matrix.os }})
+    strategy:
+      fail-fast: false
+      matrix:
+        os: [ubuntu-latest, macos-latest]
+    runs-on: ${{ matrix.os }}
     steps:
       - uses: actions/checkout@v3
         with:
@@ -37,19 +41,17 @@ jobs:
         uses: actions/setup-python@v4
         with:
           python-version: 3.9
-      - name: Get current date
-        id: date
-        run: echo "::set-output name=date::$(date +'%Y-%m-%d')"
       - name: cache-kernels
+        id: cache-kernels
         uses: actions/cache@v3
         with:
-          path: ~/ap_cache
-          key: kernel_cache-${{ steps.date.outputs.date }}
-      - name: Download and cache
-        run: |
-          pip install astropy
-          python ./ci/cache_kernels.py save
+          path: lunarsky/data/spk/planets/de430.bsp
+          key: kernel-de430-${{ runner.os }}
+      - name: Download kernels
         if: steps.cache-kernels.outputs.cache-hit != 'true'
+        run: |
+          pip install -e .
+          python -c "from lunarsky.spice_utils import download_big_kernels; download_big_kernels()"
 
   tests:
     env:
@@ -72,18 +74,11 @@ jobs:
         uses: actions/setup-python@v4
         with:
           python-version: ${{ matrix.python-version }}
-      - name: Get current date
-        id: date
-        run: echo "::set-output name=date::$(date +'%Y-%m-%d')"
       - name: load-cache
         uses: actions/cache@v3
         with:
-          path: ~/ap_cache
-          key: kernel_cache-${{ steps.date.outputs.date }}
-      - name: load kernels
-        run: |
-          pip install astropy
-          python ./ci/cache_kernels.py load
+          path: lunarsky/data/spk/planets/de430.bsp
+          key: kernel-de430-${{ runner.os }}
       - name: Install
         run: |
           pip install pytest pytest-cov

.gitignore

@@ -139,3 +139,7 @@ dmypy.json
 
 # Cython debug symbols
 cython_debug/
+
+# Downloaded SPICE kernels (fetched by download_big_kernels / CI cache)
+lunarsky/data/spk/planets/de430.bsp
+lunarsky/data/**/*.lock

CHANGELOG.md

@@ -1,5 +1,14 @@
 # Changelog
 
+## [Unreleased]
+
+## Changed
+- spice_utils now has a PCK file reader to read the moon's planetary constants
+- Transforms are now defined in Python, not using SPICE
+
+## Deprecated
+- Removed dependency on spiceypy for coordinate definitions. Kernels are now handled with jplephem
+
 ## [0.2.6] -- 2024-12-12
 ## Fixed
 - Avoid deprecation warning with numpy.broadcast_shape

README.md

@@ -9,15 +9,13 @@ An extension to `astropy`, providing selenocentric and topocentric reference fra
 for the Moon and transformations of star positions from the ICRS system to these
 frames. This is to describe the sky as observed from the surface of the Moon.
 
-Non-relativistic transformations are calculated using the SPICE toolkit. Relativistic
-corrections (stellar aberration) will be added.
+Non-relativistic transformations are calculated using lunar orientation data from JPL SPICE kernel files, read with `jplephem`.
 
 
 ## Dependencies
 * `numpy`
-* `astropy>=3.0`
+* `astropy>=6.0.0`
 * `jplephem`
-* `spiceypy`
 
 ## Installation
 
@@ -33,6 +31,14 @@ git clone https://github.com/aelanman/lunarsky
 python setup.py install
 ```
 
+Once you've installed lunarsky, you will need to ensure the relevant SPICE kernel files are downloaded
+before running. To do this, run:
+```
+from lunarsky import spice_utils
+spice_utils.download_big_kernels()
+```
+After that, the kernel files will be stored with the package data. This is around 150 MB.
+
 ## Usage
 
 ![mcmf_coords](./docs/figure.png)
@@ -53,7 +59,7 @@ The cartesian axes of the selenocentric system are those of the MCMF frame. Sele
 
 ## Credit
 
-This package makes use of the ``spiceypy`` wrapper [2] for the JPL SPICE Toolkit, produced by the NASA Navigation and Ancillary Information Facility (NAIF) [3] [4]. The transformations are defined using data in kernel files ``pck/moon_pa_de421_1900-2050.bpc``, ``moon_080317.tf``, and ``moon_assoc_me.tf``. These may be found at [the NAIF website](https://naif.jpl.nasa.gov/pub/naif/generic_kernels), and were produced by Nat Bachman (NAIF/JPL) in March 2008. Further information may be found in the comments in these files in the `data` directory.
+This package uses lunar orientation data from SPICE kernel files produced by the NASA Navigation and Ancillary Information Facility (NAIF) [3] [4]. The transformations are defined using data in kernel files ``pck/moon_pa_de421_1900-2050.bpc``, ``moon_080317.tf``, and ``moon_assoc_me.tf``. These may be found at [the NAIF website](https://naif.jpl.nasa.gov/pub/naif/generic_kernels), and were produced by Nat Bachman (NAIF/JPL) in March 2008. Binary kernel files are read using ``jplephem`` [5]. Further information may be found in the comments in the kernel files in the `data` directory.
 
 ## References
 [1]: Ye, Hanlin, et al. "Looking Vector Direction Analysis for the Moon-Based Earth Observation Optical Sensor." IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, vol. 11, no. 11, Nov. 2018, pp. 4488–99. IEEE Xplore, doi:10.1109/JSTARS.2018.2870247.
@@ -63,3 +69,5 @@ This package makes use of the ``spiceypy`` wrapper [2] for the JPL SPICE Toolkit
 [3]: Acton, C.H.; "Ancillary Data Services of NASA's Navigation and Ancillary Information Facility;" Planetary and Space Science, Vol. 44, No. 1, pp. 65-70, 1996.
 
 [4]: Charles Acton, Nathaniel Bachman, Boris Semenov, Edward Wright; A look toward the future in the handling of space science mission geometry; Planetary and Space Science (2017); DOI 10.1016/j.pss.2017.02.013; https://doi.org/10.1016/j.pss.2017.02.013
+
+[5]: Rhodes, Brandon. jplephem: Python module for evaluating JPL ephemerides. https://github.com/brandon-rhodes/python-jplephem

lunarsky/mcmf.py

@@ -17,8 +17,6 @@
 from astropy.coordinates.transformations import FunctionTransformWithFiniteDifference
 from astropy.coordinates.builtin_frames.icrs import ICRS
 
-import spiceypy as spice
-
 __all__ = ["MCMF"]
 
 _J2000 = Time("J2000", scale="tt")
@@ -56,6 +54,7 @@ def moon_location(self):
 
 
 def make_transform(coo, toframe):
+    from .spice_utils import j2000_to_moon_me, body_position
 
     ap_to_spice = {"icrs": ("J2000", 0), "mcmf": ("MOON_ME", 301)}
 
@@ -67,12 +66,12 @@ def make_transform(coo, toframe):
     obstime = toframe.obstime if to_mcmf else coo.obstime
 
     # Make arrays
-    ets = np.atleast_1d((obstime - _J2000).sec)
+    ets = np.atleast_1d((obstime.tdb - _J2000).sec)
     shape_out = np.broadcast_shapes(coo.shape, ets.shape)
 
     coo_cart = coo.cartesian
 
-    mats = np.stack([spice.pxform("J2000", "MOON_ME", et) for et in ets], axis=0)
+    mats = j2000_to_moon_me(ets)
     if not to_mcmf:
         mats = np.linalg.inv(mats)
 
@@ -85,11 +84,8 @@ def make_transform(coo, toframe):
     # If not unitspherical, shift by origin vector before rotating.
     if not is_unitspherical:
         # Make origin vector(s) in coo's frame.
-        orig_posvel = (
-            np.asarray([spice.spkgeo(to_id, et, from_name, from_id)[0] for et in ets])
-            * un.km
-        )
-        coo_cart -= CartesianRepresentation((orig_posvel.T)[:3])
+        orig_pos = body_position(to_id, ets, from_name, from_id) * un.km
+        coo_cart -= CartesianRepresentation(orig_pos.T)
 
     newrepr = coo_cart.transform(mats).reshape(shape_out)
 

lunarsky/moon.py

@@ -10,22 +10,16 @@
 from astropy.coordinates.attributes import Attribute
 from astropy.utils import minversion
 
-from .spice_utils import remove_topo
-
 COPY_IF_NEEDED = None if minversion(np, "2.0") else False
 
-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
+    _equatorial_radius = 1737.1e3 * u.m
     _flattening = 0.0
 
 
@@ -38,7 +32,7 @@ class GSFCSelenodeticRepresentation(BaseGeodeticRepresentation):
     _equatorial_radius = 1738.1e3 * u.m
     _flattening = 0.0012
 
-
+e
 class GRAIL23SelenodeticRepresentation(BaseGeodeticRepresentation):
     """Lunar ellipsoid defined by gravimetry of GRAIL data.
 
@@ -188,16 +182,6 @@ class MoonLocation(u.Quantity):
     _location_dtype = np.dtype({"names": ["x", "y", "z"], "formats": [np.float64] * 3})
     _array_dtype = np.dtype((np.float64, (3,)))
 
-    # Manage the set of defined ephemerides.
-    # Class attributes only
-    _inuse_stat_ids = []
-    _avail_stat_ids = None
-    _existing_locs = []
-    _ref_count = []
-
-    # This instance's station id(s)
-    station_ids = []
-
     info = MoonLocationInfo()
 
     def __new__(cls, *args, **kwargs):
@@ -217,56 +201,6 @@ def __new__(cls, *args, **kwargs):
                 )
         return self
 
-    @classmethod
-    def _set_site_id(cls, inst):
-        """
-        Set the station ID number and manage registry of defined stations.
-        """
-        if inst.isscalar:
-            llh_arr = [
-                (
-                    inst.lon.deg.item(),
-                    inst.lat.deg.item(),
-                    inst.height.to_value("km").item(),
-                )
-            ]
-            ncrds = 1
-        else:
-            llh_arr = zip(inst.lon.deg, inst.lat.deg, inst.height.to_value("km"))
-            ncrds = inst.lon.size
-
-        statids = []
-        if cls._avail_stat_ids is None:
-            cls._avail_stat_ids = list(range(999, 0, -1))
-
-        if len(cls._avail_stat_ids) < ncrds:
-            raise ValueError("Too many unique MoonLocation objects open at once.")
-
-        for llh in llh_arr:
-            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)
-                statids.append(new_stat_id)
-                cls._inuse_stat_ids.append(new_stat_id)
-                cls._ref_count.append(1)
-            else:
-                ind = cls._existing_locs.index(lonlatheight)
-                cls._ref_count[ind] += 1
-                statids.append(cls._inuse_stat_ids[ind])
-        inst.station_ids = statids
-        return inst
-
-    def _set_station_id(self):
-        """
-        Run classmethod for setting station IDs.
-
-        Convenience function used for testing mostly
-        """
-        self.__class__._set_site_id(self)
-
     @classmethod
     def from_selenocentric(cls, x, y, z, unit=None):
         """
@@ -389,35 +323,6 @@ def from_selenodetic(cls, lon, lat, height=0.0, ellipsoid=None):
     def __str__(self):
         return self.__repr__()
 
-    def copy(self):
-        # Necessary to preserve station_ids list
-        return self.__copy__()
-
-    def __copy__(self):
-        # Ensure that the station_ids are copied as well under shallow copy
-        obj = super().copy()
-        obj.station_ids = self.station_ids
-        return obj
-
-    def __del__(self):
-        # Remove this MoonLocation's station_ids from _inuse_stat_ids and
-        # locations from _existing_locs.
-        # Also clear the corresponding frames from spice variable pool.
-        for si, stat_id in enumerate(self.station_ids):
-            try:
-                ind = self.__class__._inuse_stat_ids.index(stat_id)
-            except ValueError:
-                continue
-            count = self.__class__._ref_count[ind]
-            if count <= 1:
-                sid = self.__class__._inuse_stat_ids.pop(ind)
-                self.__class__._existing_locs.pop(ind)
-                self.__class__._avail_stat_ids.insert(0, sid)
-                remove_topo(stat_id)
-                self.__class__._ref_count.pop(ind)
-            else:
-                self.__class__._ref_count[ind] -= 1
-
     @property
     def selenodetic(self):
         """Convert to selenodetic coordinates."""
@@ -464,7 +369,7 @@ def lon(self):
 
     @property
     def lat(self):
-        """Longitude of the location"""
+        """Latitude of the location"""
         return self.selenodetic[1]
 
     @property
@@ -506,7 +411,7 @@ def get_mcmf(self, obstime=None):
     mcmf = property(
         get_mcmf,
         doc="""An `~astropy.coordinates.MCMF` object  with
-                                     for the location of this object at the
+                                     or the location of this object at the
                                      default ``obstime``.""",
     )