Merge pull request #19 from aelanman/dimensionless_bug

addressing behavior in transforms with dimensionless positions

.github/workflows/testsuite.yaml

@@ -62,7 +62,7 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest]
-        python-version: [3.6, 3.7, 3.8]
+        python-version: [3.7, 3.8, 3.9]
     steps:
       - uses: actions/checkout@v2
         with:

.pre-commit-config.yaml

@@ -1,8 +1,8 @@
-exclude: '(^lunarsky/data/*)'
+exclude: '(^lunarsky/data/)'
 
 repos:
-    - repo: git://github.com/pre-commit/pre-commit-hooks
-      rev: v2.5.0
+    - repo: https://github.com/pre-commit/pre-commit-hooks
+      rev: v4.3.0
       hooks:
           - id: trailing-whitespace
           - id: check-ast
@@ -23,7 +23,7 @@ repos:
             - flake8-comprehensions
             - flake8-pytest
     - repo: https://github.com/psf/black
-      rev: 20.8b1
+      rev: 22.3.0
       hooks:
           - id: black
             language_version: python3.8

CHANGELOG.md

@@ -3,8 +3,15 @@
 ## Unreleased
 
 ## Fixed
+- Updated version of pre-commit-hooks used
+- Accept tuple for location in Time class (in this case, assumes EarthLocation)
+- Use newer PCK file in unit test with Earth positioning.
+- Match behavior of astropy when transforming non-unit cartesian positions without units (treat as unitspherical using direction info only)
 - Now tracking available lunar station_ids, instead of incrementing a counter naively.
 
+## Deprecated
+- Dropping support for Python 3.6
+
 ## [0.1.2] -- 2022-01-17
 
 ## Added

lunarsky/mcmf.py

@@ -94,6 +94,9 @@ def make_transform(coo, toframe):
 
     newrepr = coo_cart.transform(mats).reshape(shape_out)
 
+    if is_unitspherical:
+        newrepr = newrepr.represent_as(UnitSphericalRepresentation)
+
     return toframe.realize_frame(newrepr)
 
 

lunarsky/moon.py

@@ -10,7 +10,7 @@
 )
 from astropy.coordinates.attributes import Attribute
 
-from .spice_utils import remove_topo
+from .spice_utils import remove_topo, LUNAR_RADIUS
 
 __all__ = ["MoonLocation", "MoonLocationAttribute"]
 
@@ -113,7 +113,7 @@ class MoonLocation(u.Quantity):
     _location_dtype = np.dtype({"names": ["x", "y", "z"], "formats": [np.float64] * 3})
     _array_dtype = np.dtype((np.float64, (3,)))
 
-    _lunar_radius = 1737.1e3  # m
+    _lunar_radius = LUNAR_RADIUS
 
     # Manage the set of defined ephemerides.
     # Class attributes only
@@ -282,7 +282,7 @@ def from_selenodetic(cls, lon, lat, height=0.0):
         the mean position of the "sub-Earth" point on the lunar surface.
 
         """
-        lon = Longitude(lon, u.degree, wrap_angle=180 * u.degree, copy=False)
+        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.
         if not isinstance(height, u.Quantity):

lunarsky/spice_utils.py

@@ -12,9 +12,13 @@
 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):
     """
@@ -53,7 +57,6 @@ def furnish_kernels():
         "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)
@@ -94,9 +97,10 @@ def lunar_surface_ephem(latitude, longitude, station_num=98):
 
     lat = np.radians(latitude)
     lon = np.radians(longitude)
-    lunar_radius = 1737.1  # km
     ets = np.array([spice.str2et("1950-01-01"), spice.str2et("2150-01-01")])
-    pos_mcmf = spice.latrec(lunar_radius, lon, lat)  # TODO Use MoonLocation instead?
+    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)
@@ -207,7 +211,7 @@ def remove_topo(station_num):
 
     frame_vars = [s.format(idnum, station_name, fm_center_id) for s in fmt_vars]
 
-    # Handle a glitch in spiceypy for older versions of numpy
+    # Handle a bug in spiceypy for older versions of numpy
     if np.str_ is None:
         return
     for var in frame_vars:

lunarsky/tests/test_time.py

@@ -1,5 +1,5 @@
 import numpy as np
-from astropy.coordinates import ICRS
+from astropy.coordinates import ICRS, EarthLocation
 from astropy.time import TimeDelta
 import pytest
 
@@ -23,3 +23,15 @@ def test_sidereal_time_calculation(lat, lon):
         src = SkyCoord(alt="90d", az="0d", frame="lunartopo", obstime=tt, location=loc)
         lst = tt.sidereal_time("mean")
         assert np.isclose(lst.deg, src.transform_to(ICRS()).ra.deg, atol=1e-4)
+
+
+def test_earthloc():
+    lat = -30.0
+    lon = 127.0
+
+    t0 = Time(
+        "2020-01-01T00:00:00", location=EarthLocation.from_geodetic(lon=lon, lat=lat)
+    )
+    t1 = Time("2020-01-01T00:00:00", location=(lon, lat))
+
+    assert t0.location == t1.location

lunarsky/tests/test_transforms.py

@@ -3,7 +3,7 @@
 from lunarsky.time import Time, TimeDelta
 import astropy.coordinates as ac
 from astropy import units as un
-from astropy.utils import IncompatibleShapeError
+from astropy.utils import IncompatibleShapeError, exceptions
 from astropy.tests.helper import assert_quantity_allclose
 import lunarsky
 import pytest
@@ -25,6 +25,7 @@
 
 @pytest.mark.parametrize("time", jd_10yr)
 @pytest.mark.parametrize("lat,lon", latlons_grid)
+@pytest.mark.filterwarnings("ignore::erfa.ErfaWarning")
 def test_icrs_to_topo_long_time(time, lat, lon, grcat):
     # Check that the following transformation paths are equivalent:
     #   ICRS -> MCMF -> TOPO
@@ -48,9 +49,9 @@ def test_icrs_to_topo_long_time(time, lat, lon, grcat):
     "time, lat, lon",
     [
         (t0, 11.2, 1.4),
-        (t0, (10.3, 11.2), (0.0, 1.4)),
+        (t0, [10.3, 11.2], [0.0, 1.4]),
         (jd_4mo[:2], 10.3, 0.0),
-        (jd_4mo[:2], (10.3, 11.2), (0.0, 1.4)),
+        (jd_4mo[:2], [10.3, 11.2], [0.0, 1.4]),
     ],
 )
 def test_transform_loops(obj, path, time, lat, lon):
@@ -84,6 +85,7 @@ def test_topo_to_topo():
     assert new.az.deg == 90
 
 
+@pytest.mark.filterwarnings("ignore::erfa.ErfaWarning")
 def test_mcmf_to_mcmf():
     # Transform MCMF positions to MCMF frame half a lunar sidereal day later.
     # Assert that the new positions are roughly close to 180 deg from the original.
@@ -264,3 +266,24 @@ def test_incompatible_transform(fromframe):
     src = lunarsky.SkyCoord(coo)
     with pytest.raises(IncompatibleShapeError):
         src.transform_to(ltop)
+
+
+def test_finite_vs_spherical():
+    # Transform MCMF coordinates with distance and without
+    # Assumes infinite distance if no unit given, as astropy does.
+    R0 = 2e3
+    N = 100
+    phi = np.linspace(0, 2 * np.pi, N)
+    xyz = R0 * un.km * np.array([np.cos(phi), np.sin(phi), np.zeros(N)])
+    with_units = lunarsky.SkyCoord(lunarsky.MCMF(*xyz))
+    xyz = R0 * np.array([np.cos(phi), np.sin(phi), np.zeros(N)])
+    sans_units = lunarsky.SkyCoord(lunarsky.MCMF(*xyz))
+
+    loc = lunarsky.MoonLocation.from_selenodetic(lon=180 * un.deg, lat=0, height=0)
+    altaz_with_units = with_units.transform_to(lunarsky.LunarTopo(location=loc))
+    with pytest.warns(exceptions.AstropyUserWarning, match="Coordinates do not "):
+        altaz_sans_units = sans_units.transform_to(lunarsky.LunarTopo(location=loc))
+
+    radius = lunarsky.spice_utils.LUNAR_RADIUS * un.m
+    assert np.all(altaz_with_units.distance < radius + R0 * un.km)
+    assert_quantity_allclose(altaz_sans_units.distance, 1.0)

lunarsky/time.py

@@ -33,7 +33,7 @@ def __init__(
     ):
 
         super_loc = None
-        if isinstance(location, EarthLocation):
+        if isinstance(location, (EarthLocation, tuple)):
             super_loc = location
 
         super().__init__(
@@ -53,7 +53,7 @@ def __init__(
 
     def sidereal_time(self, kind, longitude=None, model=None):
         # Currently returns the zenith RA as the LST.
-        if self.location is None or self.location is EarthLocation:
+        if self.location is None or isinstance(self.location, EarthLocation):
             return super().sidereal_time(kind, longitude=longitude, model=model)
 
         if model is not None:

lunarsky/topo.py

@@ -1,12 +1,13 @@
 import numpy as np
+import warnings
 from astropy.utils.decorators import format_doc
 from astropy.coordinates.representation import (
     SphericalRepresentation,
     SphericalCosLatDifferential,
     UnitSphericalRepresentation,
     CartesianRepresentation,
 )
-from astropy.utils import check_broadcast
+from astropy.utils import check_broadcast, exceptions
 from astropy.coordinates.baseframe import (
     BaseCoordinateFrame,
     base_doc,
@@ -124,6 +125,9 @@ def make_transform(coo, toframe):
         obstime = toframe.obstime
         location = toframe.location
 
+    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)
@@ -180,11 +184,19 @@ def make_transform(coo, toframe):
             )
             * un.km
         )
-
         coo_cart -= CartesianRepresentation((orig_posvel.T)[:3])
 
     newrepr = coo_cart.transform(mats).reshape(shape_out)
 
+    if is_unitspherical:
+        if not np.allclose(newrepr.norm(), 1.0):
+            warnings.warn(
+                "Coordinates do not all have unit magnitude, but will be treated as unit spherical,"
+                "Define coordinates as Quantity or normalize to remove this warning.",
+                exceptions.AstropyUserWarning,
+            )
+        newrepr = newrepr.represent_as(UnitSphericalRepresentation)
+
     return toframe.realize_frame(newrepr)