From 2d6eb71ea72030a736255efc48b3eac56276ed36 Mon Sep 17 00:00:00 2001 From: aelanman Date: Wed, 7 May 2025 18:15:42 -0400 Subject: [PATCH 01/21] separate kernel file download into a new function; download to package dir instead of astropy cache --- README.md | 8 ++++++ lunarsky/spice_utils.py | 55 ++++++++++++++++++++++++++++++++--------- 2 files changed, 52 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index ee96430..83ccc6b 100644 --- a/README.md +++ b/README.md @@ -33,6 +33,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) diff --git a/lunarsky/spice_utils.py b/lunarsky/spice_utils.py index 254f2d5..fd42fc9 100644 --- a/lunarsky/spice_utils.py +++ b/lunarsky/spice_utils.py @@ -1,7 +1,10 @@ import numpy as np import os +import urllib.request +from filelock import FileLock import tempfile -from astropy.utils.data import download_files_in_parallel +import warnings +from astropy.utils.console import ProgressBar from astropy.coordinates.matrix_utilities import rotation_matrix from astropy.time import Time import astropy.units as unit @@ -46,29 +49,60 @@ def list_kernels(): ktypes.append(dat[1]) return knames, ktypes +def download_big_kernels(): + """ + Download large (~150 MB) spice kernel files that can't be + included in the python package. + + Installs them to the package data directory. + """ + # LSK and DE430 Kernels + knames = ["lsk/naif0012.tls", "spk/planets/de430.bsp"] + _naif_kernel_url = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/" + paths = [] + for kern in knames: + kurl = _naif_kernel_url + kern + kf = os.path.join(DATA_PATH, kern) + paths.append(kf) + if os.path.exists(kf): + continue + with urllib.request.urlopen(kurl) as response: + total_size = int(response.headers.get("Content-Length", 0)) + if total_size > 0: + chunk_size = max(4096, total_size // 20) + + downloaded = 0 + with ProgressBar(total_size // chunk_size + 1, ipython_widget=False) as pbar, open(kf, 'wb') as ofile, FileLock(kf + ".lock"): + while cur_buf := response.read(chunk_size): + downloaded += len(cur_buf) + if total_size > 0: + pbar.update(downloaded // chunk_size) + + return paths def furnish_kernels(): kernel_names = [ "pck/moon_pa_de421_1900-2050.bpc", "fk/satellites/moon_080317.tf", "fk/satellites/moon_assoc_me.tf", + "lsk/naif0012.tls", + "spk/planets/de430.bsp" ] kernel_paths = [os.path.join(DATA_PATH, kn) for kn in kernel_names] + missing = [] for kp in kernel_paths: - spice.furnsh(kp) + if not os.path.exists(kp): + missing.append(os.path.basename(kp)) + if len(missing) > 0: + warnings.warn(f"Missing kernel file(s). Please run " + "lunarsky.spice_utils.download_big_kernels()") + return None - # LSK and DE430 Kernels - knames = ["lsk/naif0012.tls", "spk/planets/de430.bsp"] - _naif_kernel_url = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/" - kurls = [_naif_kernel_url + kname for kname in knames] - paths = download_files_in_parallel(kurls, cache=True, show_progress=False) - for kp in paths: - kernel_paths.append(kp) + for kp in kernel_paths: spice.furnsh(kp) return kernel_paths - def lunar_surface_ephem(pos_x, pos_y, pos_z, station_num=98): """ Make an SPK for the point on the lunar surface @@ -226,5 +260,4 @@ def earth_pos_mcmf(obstimes): earthpos = unit.Quantity(earthpos.T, "km") return MCMF(*earthpos, obstime=obstimes) - KERNEL_PATHS = furnish_kernels() From 755f1bc2d6dbdb57ad6c4ad11bbc7351228dc368 Mon Sep 17 00:00:00 2001 From: aelanman Date: Thu, 5 Jun 2025 12:06:16 -0400 Subject: [PATCH 02/21] build parent dirs for kernel files on download --- lunarsky/moon.py | 4 ++-- lunarsky/spice_utils.py | 1 + setup.py | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/lunarsky/moon.py b/lunarsky/moon.py index 699fd2d..2ee8370 100644 --- a/lunarsky/moon.py +++ b/lunarsky/moon.py @@ -464,7 +464,7 @@ def lon(self): @property def lat(self): - """Longitude of the location""" + """Latitude of the location""" return self.selenodetic[1] @property @@ -506,7 +506,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``.""", ) diff --git a/lunarsky/spice_utils.py b/lunarsky/spice_utils.py index fd42fc9..91076ac 100644 --- a/lunarsky/spice_utils.py +++ b/lunarsky/spice_utils.py @@ -72,6 +72,7 @@ def download_big_kernels(): chunk_size = max(4096, total_size // 20) downloaded = 0 + os.makedirs(os.path.dirname(kf), exist_ok=True) with ProgressBar(total_size // chunk_size + 1, ipython_widget=False) as pbar, open(kf, 'wb') as ofile, FileLock(kf + ".lock"): while cur_buf := response.read(chunk_size): downloaded += len(cur_buf) diff --git a/setup.py b/setup.py index 60461b5..53a7e60 100644 --- a/setup.py +++ b/setup.py @@ -32,7 +32,7 @@ "test_suite": "pytest", "tests_require": ["pytest"], "setup_requires": ["pytest-runner", "setuptools_scm"], - "install_requires": ["numpy>=1.15", "astropy>=6.0.0", "spiceypy", "jplephem"], + "install_requires": ["numpy>=1.15", "astropy>=6.0.0", "spiceypy", "jplephem", "filelock"], "classifiers": [ "Development Status :: 3 - Alpha", "Intended Audience :: Science/Research", From e62a752086bf009f78cb9702a569116c6498dc36 Mon Sep 17 00:00:00 2001 From: aelanman Date: Thu, 5 Jun 2025 12:08:51 -0400 Subject: [PATCH 03/21] actually write downloaded data to disk --- lunarsky/spice_utils.py | 1 + 1 file changed, 1 insertion(+) diff --git a/lunarsky/spice_utils.py b/lunarsky/spice_utils.py index 91076ac..02e4854 100644 --- a/lunarsky/spice_utils.py +++ b/lunarsky/spice_utils.py @@ -78,6 +78,7 @@ def download_big_kernels(): downloaded += len(cur_buf) if total_size > 0: pbar.update(downloaded // chunk_size) + ofile.write(cur_buf) return paths From cc5f606f32aa14c488b93bc47e65086f75a45ec4 Mon Sep 17 00:00:00 2001 From: aelanman Date: Thu, 5 Jun 2025 12:16:38 -0400 Subject: [PATCH 04/21] specify package data --- setup.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/setup.py b/setup.py index 53a7e60..311b24e 100644 --- a/setup.py +++ b/setup.py @@ -29,6 +29,10 @@ "write_to": "lunarsky/version.py", }, "include_package_data": True, + "package_data":{ + "lunarsky.data" : ["lunarsky/data/*"], + "lunarsky.tests": [], + }, "test_suite": "pytest", "tests_require": ["pytest"], "setup_requires": ["pytest-runner", "setuptools_scm"], From bbc9d0d4c039d254299f9241f0f0d64217706948 Mon Sep 17 00:00:00 2001 From: aelanman Date: Thu, 5 Jun 2025 12:22:38 -0400 Subject: [PATCH 05/21] undo --- setup.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/setup.py b/setup.py index 311b24e..53a7e60 100644 --- a/setup.py +++ b/setup.py @@ -29,10 +29,6 @@ "write_to": "lunarsky/version.py", }, "include_package_data": True, - "package_data":{ - "lunarsky.data" : ["lunarsky/data/*"], - "lunarsky.tests": [], - }, "test_suite": "pytest", "tests_require": ["pytest"], "setup_requires": ["pytest-runner", "setuptools_scm"], From 7a9a893ab8911ebd29d7d8653b834bc912571c31 Mon Sep 17 00:00:00 2001 From: aelanman Date: Tue, 17 Mar 2026 12:55:54 -0400 Subject: [PATCH 06/21] use jplephem and python instead of SPICE toolkit --- lunarsky/mcmf.py | 12 +- lunarsky/moon.py | 91 -------- lunarsky/spice_utils.py | 416 +++++++++++++++++++++-------------- lunarsky/tests/conftest.py | 18 -- lunarsky/tests/test_moon.py | 149 ------------- lunarsky/tests/test_spice.py | 149 +++++-------- lunarsky/time.py | 4 +- lunarsky/topo.py | 127 ++++++----- setup.py | 2 +- 9 files changed, 373 insertions(+), 595 deletions(-) diff --git a/lunarsky/mcmf.py b/lunarsky/mcmf.py index e462f6a..f497387 100644 --- a/lunarsky/mcmf.py +++ b/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)} @@ -72,7 +71,7 @@ def make_transform(coo, toframe): 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) diff --git a/lunarsky/moon.py b/lunarsky/moon.py index 699fd2d..a8229ce 100644 --- a/lunarsky/moon.py +++ b/lunarsky/moon.py @@ -10,8 +10,6 @@ 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 @@ -188,16 +186,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 +205,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 +327,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.""" diff --git a/lunarsky/spice_utils.py b/lunarsky/spice_utils.py index 254f2d5..45e749d 100644 --- a/lunarsky/spice_utils.py +++ b/lunarsky/spice_utils.py @@ -1,230 +1,304 @@ import numpy as np import os -import tempfile from astropy.utils.data import download_files_in_parallel from astropy.coordinates.matrix_utilities import rotation_matrix from astropy.time import Time import astropy.units as unit -import spiceypy as spice - from .mcmf import MCMF from .data import DATA_PATH _J2000 = Time("J2000") -TEMPORARY_KERNEL_DIR = tempfile.TemporaryDirectory() +# --------------------- +# Binary PCK reader +# --------------------- + +_PCK_DATA = None -def check_is_loaded(search): +def _read_bpc(filepath): """ - Search the kernel pool variable names for a given string. + Read a DAF-format binary PCK file (Type 2 Chebyshev). + + Returns a dict with segment metadata and coefficient arrays. """ - try: - spice.gnpool(search, 0, 100) - except (spice.support_types.SpiceyError): - return False - return True + from jplephem.daf import DAF + daf = DAF(open(filepath, "rb")) -def list_kernels(): - """ - List loaded kernels. + segments = list(daf.summaries()) + if len(segments) != 1: + raise ValueError(f"Expected 1 segment in binary PCK, got {len(segments)}") + + name, values = segments[0] + data_type = int(values[4]) + start_addr = int(values[5]) + end_addr = int(values[6]) + + if data_type != 2: + raise ValueError(f"Unsupported binary PCK data type {data_type}") + + arr = daf.read_array(start_addr, end_addr) + + init_epoch = arr[-4] + interval = arr[-3] + rsize = int(arr[-2]) + n_records = int(arr[-1]) + n_coeffs = (rsize - 2) // 3 + + data = arr[: n_records * rsize].reshape(n_records, rsize) + + return { + "init_epoch": init_epoch, + "interval": interval, + "n_coeffs": n_coeffs, + "n_records": n_records, + "data": data, + } + + +def _ensure_pck(): + global _PCK_DATA + if _PCK_DATA is None: + _PCK_DATA = _read_bpc( + os.path.join(DATA_PATH, "pck", "moon_pa_de421_1900-2050.bpc") + ) + return _PCK_DATA + + +def _cheby_eval(coeffs, tau): + """Evaluate Chebyshev polynomial at tau in [-1, 1]. coeffs shape: (..., n_coeffs).""" + n = coeffs.shape[-1] + if n == 0: + return np.zeros(coeffs.shape[:-1]) + T_prev = np.ones(coeffs.shape[:-1]) + if n == 1: + return coeffs[..., 0] * T_prev + T_curr = tau + result = coeffs[..., 0] * T_prev + coeffs[..., 1] * T_curr + for i in range(2, n): + T_prev, T_curr = T_curr, 2 * tau * T_curr - T_prev + result = result + coeffs[..., i] * T_curr + return result + + +def _R1(a): + c, s = np.cos(a), np.sin(a) + return np.array([[1, 0, 0], [0, c, s], [0, -s, c]]) - Returns - ------- - list of str - Kernel names (file paths) - list of str - Corresponding kernel types - """ - knames, ktypes = [], [] - for typ in ["spk", "fk", "tk", "pck", "lsk"]: - for ii in range(spice.ktotal(typ)): - dat = spice.kdata(ii, typ) - knames.append(dat[0]) - ktypes.append(dat[1]) - return knames, ktypes - - -def furnish_kernels(): - kernel_names = [ - "pck/moon_pa_de421_1900-2050.bpc", - "fk/satellites/moon_080317.tf", - "fk/satellites/moon_assoc_me.tf", - ] - kernel_paths = [os.path.join(DATA_PATH, kn) for kn in kernel_names] - for kp in kernel_paths: - spice.furnsh(kp) - - # LSK and DE430 Kernels - knames = ["lsk/naif0012.tls", "spk/planets/de430.bsp"] - _naif_kernel_url = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/" - kurls = [_naif_kernel_url + kname for kname in knames] - paths = download_files_in_parallel(kurls, cache=True, show_progress=False) - for kp in paths: - kernel_paths.append(kp) - spice.furnsh(kp) - - return kernel_paths - - -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. +def _R2(a): + c, s = np.cos(a), np.sin(a) + return np.array([[c, 0, -s], [0, 1, 0], [s, 0, c]]) + + +def _R3(a): + c, s = np.cos(a), np.sin(a) + return np.array([[c, s, 0], [-s, c, 0], [0, 0, 1]]) + + +# PA to ME constant rotation from moon_080317.tf +# TKFRAME_31007_ANGLES = (67.92, 78.56, 0.30) arcseconds, AXES = (3, 2, 1) +# SPICE convention: R_PA_to_ME = R1(-a3) @ R2(-a2) @ R3(-a1) +_a1 = np.radians(67.92 / 3600) +_a2 = np.radians(78.56 / 3600) +_a3 = np.radians(0.30 / 3600) +_PA_TO_ME = _R1(-_a3) @ _R2(-_a2) @ _R3(-_a1) + + +def j2000_to_moon_me(ets): + """ + Compute J2000 -> MOON_ME rotation matrices. Parameters ---------- - pos_x, pos_y, pos_z: float - MCMF frame cartesian position in km - station_num: int - Station number + ets : float or array_like + TDB seconds past J2000. Returns ------- - int: - Ephemeris ID number + mats : ndarray, shape (..., 3, 3) """ - point_id = 301000 + station_num + pck = _ensure_pck() + ets = np.atleast_1d(np.asarray(ets, dtype=float)) - ets = np.array([spice.str2et("1950-01-01"), spice.str2et("2150-01-01")]) + record_idx = ((ets - pck["init_epoch"]) / pck["interval"]).astype(int) + record_idx = np.clip(record_idx, 0, pck["n_records"] - 1) - states = np.zeros((len(ets), 6)) - states[:, :3] = np.repeat([[pos_x], [pos_y], [pos_z]], len(ets), axis=1).T + records = pck["data"][record_idx] + mid = records[:, 0] + half = records[:, 1] + nc = pck["n_coeffs"] - center = 301 - frame = "MOON_ME" - degree = 1 + tau = (ets - mid) / half - fname = os.path.join(TEMPORARY_KERNEL_DIR.name, "lunar_points.bsp") - if os.path.exists(fname): - spice.unload(fname) - handle = spice.spkopa(fname) - else: - handle = spice.spkopn(fname, "SPK_FILE", 0) - spice.spkw09( - handle, - point_id, - center, - frame, - ets[0], - ets[-1], - "0", - degree, - len(ets), - states.tolist(), - ets.tolist(), - ) - spice.spkcls(handle) - spice.furnsh(fname) - - return point_id - - -def topo_frame_def(latitude, longitude, station_num=98, moon=True): - """ - Make a list of strings defining a topocentric frame. This can then be loaded - with spiceypy.lmpool. + ra_coeffs = records[:, 2 : 2 + nc] + dec_coeffs = records[:, 2 + nc : 2 + 2 * nc] + w_coeffs = records[:, 2 + 2 * nc : 2 + 3 * nc] + + ra = _cheby_eval(ra_coeffs, tau) + dec = _cheby_eval(dec_coeffs, tau) + w = _cheby_eval(w_coeffs, tau) + + # J2000 -> MOON_PA: R3(W) @ R1(Dec) @ R3(RA) + # Then apply PA -> ME constant rotation + mats = np.zeros((len(ets), 3, 3)) + for i in range(len(ets)): + j2000_to_pa = _R3(w[i]) @ _R1(dec[i]) @ _R3(ra[i]) + mats[i] = _PA_TO_ME @ j2000_to_pa + + return mats + + +def moon_me_to_j2000(ets): + """Transpose of j2000_to_moon_me.""" + return np.swapaxes(j2000_to_moon_me(ets), -2, -1) + + +# --------------------- +# Ephemeris (jplephem) +# --------------------- + +_SPK = None + + +def _ensure_spk(): + global _SPK + if _SPK is None: + from jplephem.spk import SPK + + spk_name = "spk/planets/de430.bsp" + _naif_kernel_url = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/" + kurl = [_naif_kernel_url + spk_name] + paths = download_files_in_parallel(kurl, cache=True, show_progress=False) + _SPK = SPK.open(paths[0]) + return _SPK + + +def _et_to_jd(ets): + """Convert ET (TDB seconds past J2000) to Julian date.""" + return np.asarray(ets) / 86400.0 + 2451545.0 + + +def _pos_ssb_j2000(body_id, ets): """ - if moon: - idnum = 1301000 - station_name = f"LUNAR-TOPO-{station_num}" - relative = "MOON_ME" - else: - # Used in tests only - idnum = 1399000 - station_name = f"EARTH-TOPO-{station_num}" - relative = "ITRF93" + Get position of a solar system body relative to SSB in J2000. - # The DSS stations are built into SPICE, and they number up to 66. - # We will call this station number 98. - idnum += station_num - fm_center_id = idnum - 1000000 + Parameters + ---------- + body_id : int + 0=SSB, 301=Moon, 399=Earth + ets : ndarray + ET values - ecef_to_enu = np.matmul( - rotation_matrix(-longitude, "z", unit="deg"), - rotation_matrix(latitude, "y", unit="deg"), - ).T - # Reorder the axes so that X,Y,Z = E,N,U - ecef_to_enu = ecef_to_enu[[2, 1, 0]] + Returns + ------- + pos : ndarray, shape (N, 3) + Position in km. + """ + if body_id == 0: + return np.zeros((len(ets), 3)) + + spk = _ensure_spk() + jd = _et_to_jd(ets) + + if body_id == 301: + emb = np.asarray(spk[0, 3].compute(jd)).T + moon_emb = np.asarray(spk[3, 301].compute(jd)).T + return emb + moon_emb + elif body_id == 399: + emb = np.asarray(spk[0, 3].compute(jd)).T + earth_emb = np.asarray(spk[3, 399].compute(jd)).T + return emb + earth_emb + else: + raise ValueError(f"Unsupported body ID {body_id}") - mat = " ".join(map("{:.7f}".format, ecef_to_enu.flatten())) - fmt_strs = [ - "FRAME_{1} = {0}", - "FRAME_{0}_NAME = '{1}'", - "FRAME_{0}_CLASS = 4", - "FRAME_{0}_CLASS_ID = {0}", - "FRAME_{0}_CENTER = {2}", - "OBJECT_{2}_FRAME = '{1}'", - "TKFRAME_{0}_RELATIVE = '{3}'", - "TKFRAME_{0}_SPEC = 'MATRIX'", - "TKFRAME_{0}_MATRIX = {4}", - ] +def body_position(target_id, ets, frame, observer_id): + """ + Get position of target relative to observer in given frame. - frame_specs = [ - s.format(idnum, station_name, fm_center_id, relative, mat) for s in fmt_strs - ] + Parameters + ---------- + target_id : int + NAIF body ID of target (0=SSB, 301=Moon, 399=Earth) + ets : array_like + ET values (TDB seconds past J2000) + frame : str + "J2000" or "MOON_ME" + observer_id : int + NAIF body ID of observer - frame_dict = {} + Returns + ------- + pos : ndarray, shape (N, 3) + Position in km. + """ + ets = np.atleast_1d(np.asarray(ets, dtype=float)) - for spec in frame_specs: - k, v = map(str.strip, spec.split("=")) - frame_dict[k] = v + pos_j2000 = _pos_ssb_j2000(target_id, ets) - _pos_ssb_j2000(observer_id, ets) - latlon = ["{:.4f}".format(coord) for coord in [latitude, longitude]] + if frame == "J2000": + return pos_j2000 + elif frame == "MOON_ME": + mats = j2000_to_moon_me(ets) + return np.einsum("nij,nj->ni", mats, pos_j2000) + else: + raise ValueError(f"Unsupported frame {frame}") - return station_name, idnum, frame_dict, latlon +def station_pos_ssb_j2000(pos_me_km, ets): + """ + Get position of a lunar surface station relative to SSB in J2000. -def remove_topo(station_num): - """Remove a lunar station, by number, from variable pool.""" + Parameters + ---------- + pos_me_km : ndarray, shape (3,) + Station position in MOON_ME frame, in km. + ets : ndarray + ET values. - idnum = 1301000 + station_num - fm_center_id = idnum - 1000000 - station_name = f"LUNAR-TOPO-{station_num}" + Returns + ------- + pos : ndarray, shape (N, 3) + Position in km. + """ + ets = np.atleast_1d(np.asarray(ets, dtype=float)) + moon_ssb = _pos_ssb_j2000(301, ets) + me_to_j2000 = moon_me_to_j2000(ets) + station_j2000 = np.einsum("nij,j->ni", me_to_j2000, pos_me_km) + return moon_ssb + station_j2000 - fmt_vars = [ - "FRAME_{1}", - "FRAME_{0}_NAME", - "FRAME_{0}_CLASS", - "FRAME_{0}_CLASS_ID", - "FRAME_{0}_CENTER", - "OBJECT_{2}_FRAME", - "TKFRAME_{0}_RELATIVE", - "TKFRAME_{0}_SPEC", - "TKFRAME_{0}_MATRIX", - ] - frame_vars = [s.format(idnum, station_name, fm_center_id) for s in fmt_vars] +def topo_rotation_matrix(lat_deg, lon_deg): + """ + Compute the MOON_ME -> topocentric (E/N/U) rotation matrix for a surface location. - # Handle a bug in spiceypy for older versions of numpy - if np.str_ is None: - return - for var in frame_vars: - spice.dvpool(var) + Parameters + ---------- + lat_deg, lon_deg : float + Selenodetic latitude and longitude in degrees. - # Ideally, one would also remove the ephemeris data for this station from - # the lunar_points.bsp file. This doesn't seem to be possible. However, - # the ephemeris _should_ be overwritten if the station_id is reused. + Returns + ------- + matrix : ndarray, shape (3, 3) + """ + ecef_to_enu = np.matmul( + rotation_matrix(-lon_deg, "z", unit="deg"), + rotation_matrix(lat_deg, "y", unit="deg"), + ).T + ecef_to_enu = ecef_to_enu[[2, 1, 0]] + return np.asarray(ecef_to_enu) def earth_pos_mcmf(obstimes): """ Get the position of the Earth in the MCMF frame. - Using SPICE. - Used for tests. """ ets = (obstimes - Time("J2000")).sec - earthpos = np.stack( - [spice.spkpos("399", et, "MOON_ME", "None", "301")[0] for et in ets] - ) + earthpos = body_position(399, ets, "MOON_ME", 301) earthpos = unit.Quantity(earthpos.T, "km") return MCMF(*earthpos, obstime=obstimes) - - -KERNEL_PATHS = furnish_kernels() diff --git a/lunarsky/tests/conftest.py b/lunarsky/tests/conftest.py index e79153d..89e3608 100644 --- a/lunarsky/tests/conftest.py +++ b/lunarsky/tests/conftest.py @@ -3,11 +3,7 @@ import pytest import warnings -from lunarsky.spice_utils import remove_topo -from lunarsky.moon import MoonLocation - -# Ignore deprecation warnings from spiceypy @pytest.fixture(autouse=True) def ignore_representation_deprecation(): warnings.filterwarnings( @@ -30,17 +26,3 @@ def grcat(): stars = SkyCoord(ra=ras, dec=decs, unit="deg", frame="icrs") return stars - - -@pytest.fixture -def cleanup_moonlocs(): - # Reset the existing set of station IDs and corresponding kernel vars. - yield - - for sid in MoonLocation._inuse_stat_ids: - remove_topo(sid) - - MoonLocation._inuse_stat_ids = [] - MoonLocation._avail_stat_ids = None - MoonLocation._existing_locs = [] - MoonLocation._ref_count = [] diff --git a/lunarsky/tests/test_moon.py b/lunarsky/tests/test_moon.py index 08b4cbd..a444ade 100644 --- a/lunarsky/tests/test_moon.py +++ b/lunarsky/tests/test_moon.py @@ -1,5 +1,4 @@ import pytest -import copy import numpy as np import astropy.units as unit from astropy.coordinates import ( @@ -107,151 +106,3 @@ def test_moonlocation_attribute(): attr2, boo = mlattr.convert_input(moonloc.mcmf) assert np.all(attr2 == moonloc) - - -def test_moonlocation_copy(): - # Check that station_ids are copied properly - loc0 = MoonLocation.from_selenodetic(lat=["-15d", "25d"], lon=["97d", "0d"]) - loc0._set_station_id() - lcop = loc0.copy() - assert lcop.station_ids == loc0.station_ids - lcop2 = copy.copy(loc0) - assert lcop2.station_ids == loc0.station_ids - - -def test_moonlocation_delete(): - # Check that making multiple instances of the same location raises the _ref_count - # appropriately, and deleting them removes them correctly. - before = copy.copy(MoonLocation._inuse_stat_ids) - if MoonLocation._avail_stat_ids is not None: - exp = before + MoonLocation._avail_stat_ids[-2:][::-1] - else: - exp = before + [1, 2] - - locs = [] - for ii in range(5): - locs.append(MoonLocation.from_selenodetic(lat=["-15d", "25d"], lon=["97d", "0d"])) - locs[-1]._set_station_id() - - check0 = copy.copy(MoonLocation._inuse_stat_ids) - for ii in range(5, 0, -1): - assert MoonLocation._ref_count[-1] == ii - assert MoonLocation._inuse_stat_ids[-1] == check0[-1] - locs.pop(ii - 1) - - check1 = copy.copy(MoonLocation._inuse_stat_ids) - - assert exp == check0 - assert check1 == before - - -@pytest.mark.parametrize("ell", SELENOIDS) -def test_station_ids(ell): - # Check that when a MoonLocations are made, the appropriate station_ids are assigned. - - # Get whatever are already recorded in the class. - orig_statids = copy.copy(MoonLocation._inuse_stat_ids) - - # Random positions with some repeats, in five groups - lonlatheights = [ - [ - (294.67, -67.68, 15.23), - (78.31, 51.60, 16.59), - (335.27, 57.45, 27.32), - ], - [ - (335.27, 57.45, 27.32), - (133.29, -80.63, 25.92), - (323.28, -2.46, 4.74), - (45.94, -25.17, 20.22), - (326.40, 29.21, 9.67), - (197.94, 3.92, 16.34), - ], - [(197.94, 3.92, 16.34)], - [ - (197.94, 3.92, 16.34), - (242.63, -20.86, 4.90), - (70.30, -40.22, 29.91), - (187.70, -85.44, 9.35), - (129.55, 14.11, 0.24), - (140.46, 4.19, 28.03), - (232.22, 56.89, 15.22), - (51.31, -30.92, 11.53), - (55.59, -61.90, 3.98), - ], - [ - (307.76, 32.31, 29.88), - (248.07, 2.89, 22.41), - (349.94, -77.10, 14.27), - (315.41, -72.71, 13.26), - (173.39, 14.43, 3.37), - (183.19, -83.01, 16.13), - (258.53, -56.66, 29.75), - (115.36, -44.42, 11.34), - (270.61, 57.85, 26.16), - (320.07, -75.06, 21.15), - (216.74, 54.75, 1.92), - ], - ] - - all_pos = np.array(sum(lonlatheights, [])) - n_unique = np.unique(all_pos, axis=0).shape[0] - - locs = [] - locstrs = [] - - for gp in lonlatheights: - lons, lats, heights = np.array(gp).T - locs.append( - MoonLocation.from_selenodetic( - lat=lats, lon=lons, height=heights, ellipsoid=ell - ) - ) - locs[-1]._set_station_id() - - # Check that only unique positions got added - added = len(MoonLocation._inuse_stat_ids) - len(orig_statids) - assert added == n_unique - - statids = [loc.station_ids for loc in locs] - - locstrs = [] - for inst in locs: - llhs = [] - if inst.isscalar: - llh_arr = [ - ( - inst.lon.deg.item(), - inst.lat.deg.item(), - inst.height.to_value("km").item(), - ) - ] - else: - llh_arr = zip(inst.lon.deg, inst.lat.deg, inst.height.to_value("km")) - for llh in llh_arr: - llhs.append("_".join(["{:.4f}".format(ll) for ll in llh] + [ell])) - locstrs.append(llhs) - - # Check that saved location strings correspond with station IDs in each instance - for gi, gp in enumerate(statids): - for sti, sid in enumerate(gp): - ind = MoonLocation._inuse_stat_ids.index(sid) - assert locstrs[gi][sti] == MoonLocation._existing_locs[ind] - - -def test_statid_management(cleanup_moonlocs): - lats = np.random.uniform(-90, 90, 1000) - lons = np.random.uniform(0, 360, 1000) - - for ii in range(2): - loc = MoonLocation.from_selenodetic( - lon=lons[ii * 500 : (ii + 1) * 500], - lat=lats[ii * 500 : (ii + 1) * 500], - ) - loc._set_station_id() - assert len(loc.station_ids) == 500 - del loc - - bigloc = MoonLocation.from_selenodetic(lon=lons, lat=lats) - with pytest.raises(ValueError, match="Too many unique MoonLocation"): - bigloc._set_station_id() diff --git a/lunarsky/tests/test_spice.py b/lunarsky/tests/test_spice.py index fcdfc51..6bf0104 100644 --- a/lunarsky/tests/test_spice.py +++ b/lunarsky/tests/test_spice.py @@ -1,120 +1,73 @@ -from astropy.coordinates.baseframe import frame_transform_graph -from astropy.coordinates.transformations import FunctionTransformWithFiniteDifference -from astropy.coordinates import AltAz, ICRS, EarthLocation, Angle -from astropy.utils.data import download_files_in_parallel -import lunarsky +import numpy as np +import os +import pytest import lunarsky.spice_utils as spice_utils -from lunarsky.time import Time -import spiceypy as spice +FIXTURE_PATH = os.path.join(os.path.dirname(__file__), "data", "spice_fixtures.npz") -def test_topo_frame_setup(): - # Check that the frame setup puts all the correct values in the kernel pool. +@pytest.fixture +def fixtures(): + return np.load(FIXTURE_PATH) - latitude, longitude = 30, 25 - name, idnum, frame_dict, latlon = spice_utils.topo_frame_def(latitude, longitude) - frame_strs = ["{}={}".format(k, v) for (k, v) in frame_dict.items()] - spice.lmpool(frame_strs) - for k, v in frame_dict.items(): - N, typecode = spice.dtpool(k) - if typecode == "N": - res = spice.gdpool(k, 0, 100) - if len(res) == 1: - res = [int(res[0])] - if N > 1: - v = [float(it) for it in v.split(" ")] - res = res.tolist() - else: - res = res[0] - v = int(v) - assert v == res - else: - res = spice.gcpool(k, 0, 100)[0] - v = v.replace("'", "") - assert v == res +def test_j2000_to_moon_me(fixtures): + ets = fixtures["ets"] + expected = fixtures["j2000_to_me"] + computed = spice_utils.j2000_to_moon_me(ets) + np.testing.assert_allclose(computed, expected, atol=1e-10) -def test_kernel_paths(): - # Check that the correct kernel files are downloaded - # Need to unhash the file names +def test_body_position_moon_earth(fixtures): + ets = fixtures["ets"] + expected = fixtures["moon_earth_j2000"][:, :3] + computed = spice_utils.body_position(301, ets, "J2000", 399) + np.testing.assert_allclose(computed, expected, atol=1e-6) - assert len(lunarsky.spice_utils.KERNEL_PATHS) == 5 +def test_body_position_moon_ssb(fixtures): + ets = fixtures["ets"] + expected = fixtures["moon_ssb_j2000"][:, :3] + computed = spice_utils.body_position(301, ets, "J2000", 0) + np.testing.assert_allclose(computed, expected, atol=1e-6) -def test_spice_earth(grcat): - # Replace the ICRS->AltAz transform in astropy with one using SPICE. - # Confirm that star positions are the same as with the original transform - # to within the error due to relativistic aberration (~ 21 arcsec) - # DSS-15, 399015 - lat, lon, height = 35.42185, -116.8871951, 973.211 +def test_body_position_ssb_moon_me(fixtures): + ets = fixtures["ets"] + expected = fixtures["ssb_moon_me"][:, :3] + computed = spice_utils.body_position(0, ets, "MOON_ME", 301) + np.testing.assert_allclose(computed, expected, atol=1e-4) - loc = EarthLocation.from_geodetic(lon, lat, height) - t0 = Time.now() - _J2000 = Time("J2000") - altaz = grcat.transform_to(AltAz(location=loc, obstime=t0)) +def test_earth_pos_mcmf(fixtures): + ets = fixtures["ets"] + expected = fixtures["earth_me"] + computed = spice_utils.body_position(399, ets, "MOON_ME", 301) + np.testing.assert_allclose(computed, expected, atol=1e-4) - # Make the Earth topo frame in spice. - framename, idnum, frame_dict, latlon = lunarsky.spice_utils.topo_frame_def( - lat, lon, moon=False - ) - # One more kernel is needed for the ITRF93 frame. - kname = "pck/earth_latest_high_prec.bpc" - _naif_kernel_url = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels" - kurl = [_naif_kernel_url + "/" + kname] - kernpath = download_files_in_parallel( - kurl, - cache=True, - show_progress=False, +def test_topo_rotation_matrix(fixtures): + expected = fixtures["topo_to_me"] + topo_matrix = spice_utils.topo_rotation_matrix( + float(fixtures["topo_station_lat_deg"]), + float(fixtures["topo_station_lon_deg"]), ) - spice.furnsh(kernpath) - - frame_strs = ["{}={}".format(k, v) for (k, v) in frame_dict.items()] - spice.lmpool(frame_strs) - et = (t0 - _J2000).sec - - @frame_transform_graph.transform(FunctionTransformWithFiniteDifference, ICRS, AltAz) - def icrs_to_ecef(icrs_coo, ecef_frame): - mat = spice.pxform("J2000", framename, et) - newrepr = icrs_coo.cartesian.transform(mat) - - return ecef_frame.realize_frame(newrepr) - - tr = frame_transform_graph.get_transform(ICRS, AltAz).transforms - assert tr[0].func.__name__ == "icrs_to_ecef" # Confirm new transform added correctly - - altaz2 = grcat.transform_to(AltAz(location=loc, obstime=t0)) - - # Having done the transform, remove the spice transform from the graph - frame_transform_graph.remove_transform(ICRS, AltAz, None) - - # Compare astropy topocentric coordinate transform results to - assert all(altaz.separation(altaz2) < Angle("25arcsec")) - # NOTE Large deviations in azimuth high altitudes. - + # TOPO -> MOON_ME should be topo_matrix.T (constant) + for i in range(len(expected)): + np.testing.assert_allclose(topo_matrix.T, expected[i], atol=1e-7) -def test_topo_kernel_setup(): - # Tests a single function in topo.py - # Need to clear the kernel pool first + # Verify full precision matrix is orthogonal + np.testing.assert_allclose(topo_matrix @ topo_matrix.T, np.eye(3), atol=1e-14) - spice.clpool() - # Confirm no variables are loaded. - assert not lunarsky.spice_utils.check_is_loaded("*") - for filepath in lunarsky.spice_utils.KERNEL_PATHS: - spice.furnsh(filepath) - loc = lunarsky.MoonLocation.from_selenodetic(lat="30d", lon="20d") - station_name, idnum, frame_specs, latlon = lunarsky.spice_utils.topo_frame_def( - loc.lat.deg, loc.lon.deg, moon=True +def test_topo_to_j2000(fixtures): + ets = fixtures["ets"] + expected_to_j2000 = fixtures["topo_to_j2000"] + topo_matrix = spice_utils.topo_rotation_matrix( + float(fixtures["topo_station_lat_deg"]), + float(fixtures["topo_station_lon_deg"]), ) - statnum = idnum - 1301000 - lunarsky.topo._spice_setup(loc, [statnum]) - try: - assert lunarsky.spice_utils.check_is_loaded("*{}*".format(idnum)) - finally: - lunarsky.spice_utils.remove_topo(statnum) + me_to_j2000 = spice_utils.moon_me_to_j2000(ets) + computed_to_j2000 = np.einsum("nij,jk->nik", me_to_j2000, topo_matrix.T) + np.testing.assert_allclose(computed_to_j2000, expected_to_j2000, atol=1e-10) diff --git a/lunarsky/time.py b/lunarsky/time.py index 1372a77..2461a14 100644 --- a/lunarsky/time.py +++ b/lunarsky/time.py @@ -6,7 +6,7 @@ from astropy.coordinates import EarthLocation, Longitude from .moon import MoonLocation, COPY_IF_NEEDED -import spiceypy as spice +from .spice_utils import moon_me_to_j2000 __all__ = ["Time", "TimeDelta"] @@ -71,7 +71,7 @@ def sidereal_time(self, kind, longitude=None, model=None): # be defined in order to get here. et = np.atleast_1d((self - Time("J2000")).sec) - mats = np.array([spice.pxform("MOON_ME", "J2000", t) for t in et]) + mats = moon_me_to_j2000(et) # Zenith vector zvec = self.location.mcmf.cartesian.xyz diff --git a/lunarsky/topo.py b/lunarsky/topo.py index 3146481..b4c20bd 100644 --- a/lunarsky/topo.py +++ b/lunarsky/topo.py @@ -22,9 +22,13 @@ from .mcmf import MCMF from .moon import MoonLocationAttribute -from .spice_utils import check_is_loaded, topo_frame_def, lunar_surface_ephem - -import spiceypy as spice +from .spice_utils import ( + j2000_to_moon_me, + moon_me_to_j2000, + body_position, + station_pos_ssb_j2000, + topo_rotation_matrix, +) _90DEG = 90 * un.deg __all__ = ["LunarTopo"] @@ -88,22 +92,26 @@ def zen(self): # ----------------- -def _spice_setup(locations, station_ids): - for li, loc in enumerate(np.atleast_1d(locations)): - sid = int(station_ids[li]) - frameloaded = check_is_loaded(f"FRAME_LUNAR-TOPO-{sid}") - if not frameloaded: - lunar_surface_ephem( - loc.x.to_value("km"), - loc.y.to_value("km"), - loc.z.to_value("km"), - station_num=sid, - ) # Furnishes SPK for lunar surface point - station_name, idnum, frame_specs, latlon = topo_frame_def( - loc.lat.deg, loc.lon.deg, moon=True, station_num=sid - ) - frame_strs = ["{}={}".format(k, v) for (k, v) in frame_specs.items()] - spice.lmpool(frame_strs) +def _get_topo_data(location): + """ + Compute topo rotation matrices and MOON_ME positions for each location element. + + Returns + ------- + topo_mats : ndarray, shape (M, 3, 3) + ME-to-topo rotation matrices. + pos_mes : ndarray, shape (M, 3) + Station positions in MOON_ME, in km. + """ + locations = np.atleast_1d(location) + topo_mats = [] + pos_mes = [] + for loc in locations: + topo_mats.append(topo_rotation_matrix(loc.lat.deg, loc.lon.deg)) + pos_mes.append( + [loc.x.to_value("km"), loc.y.to_value("km"), loc.z.to_value("km")] + ) + return np.array(topo_mats), np.array(pos_mes) def make_transform(coo, toframe): @@ -123,29 +131,28 @@ def make_transform(coo, toframe): if location is None: raise ValueError("location must be defined for LunarTopo transformations") - # Initialize station_ids if not defined. - if location.station_ids == []: - location.__class__._set_site_id(location) - # Make arrays ets = (obstime - _J2000).sec - stat_ids = np.asarray(location.station_ids) + topo_mats, pos_mes = _get_topo_data(location) + loc_indices = np.arange(len(topo_mats)) shape_out = np.broadcast_shapes(coo.shape, ets.shape, location.shape) - # Set up SPICE ephemerides and frame details - _spice_setup(location, stat_ids) - - ets_ids = np.atleast_2d(np.stack(np.broadcast_arrays(ets, stat_ids)).T) + ets_locs = np.atleast_2d(np.stack(np.broadcast_arrays(ets, loc_indices)).T) coo_cart = coo.cartesian # Make rotation matrices - mats = np.asarray( - [ - spice.pxform(f"LUNAR-TOPO-{int(stat_id)}", frame_spice_name, et) - for (et, stat_id) in ets_ids - ] - ) + # TOPO->ME is topo_matrix.T (constant); TOPO->J2000 = R(ME->J2000) @ topo_matrix.T + mats_list = [] + for et, loc_idx in ets_locs: + loc_idx = int(loc_idx) + topo_mat = topo_mats[loc_idx] + if frame_spice_name == "MOON_ME": + mats_list.append(topo_mat.T) + else: + me_to_j2000 = moon_me_to_j2000(np.array([et]))[0] + mats_list.append(me_to_j2000 @ topo_mat.T) + mats = np.asarray(mats_list) if totopo: mats = np.linalg.inv(mats) @@ -157,29 +164,35 @@ 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. - if totopo: - origin_id = lambda n: int(frame_id) # MCMF or ICRS frame origin - target_id = lambda n: int(n) + 301000 # Station ID - frame_name = lambda n: frame_spice_name - - else: - origin_id = lambda n: int(n) + 301000 - target_id = lambda n: int(frame_id) - frame_name = lambda n: f"LUNAR-TOPO-{int(n)}" - - orig_posvel = ( - np.asarray( - [ - spice.spkgeo( - target_id(stat_id), et, frame_name(stat_id), origin_id(stat_id) - )[0] - for (et, stat_id) in ets_ids - ] - ) - * un.km - ) - coo_cart -= CartesianRepresentation((orig_posvel.T)[:3]) + orig_pos_list = [] + for et, loc_idx in ets_locs: + loc_idx = int(loc_idx) + pos_me = pos_mes[loc_idx] + topo_mat = topo_mats[loc_idx] + et_arr = np.array([et]) + + if totopo: + # Origin = station position relative to frame origin, in coo's frame + if frame_id == 0: # SSB / ICRS + orig_pos_list.append( + station_pos_ssb_j2000(pos_me, et_arr)[0] + ) + else: # Moon center / MCMF + orig_pos_list.append(pos_me) + else: + # Origin = frame origin relative to station, in coo's frame (topo) + if frame_id == 0: # SSB + station_ssb = station_pos_ssb_j2000(pos_me, et_arr)[0] + ssb_station_j2000 = -station_ssb + # Rotate to topo + j2000_me = j2000_to_moon_me(et_arr)[0] + ssb_station_me = j2000_me @ ssb_station_j2000 + orig_pos_list.append(topo_mat @ ssb_station_me) + else: # Moon center + orig_pos_list.append(-(topo_mat @ pos_me)) + + orig_pos = np.asarray(orig_pos_list) * un.km + coo_cart -= CartesianRepresentation(orig_pos.T) newrepr = coo_cart.transform(mats).reshape(shape_out) diff --git a/setup.py b/setup.py index 60461b5..93bd4a5 100644 --- a/setup.py +++ b/setup.py @@ -32,7 +32,7 @@ "test_suite": "pytest", "tests_require": ["pytest"], "setup_requires": ["pytest-runner", "setuptools_scm"], - "install_requires": ["numpy>=1.15", "astropy>=6.0.0", "spiceypy", "jplephem"], + "install_requires": ["numpy>=1.15", "astropy>=6.0.0", "jplephem"], "classifiers": [ "Development Status :: 3 - Alpha", "Intended Audience :: Science/Research", From 506873675e63efe705e26325d6b45cc8682d1584 Mon Sep 17 00:00:00 2001 From: aelanman Date: Tue, 17 Mar 2026 13:38:31 -0400 Subject: [PATCH 07/21] linting --- lunarsky/tests/test_moon.py | 1 - lunarsky/topo.py | 9 ++------- 2 files changed, 2 insertions(+), 8 deletions(-) diff --git a/lunarsky/tests/test_moon.py b/lunarsky/tests/test_moon.py index a444ade..fd605da 100644 --- a/lunarsky/tests/test_moon.py +++ b/lunarsky/tests/test_moon.py @@ -7,7 +7,6 @@ ) from astropy.tests.helper import quantity_allclose from lunarsky import MoonLocation, MoonLocationAttribute, MCMF -from lunarsky.moon import SELENOIDS class TestsWithObject: diff --git a/lunarsky/topo.py b/lunarsky/topo.py index b4c20bd..27e9a7c 100644 --- a/lunarsky/topo.py +++ b/lunarsky/topo.py @@ -25,7 +25,6 @@ from .spice_utils import ( j2000_to_moon_me, moon_me_to_j2000, - body_position, station_pos_ssb_j2000, topo_rotation_matrix, ) @@ -108,9 +107,7 @@ def _get_topo_data(location): pos_mes = [] for loc in locations: topo_mats.append(topo_rotation_matrix(loc.lat.deg, loc.lon.deg)) - pos_mes.append( - [loc.x.to_value("km"), loc.y.to_value("km"), loc.z.to_value("km")] - ) + pos_mes.append([loc.x.to_value("km"), loc.y.to_value("km"), loc.z.to_value("km")]) return np.array(topo_mats), np.array(pos_mes) @@ -174,9 +171,7 @@ def make_transform(coo, toframe): if totopo: # Origin = station position relative to frame origin, in coo's frame if frame_id == 0: # SSB / ICRS - orig_pos_list.append( - station_pos_ssb_j2000(pos_me, et_arr)[0] - ) + orig_pos_list.append(station_pos_ssb_j2000(pos_me, et_arr)[0]) else: # Moon center / MCMF orig_pos_list.append(pos_me) else: From 56c744e1db94a8aacfad1265b161249dff40b94d Mon Sep 17 00:00:00 2001 From: aelanman Date: Tue, 17 Mar 2026 13:41:16 -0400 Subject: [PATCH 08/21] add spice-fixtures data --- lunarsky/tests/data/spice_fixtures.npz | Bin 0 -> 9416 bytes 1 file changed, 0 insertions(+), 0 deletions(-) create mode 100644 lunarsky/tests/data/spice_fixtures.npz diff --git a/lunarsky/tests/data/spice_fixtures.npz b/lunarsky/tests/data/spice_fixtures.npz new file mode 100644 index 0000000000000000000000000000000000000000..1d3e4e8c8fc64d613dd7e06891035e0bc8c66e4a GIT binary patch literal 9416 zcmeI22T&B*)_?~UkRTvA2ye`=U;>PU+l`7SieToUA_xc)B?}Tn7X(MZRRjz;ilPKT zVO3D^cCU&cC^_de5@rMh1VP~MQ5tvV?N+^7|Ejn3|F5TRox0O~&fIhEIo!^#hZPqMSqb{~_c6E2Pa58sxwX$Y^uW#YtW=(wWW^3VMP4tyltWZ-?qxw<*`XOi% z@}RpuFCV|}YP+b#;y3(Idy&4<-~)VkXDhF=nLP?t9t_-f`W8x;U#1t3n8B9oOT@BK z<8#rLPqkmLWs6ns>#u#0TLtLd9K+eXe5I%{M?ge&#Tt_V5hVh0RDJpK<>v0r=8o26 zU}C0!iwr4$G%!>RHLAv6kB;ExQ-+IND>}fPH?(+7bl!Ad=|BouH)1dF!*!ELjD&b}=+Fqw57WjNQq)>HR6$r_@aMtom z7wl`Q7TK}43hw>f-@87R1@adfYFMY11FP`*{DN6IKuXZJB&<3T=GkTC$7a+(e0lKo z=MpRs7=O%g+oK*hq;_cI_`z&o5!B@`ck2Prr_mpF^MY=xTk1Dh*#mHQ(aWMo zAAu_;f~)x~?*J`Np&F)g7JSrA^-d0Bf$#`*2i@FCp#8I8#)C(v;Is=VdVB+w@RI7l z8TAwv@bfgj5O${pbVQkZVOj5C#$#&KnvLbKE&E8^!!Mm6*Cf><-?15l?o*3@`t}Ww zbO?WDb?7oMDtwwIC0z}7(&W?p9a$hwwM7WqT?21kaH)9W)&)-uCD=uuuY$tr&AHot zVga#AAFD@2mGHQZRKX|T9O#h8m-zGM!{E|IyH!auYT@02_BQot#C(}icyH@{fV*X1 zT=G0r0xz2MN?oG8fZ>N^O^;=l0VngD+l|}0p+SE5hx?(GP%iGXC2h%7*x2^@5|3jw z6zvYuJMvE!@OpF>TQj!<@+>qsZhgEJynMg5Q;=Q=HWy@Kp)=YcC_@#_*)ibl)e;TW z&BbtE+)L)zHCv#B%HMm~PzjyPiuvm9v%m|l@@H%ObAfA6vUGVq6D)s_l1G)Q2FJ~K zZnRkUz~?6@gL;auK&@<3Vfwu+u=>^gQi|4PQ0D0``nIAPp1a7j%399?^L*-MM7J^l zpP>HeBJn!7&x>W^aPlQQpcL~&w5lGMO}{GTGu#OtW{bSd53T?PsN_w>{uo$gUwzIY zqY7?{^-#JM!~&d!om~`gk1T9SqOcvEot@0BEnMAg&A-+-S>2J-zEw9Vf4I6={q@o& z?v{da=1*T)mP+6Sj`Y3CHN5zy$T54NlOx!KyzPNlgI1tux?+K(%X@fC@#@9XE|u{7 zsoW*rv-$9q0hJSos{|TkrgB1sGLDFjY@P8qmIW8f(KFSHY9LHfcGI4j35`lvtGivL z(aov%NAK(%!rqi!3(!4Sf;BF@YN_Xy2TZJ81TWoZg7$%{m%ilJgJWi+O7bZ@c*llU z??)1(@ZBq~w3&pCV>@75P(>{hJRG_C@$&gH7-MIBB~befJk-fNU8E?B^I)|%e_0lV zNwQuKZcoiaqJ6DHx^CYFPX+2H81>!oSkZN=$7nrV8kCy4%tZi?3JX-tT{s(+{2&l? zRd)nYuz?qaM46C!V*go2Llbt_7Od!E>Mm_=2lIX8lZ!6}#8%18- z$biD1Kfuv_+1GMQD*>%4DZFWZHW)C&q?m{dnvYOMnyo0POli3Jk~^hH->k5rt%u8C zQ}~ghY4%-Ums;4sG>K~9H6VEHrKA{I7o0Z8>o0<@_7Z-eB22-9iW$Y})eI1zkeDuE z*anukxmpxORe*V3A#o-1K42c1ar2VPBye86d$szF6nyU@xyDh=J}^USzUIV<0+<(b zT=uFc4Zy}zY8Byq(KJ1Yw`*yH2DNsgyZ)A4C$>i#ClG@SL>M_d%sK~63% zU6|A>gguZ-^#3$I4EGETxofPo#m7uif5{G{mlj6HC+vSdXq!sIk-2IgT=(lB#ziNVjw?(+&Dd4v z2524de1mJJK(hlbU+z*e?=BsGMXPXqP)ozh-rUjUk<>+!Xv~tdbbf5k+D@0H(Sva5 zn+2O7ZjZa#UQw$JrsI4yd{1bJG~934;z5#z4)XGeO|p-G z8GkwL^}2R1IzGGP#cJJ18Xg{@8TMq8HqvT!dra$q5N6cKOJ#2U4C_@FQYeod@x6{c z<0pOTxJGkPz>6z1T>sO2*4$}K5qD`gJ$Jfa zRkAkHHY+SX{EHA4^rk~0@uy*~oR(Jn#|bj0*-J2XI)A-ee_w);H)m4*r*G1xJigvZ zRC45?sb#bwMKm2ZK?{Yb4QpiBGtT)^MU4*tZmnB??Q7U?D1=I zH*xIIMXi-UA8kB58rfs-g{``0&_`dPf?bKF#K>fH!z0&<+;{5VK$*%p*~gEDpuJi{ zMx|mX9^04@+GTtYiFi>*?b+^$O-peaNRDm;0b$9wCeI19b#Pp<0^UF)(@#XqG}?-e z(`5p5Mm!M#tux1+%1yAy=%fhD@Ex$IidZSYxQv?|IlxyX9fM-EI@kaCNE6LUUfpe_ z>Wp*;s9osvI)+&FUI@w#?1txddSrObyMY(2Z%Tij7mOcr{kWz_&Ik`x#FeEY4q#(Y zZn%G)C$izS(NE>M#W1LrCqba;99}AJK6ucKj#GZwa7F`2;z}wiO{FRxnB>v;m}2L{ z$h3^?>kZ6iFh9F&*?Qq9ly-A5B3Ez?H{MfzSug#ZJ=7fP8R@NPT&0&E(=CJf} z(dBLBu;N)^yt4RBv@Mx;*GgYIG>}TW_P%`{PPba((H)K;g|1ItMeXpy_zlgplBGUE z)G<`&{Qa|dPIubQ@W*ttNhQt2tYkGBequan`ArYxV^(vYzPKe4FfP65{Gw*ArANl6 zb-$4bEleczw^hgYY(B&DZOK8(AI|0<*iEpP9HI)r)X^O>Xm5P6R$xRsw*0;Dhd4?L z(jgptdFC-0bOAL}3%}BWAiT3`z3Yn+Cy8fb4dc?Nzrklix04KPW9p{`NMkjoEf}fO z?j?)XTeK)j&8|d#eUiB9m)A{L8AW)_oM9QXY)($^(PwQKv$OnNzG4w}vUNgDe7iJi zdvbH3esm^Q)){m+O1B7ktds0e%aBHI7z|wt*-?kk`Rr4)zqB9;&RNUrLS)gqHm~Ot zEQmlPF*MphnSms0D9qtMEQ8)oFWN0s--3OK)ToNxUyEG-bTMXUqZBHCOWg?@e1aTG z9vpD*zl$ZEA8y^hNfr%z7uM2L)r9%P&#_5g{vP9e)ZSPgc9NV$Sz;F5on4&G-P|qQ z?TCe>1JSIkkC2PVGn0E~QvT44{?PjI>xNr!ayS%Dk1``ETmIZ53-`aXM_IH=ESiv zC22{Br(&6mClqc`OUy!xFuL1r+o@L<&-tPTa3H_GkPoo zA4fM%VW>~F3B#L}y8zF+Gmc5D;wj9eVxFvE<5C1iw2BQ8n1;ko4R7o9_ulS= ziFXcIs_bG;VJM_tU}@Zi8t^hlrpE3`2_P{O{erqv-Ki27KkPkLy{%*lJ4M&LmV_u1 zhMQ>Hl@R!h#89(~hkT0@DnT)R#oH*S4Um}WaiOETJTC+0>fcJ8(VGED?4&KXFkk6& z1C)O{!nbr)2Ou$2LTP2xo2|{@%}-Tu|J>#&%w*_4`!H?13Fa#3xhAhw?uuz=dClq`>kwkgcS=yd#TnMW(AnlY2&t>l?9qzT?rQL5T$kG=1U1@(`{gL%d%Kvm-ztsw``TZ@f z;9Gy{md1693!J=u>O6VEW?XDBTs0Nmpna$_M-2bmkeFsbhr!ge@?SbP~7hJ*7;Q%Bd> z;xtd=c=t}4!s*+m<OgIz}ILOo_N zAzyLg$8`jN3*T~mSsel3lt*f;{RjY;AYyyIWQ8uW)T8t5v#kQybDvLQlobPTePe8+ z0|DS{ezmzyB>>!e9ln}l1b|bpJbP3-MH|Vf5|pM631M>j2S@rBasxPlfskr4x7j;g zj@qz-$+d%R>x8%0R-RHjk z|C@bdL!N0%rgOgVz<#5F-S~eOoWQRPZu0d9ZYB0h0qo@&yNPlr%_h{iUJ>9{X8-ZX z=G5#a$)VgegYt)eLvmYS|H(sI5a+Of5iWc!NRB5Cw>}%Fllq*)zetcFF1*gI%!Ze2 zFv)Jh96xS7i^732xi#6qj13^zO_oDbl^DgJ;WKV4Y)JLTR>uCt%6DYNZG{bnVu&%Y zn+(TT3dOj_@*NX#8)3s2awrlUM()q%GV%?@a2sF)k6c2L-8d)0kC<03{K2ivh6H31 z5LOFh6iOhE1rZQZ{oA<bI4qArnZRLwBtLkwB9S q^k0`g+;-Rx7Sax9Mo1}f+4 Date: Tue, 17 Mar 2026 16:45:00 -0400 Subject: [PATCH 09/21] comment on new test setup --- lunarsky/tests/test_spice.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/lunarsky/tests/test_spice.py b/lunarsky/tests/test_spice.py index 6bf0104..13d7824 100644 --- a/lunarsky/tests/test_spice.py +++ b/lunarsky/tests/test_spice.py @@ -4,6 +4,8 @@ import lunarsky.spice_utils as spice_utils +# spice_fixtures.npz contains saved results of the earlier version of lunarsky that was based on the SPICE toolkit +# Tests here are meant to verify that the new code is consistent with older results. FIXTURE_PATH = os.path.join(os.path.dirname(__file__), "data", "spice_fixtures.npz") From 7985baa044953a9af22010fb429b23e051660691 Mon Sep 17 00:00:00 2001 From: aelanman Date: Wed, 18 Mar 2026 11:10:35 -0400 Subject: [PATCH 10/21] update README and comment --- CHANGELOG.md | 9 +++++++++ README.md | 10 +++++----- lunarsky/spice_utils.py | 3 ++- lunarsky/tests/test_spice.py | 5 +++-- 4 files changed, 19 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c88945a..e7a3faf 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/README.md b/README.md index ee96430..6a72717 100644 --- a/README.md +++ b/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 @@ -53,7 +51,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 +61,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 diff --git a/lunarsky/spice_utils.py b/lunarsky/spice_utils.py index 45e749d..efeb0ad 100644 --- a/lunarsky/spice_utils.py +++ b/lunarsky/spice_utils.py @@ -98,7 +98,8 @@ def _R3(a): return np.array([[c, s, 0], [-s, c, 0], [0, 0, 1]]) -# PA to ME constant rotation from moon_080317.tf +# PA to ME constant rotation from moon_080317.tf (FRAME_31007 = MOON_ME_DE421). +# These angles are specific to DE421; they will differ for other ephemerides. # TKFRAME_31007_ANGLES = (67.92, 78.56, 0.30) arcseconds, AXES = (3, 2, 1) # SPICE convention: R_PA_to_ME = R1(-a3) @ R2(-a2) @ R3(-a1) _a1 = np.radians(67.92 / 3600) diff --git a/lunarsky/tests/test_spice.py b/lunarsky/tests/test_spice.py index 13d7824..09c401b 100644 --- a/lunarsky/tests/test_spice.py +++ b/lunarsky/tests/test_spice.py @@ -4,8 +4,9 @@ import lunarsky.spice_utils as spice_utils -# spice_fixtures.npz contains saved results of the earlier version of lunarsky that was based on the SPICE toolkit -# Tests here are meant to verify that the new code is consistent with older results. +# spice_fixtures.npz contains saved results of the earlier version of lunarsky that was based +# on the SPICE toolkit. Tests here are meant to verify that the new code is consistent with +# older results. FIXTURE_PATH = os.path.join(os.path.dirname(__file__), "data", "spice_fixtures.npz") From f3d5952824af6a9220e0d854f7bc3249aec006e9 Mon Sep 17 00:00:00 2001 From: aelanman Date: Thu, 19 Mar 2026 11:14:33 -0400 Subject: [PATCH 11/21] ensure TDB used consistently --- lunarsky/mcmf.py | 2 +- lunarsky/spice_utils.py | 4 ++-- lunarsky/time.py | 2 +- lunarsky/topo.py | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/lunarsky/mcmf.py b/lunarsky/mcmf.py index f497387..9371b10 100644 --- a/lunarsky/mcmf.py +++ b/lunarsky/mcmf.py @@ -66,7 +66,7 @@ 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 diff --git a/lunarsky/spice_utils.py b/lunarsky/spice_utils.py index efeb0ad..cc2aeec 100644 --- a/lunarsky/spice_utils.py +++ b/lunarsky/spice_utils.py @@ -205,7 +205,7 @@ def _pos_ssb_j2000(body_id, ets): jd = _et_to_jd(ets) if body_id == 301: - emb = np.asarray(spk[0, 3].compute(jd)).T + emb = np.asarray(spk[0, 3].compute(jd)).T # spk keys are tuples (center, target) moon_emb = np.asarray(spk[3, 301].compute(jd)).T return emb + moon_emb elif body_id == 399: @@ -299,7 +299,7 @@ def earth_pos_mcmf(obstimes): Used for tests. """ - ets = (obstimes - Time("J2000")).sec + ets = (obstimes.tdb - Time("J2000")).sec earthpos = body_position(399, ets, "MOON_ME", 301) earthpos = unit.Quantity(earthpos.T, "km") return MCMF(*earthpos, obstime=obstimes) diff --git a/lunarsky/time.py b/lunarsky/time.py index 2461a14..23b46a4 100644 --- a/lunarsky/time.py +++ b/lunarsky/time.py @@ -70,7 +70,7 @@ def sidereal_time(self, kind, longitude=None, model=None): # to the selenodetic coordinate system. "self.location" must # be defined in order to get here. - et = np.atleast_1d((self - Time("J2000")).sec) + et = np.atleast_1d((self.tdb - Time("J2000")).sec) mats = moon_me_to_j2000(et) # Zenith vector diff --git a/lunarsky/topo.py b/lunarsky/topo.py index 27e9a7c..2fa473f 100644 --- a/lunarsky/topo.py +++ b/lunarsky/topo.py @@ -129,7 +129,7 @@ def make_transform(coo, toframe): raise ValueError("location must be defined for LunarTopo transformations") # Make arrays - ets = (obstime - _J2000).sec + ets = (obstime.tdb - _J2000).sec topo_mats, pos_mes = _get_topo_data(location) loc_indices = np.arange(len(topo_mats)) shape_out = np.broadcast_shapes(coo.shape, ets.shape, location.shape) From bdb823e76743a69694617f2e1f3c2f338783fe87 Mon Sep 17 00:00:00 2001 From: aelanman Date: Thu, 19 Mar 2026 11:15:58 -0400 Subject: [PATCH 12/21] linting --- lunarsky/spice_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lunarsky/spice_utils.py b/lunarsky/spice_utils.py index cc2aeec..47b3fde 100644 --- a/lunarsky/spice_utils.py +++ b/lunarsky/spice_utils.py @@ -205,7 +205,7 @@ def _pos_ssb_j2000(body_id, ets): jd = _et_to_jd(ets) if body_id == 301: - emb = np.asarray(spk[0, 3].compute(jd)).T # spk keys are tuples (center, target) + emb = np.asarray(spk[0, 3].compute(jd)).T # spk keys are tuples (center, target) moon_emb = np.asarray(spk[3, 301].compute(jd)).T return emb + moon_emb elif body_id == 399: From ce28ecb6bd80afadc4fe1848e5af95ac1f8dd026 Mon Sep 17 00:00:00 2001 From: aelanman Date: Mon, 13 Apr 2026 15:11:38 -0400 Subject: [PATCH 13/21] linting --- lunarsky/spice_utils.py | 1 - 1 file changed, 1 deletion(-) diff --git a/lunarsky/spice_utils.py b/lunarsky/spice_utils.py index 16519dd..6ea3302 100644 --- a/lunarsky/spice_utils.py +++ b/lunarsky/spice_utils.py @@ -1,7 +1,6 @@ import numpy as np import os import urllib.request -import warnings from filelock import FileLock from astropy.utils.console import ProgressBar from astropy.coordinates.matrix_utilities import rotation_matrix From 273ca584b8d37ae01bedbc081fb6560a838833d9 Mon Sep 17 00:00:00 2001 From: aelanman Date: Mon, 11 May 2026 17:05:10 -0400 Subject: [PATCH 14/21] clean up namespace --- lunarsky/moon.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/lunarsky/moon.py b/lunarsky/moon.py index a2f3b44..5843281 100644 --- a/lunarsky/moon.py +++ b/lunarsky/moon.py @@ -12,18 +12,14 @@ 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 @@ -36,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. From 7a7987d5876a6d86476ad2a99fc05f2f50be5de2 Mon Sep 17 00:00:00 2001 From: Adam Lanman Date: Wed, 13 May 2026 10:48:15 -0400 Subject: [PATCH 15/21] Update lunarsky/moon.py Co-authored-by: Leo Singer --- lunarsky/moon.py | 1 - 1 file changed, 1 deletion(-) diff --git a/lunarsky/moon.py b/lunarsky/moon.py index 5843281..18144af 100644 --- a/lunarsky/moon.py +++ b/lunarsky/moon.py @@ -32,7 +32,6 @@ class GSFCSelenodeticRepresentation(BaseGeodeticRepresentation): _equatorial_radius = 1738.1e3 * u.m _flattening = 0.0012 -e class GRAIL23SelenodeticRepresentation(BaseGeodeticRepresentation): """Lunar ellipsoid defined by gravimetry of GRAIL data. From 74e34bd0e1a5c58e72d1f188f410a271d02bc3dc Mon Sep 17 00:00:00 2001 From: aelanman Date: Wed, 13 May 2026 10:57:44 -0400 Subject: [PATCH 16/21] fix lost dependencies --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index eaf220e..cc26451 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,8 +24,8 @@ keywords = [ dependencies = [ "numpy>=1.15", "astropy>=6.0.0", - "spiceypy", "jplephem", + "filelock", "fastkml[lxml]", ] requires-python = ">=3.8" From 574e5c136d5b1adce789ad5eb0876bb658cbbba0 Mon Sep 17 00:00:00 2001 From: aelanman Date: Wed, 13 May 2026 10:59:12 -0400 Subject: [PATCH 17/21] linting --- lunarsky/moon.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lunarsky/moon.py b/lunarsky/moon.py index abaf8cc..85e0b40 100644 --- a/lunarsky/moon.py +++ b/lunarsky/moon.py @@ -18,8 +18,7 @@ class SPHERESelenodeticRepresentation(BaseGeodeticRepresentation): - """Lunar ellipsoid as a sphere - """ + """Lunar ellipsoid as a sphere""" _equatorial_radius = 1737.1e3 * u.m _flattening = 0.0 @@ -34,6 +33,7 @@ class GSFCSelenodeticRepresentation(BaseGeodeticRepresentation): _equatorial_radius = 1738.1e3 * u.m _flattening = 0.0012 + class GRAIL23SelenodeticRepresentation(BaseGeodeticRepresentation): """Lunar ellipsoid defined by gravimetry of GRAIL data. From e7c476094d80447666eb7c4b2ebdac40a9f34caa Mon Sep 17 00:00:00 2001 From: aelanman Date: Wed, 13 May 2026 13:24:47 -0400 Subject: [PATCH 18/21] use astropy for downloading/caching kernels --- .github/workflows/testsuite.yaml | 8 ++--- .gitignore | 4 --- README.md | 6 ++-- lunarsky/spice_utils.py | 52 +++++++++----------------------- lunarsky/tests/test_spice.py | 15 ++++----- pyproject.toml | 1 - 6 files changed, 28 insertions(+), 58 deletions(-) diff --git a/.github/workflows/testsuite.yaml b/.github/workflows/testsuite.yaml index 71faf31..fb28529 100644 --- a/.github/workflows/testsuite.yaml +++ b/.github/workflows/testsuite.yaml @@ -45,13 +45,13 @@ jobs: id: cache-kernels uses: actions/cache@v3 with: - path: lunarsky/data/spk/planets/de430.bsp + path: ~/.astropy/cache/download 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()" + pip install . + python -c "from lunarsky.spice_utils import download_big_kernels; download_big_kernels(show_progress=False)" tests: env: @@ -77,7 +77,7 @@ jobs: - name: load-cache uses: actions/cache@v3 with: - path: lunarsky/data/spk/planets/de430.bsp + path: ~/.astropy/cache/download key: kernel-de430-${{ runner.os }} - name: Install run: | diff --git a/.gitignore b/.gitignore index ff06c0d..77e492f 100644 --- a/.gitignore +++ b/.gitignore @@ -139,7 +139,3 @@ 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 diff --git a/README.md b/README.md index 56a0105..f7ef392 100644 --- a/README.md +++ b/README.md @@ -31,13 +31,13 @@ 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: +The first time lunarsky is used, it will automatically download the JPL DE430 ephemeris +(~150 MB) into astropy's data cache (`~/.astropy/cache/`). Subsequent uses read from the +cache. To pre-fetch the kernel (e.g. before going offline, or in CI), 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 diff --git a/lunarsky/spice_utils.py b/lunarsky/spice_utils.py index 6ea3302..fc79112 100644 --- a/lunarsky/spice_utils.py +++ b/lunarsky/spice_utils.py @@ -1,8 +1,6 @@ import numpy as np import os -import urllib.request -from filelock import FileLock -from astropy.utils.console import ProgressBar +from astropy.utils.data import download_file from astropy.coordinates.matrix_utilities import rotation_matrix from astropy.time import Time import astropy.units as unit @@ -168,36 +166,18 @@ def moon_me_to_j2000(ets): _NAIF_KERNEL_URL = "https://naif.jpl.nasa.gov/pub/naif/generic_kernels/" -def download_big_kernels(): +def download_big_kernels(show_progress=True): """ - Download large (~150 MB) SPICE kernel files that can't be - included in the python package. + Pre-fetch large (~150 MB) SPICE kernel files into astropy's data cache. - Installs them to the package data directory. Must be run once - before using lunarsky transformations that depend on these kernels. + Calling this is optional — the kernels are auto-downloaded on first use. + It is useful for offline preparation or to control when the download + happens (e.g. in CI, or before going offline). """ - paths = [] - for kern in _BIG_KERNELS: - kf = os.path.join(DATA_PATH, kern) - paths.append(kf) - if os.path.exists(kf): - continue - kurl = _NAIF_KERNEL_URL + kern - os.makedirs(os.path.dirname(kf), exist_ok=True) - with urllib.request.urlopen(kurl) as response: - total_size = int(response.headers.get("Content-Length", 0)) - chunk_size = max(4096, total_size // 20) if total_size > 0 else 4096 - downloaded = 0 - with ProgressBar( - total_size // chunk_size + 1, ipython_widget=False - ) as pbar, open(kf, "wb") as ofile, FileLock(kf + ".lock"): - while cur_buf := response.read(chunk_size): - downloaded += len(cur_buf) - if total_size > 0: - pbar.update(downloaded // chunk_size) - ofile.write(cur_buf) - - return paths + return [ + download_file(_NAIF_KERNEL_URL + kern, cache=True, show_progress=show_progress) + for kern in _BIG_KERNELS + ] def _ensure_spk(): @@ -205,13 +185,11 @@ def _ensure_spk(): if _SPK is None: from jplephem.spk import SPK - spk_path = os.path.join(DATA_PATH, "spk/planets/de430.bsp") - if not os.path.exists(spk_path): - raise FileNotFoundError( - f"Required kernel file not found: {spk_path}. " - "Please run lunarsky.spice_utils.download_big_kernels() " - "to download it." - ) + spk_path = download_file( + _NAIF_KERNEL_URL + "spk/planets/de430.bsp", + cache=True, + show_progress=True, + ) _SPK = SPK.open(spk_path) return _SPK diff --git a/lunarsky/tests/test_spice.py b/lunarsky/tests/test_spice.py index fc5d90e..2c90c5d 100644 --- a/lunarsky/tests/test_spice.py +++ b/lunarsky/tests/test_spice.py @@ -2,7 +2,6 @@ import os import pytest import lunarsky.spice_utils as spice_utils -from lunarsky.data import DATA_PATH # spice_fixtures.npz contains saved results of the earlier version of lunarsky that was based @@ -18,16 +17,14 @@ def fixtures(): def test_kernels_available(): """ - Ensure the large SPK kernel is available, downloading it if not. + Ensure the large SPK kernel is downloadable and cached. - This is a no-op when the kernel is already present (e.g. via CI cache or - a prior download), so it won't hit the JPL server unnecessarily. It does - double duty as a smoke test of download_big_kernels() on fresh environments. + On a fresh environment this hits the JPL server; on subsequent runs + (and in CI with a populated astropy cache) it is a no-op. """ - spk_path = os.path.join(DATA_PATH, "spk", "planets", "de430.bsp") - if not os.path.exists(spk_path): - spice_utils.download_big_kernels() - assert os.path.exists(spk_path) + paths = spice_utils.download_big_kernels(show_progress=False) + for p in paths: + assert os.path.exists(p) def test_j2000_to_moon_me(fixtures): diff --git a/pyproject.toml b/pyproject.toml index cc26451..57c857c 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,7 +25,6 @@ dependencies = [ "numpy>=1.15", "astropy>=6.0.0", "jplephem", - "filelock", "fastkml[lxml]", ] requires-python = ">=3.8" From 17f428626d379afba04d25f68067dde2200b8b1a Mon Sep 17 00:00:00 2001 From: aelanman Date: Wed, 13 May 2026 14:30:28 -0400 Subject: [PATCH 19/21] remove unnecessary macos patch --- lunarsky/__init__.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/lunarsky/__init__.py b/lunarsky/__init__.py index 35c308e..9eeef94 100644 --- a/lunarsky/__init__.py +++ b/lunarsky/__init__.py @@ -1,11 +1,3 @@ -import sys - -if sys.platform == "darwin": - # On MacOS, need to set multiprocessing to fork by default - import multiprocessing - - multiprocessing.set_start_method("fork", force=True) - from .moon import * # noqa from .mcmf import * # noqa from .topo import * # noqa From 5d3670440e669fb75c91f46962e31a49df2d4615 Mon Sep 17 00:00:00 2001 From: aelanman Date: Fri, 15 May 2026 09:49:20 -0400 Subject: [PATCH 20/21] fix typo --- lunarsky/moon.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lunarsky/moon.py b/lunarsky/moon.py index 46460c1..c9166e7 100644 --- a/lunarsky/moon.py +++ b/lunarsky/moon.py @@ -527,8 +527,8 @@ def get_mcmf(self, obstime=None): mcmf = property( get_mcmf, - doc="""An `~astropy.coordinates.MCMF` object with - or the location of this object at the + doc="""An `~astropy.coordinates.MCMF` object + for the location of this object at the default ``obstime``.""", ) From 5fcd27ae0f5c32793c32730d8fc15deda12a721a Mon Sep 17 00:00:00 2001 From: aelanman Date: Thu, 28 May 2026 15:12:14 -0400 Subject: [PATCH 21/21] remove cruft; properly handle file open/close --- lunarsky/moon.py | 50 ----------------------------------------- lunarsky/spice_utils.py | 23 ++++++++++--------- 2 files changed, 12 insertions(+), 61 deletions(-) diff --git a/lunarsky/moon.py b/lunarsky/moon.py index c9166e7..fa7ffe8 100644 --- a/lunarsky/moon.py +++ b/lunarsky/moon.py @@ -267,56 +267,6 @@ def of_site(cls, site_name: str): lon_lat = sites[site_name.lower()] return cls.from_selenodetic(*lon_lat, ellipsoid="IAU2000") - @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): """ diff --git a/lunarsky/spice_utils.py b/lunarsky/spice_utils.py index fc79112..d312ab4 100644 --- a/lunarsky/spice_utils.py +++ b/lunarsky/spice_utils.py @@ -25,21 +25,22 @@ def _read_bpc(filepath): """ from jplephem.daf import DAF - daf = DAF(open(filepath, "rb")) + with open(filepath, "rb") as f: + daf = DAF(f) - segments = list(daf.summaries()) - if len(segments) != 1: - raise ValueError(f"Expected 1 segment in binary PCK, got {len(segments)}") + segments = list(daf.summaries()) + if len(segments) != 1: + raise ValueError(f"Expected 1 segment in binary PCK, got {len(segments)}") - name, values = segments[0] - data_type = int(values[4]) - start_addr = int(values[5]) - end_addr = int(values[6]) + name, values = segments[0] + data_type = int(values[4]) + start_addr = int(values[5]) + end_addr = int(values[6]) - if data_type != 2: - raise ValueError(f"Unsupported binary PCK data type {data_type}") + if data_type != 2: + raise ValueError(f"Unsupported binary PCK data type {data_type}") - arr = daf.read_array(start_addr, end_addr) + arr = daf.read_array(start_addr, end_addr) init_epoch = arr[-4] interval = arr[-3]