diff --git a/.github/workflows/testsuite.yaml b/.github/workflows/testsuite.yaml
index 1452365..cc6ffb3 100644
--- a/.github/workflows/testsuite.yaml
+++ b/.github/workflows/testsuite.yaml
@@ -15,10 +15,10 @@ jobs:
       - uses: actions/checkout@v3
         with:
          fetch-depth: 1
-      - name: Set up Python 3.8
+      - name: Set up Python 3.9
         uses: actions/setup-python@v4
         with:
-          python-version: 3.8
+          python-version: 3.9
       - name: Setup Environment
         run: |
           pip install pre-commit
@@ -36,7 +36,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v4
         with:
-          python-version: 3.8
+          python-version: 3.9
       - name: Get current date
         id: date
         run: echo "::set-output name=date::$(date +'%Y-%m-%d')"
@@ -63,7 +63,7 @@ jobs:
       fail-fast: false
       matrix:
         os: [ubuntu-latest, macos-latest]
-        python-version: [3.7, 3.8, 3.9]
+        python-version: ["3.9", "3.10", "3.11", "3.12"]
     steps:
       - uses: actions/checkout@v3
         with:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 8955147..2866e1b 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -26,6 +26,6 @@ repos:
       rev: 22.3.0
       hooks:
           - id: black
-            language_version: python3.8
+            language_version: python3.9
             args:
               - --line-length=90
diff --git a/CHANGELOG.md b/CHANGELOG.md
index c3bf968..75218fd 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,13 @@
 # Changelog
 
+## [0.2.2] -- 2022-03-19
+
+## Added
+- Support for ellipsoid models for selenodetic coordinates (non-spherical)
+
+## Deprecated
+- Removed support for Python < 3.9. Astropy >= 6.0 is required for ellipsoid support
+
 ## [0.2.0] -- 2022-10-12
 
 ## Fixed
diff --git a/lunarsky/moon.py b/lunarsky/moon.py
index d5ee239..f3df176 100644
--- a/lunarsky/moon.py
+++ b/lunarsky/moon.py
@@ -4,17 +4,83 @@
 from astropy.units.quantity import QuantityInfoBase
 from astropy.coordinates.angles import Longitude, Latitude
 from astropy.coordinates.earth import GeodeticLocation
+from astropy.coordinates.representation.geodetic import BaseGeodeticRepresentation
 from astropy.coordinates.representation import (
     CartesianRepresentation,
-    SphericalRepresentation,
 )
 from astropy.coordinates.attributes import Attribute
 
-from .spice_utils import remove_topo, LUNAR_RADIUS
+from .spice_utils import remove_topo
+
+LUNAR_RADIUS = 1737.1e3  # m
 
 __all__ = ["MoonLocation", "MoonLocationAttribute"]
 
 
+class SPHERESelenodeticRepresentation(BaseGeodeticRepresentation):
+    """Lunar ellipsoid as a sphere
+
+    Radius defined by lunarsky.spice_utils.LUNAR_RADIUS
+    """
+
+    _equatorial_radius = LUNAR_RADIUS * u.m
+    _flattening = 0.0
+
+
+class GSFCSelenodeticRepresentation(BaseGeodeticRepresentation):
+    """Lunar ellipsoid from NASA/GSFC "Planetary Fact Sheet"
+
+    https://nssdc.gsfc.nasa.gov/planetary/factsheet/moonfact.html
+    """
+
+    _equatorial_radius = 1738.1e3 * u.m
+    _flattening = 0.0012
+
+
+class GRAIL23SelenodeticRepresentation(BaseGeodeticRepresentation):
+    """Lunar ellipsoid defined by gravimetry of GRAIL data.
+
+    https://doi.org/10.1007/s40328-023-00415-w
+    """
+
+    _equatorial_radius = 1737576.6 * u.m
+    _flattening = 0.000305
+
+
+class CE1LAM10SelenodeticRepresentation(BaseGeodeticRepresentation):
+    """Lunar ellipsoid from Chang'e 1 laser altimetry.
+
+    Rotation ellipsoid = CE-1-LAM-GEO
+
+    https://doi.org/10.1007/s11430-010-4060-6
+    """
+
+    _equatorial_radius = 1737.632 * u.km
+    _flattening = 1 / 973.463
+
+
+# Define reference ellipsoids
+SELENOIDS = {
+    "SPHERE": SPHERESelenodeticRepresentation,
+    "GSFC": GSFCSelenodeticRepresentation,
+    "GRAIL23": GRAIL23SelenodeticRepresentation,
+    "CE-1-LAM-GEO": CE1LAM10SelenodeticRepresentation,
+}
+
+
+class SelenodeticLocation(GeodeticLocation):
+    """Rename GeodeticLocation class for clarity"""
+
+
+def _check_ellipsoid(ellipsoid=None, default="SPHERE"):
+    """Defaulting lunar ellipsoid"""
+    if ellipsoid is None:
+        ellipsoid = default
+    if ellipsoid not in SELENOIDS:
+        raise ValueError(f"Ellipsoid {ellipsoid} not among known ones ({SELENOIDS})")
+    return ellipsoid
+
+
 class MoonLocationInfo(QuantityInfoBase):
     """
     Container for meta information like name, description, format.  This is
@@ -22,10 +88,12 @@ class MoonLocationInfo(QuantityInfoBase):
     be used as a general way to store meta information.
     """
 
-    _represent_as_dict_attrs = ("x", "y", "z")
+    _represent_as_dict_attrs = ("x", "y", "z", "ellipsoid")
 
     def _construct_from_dict(self, map):
+        ellipsoid = map.pop("ellipsoid")
         out = self._parent_cls(**map)
+        out.ellipsoid = ellipsoid
         return out
 
     def new_like(self, cols, length, metadata_conflicts="warn", name=None):
@@ -53,7 +121,7 @@ def new_like(self, cols, length, metadata_conflicts="warn", name=None):
             Empty instance of this class consistent with ``cols``
         """
         # Very similar to QuantityInfo.new_like, but the creation of the
-        # map is different enough that this needs its own rouinte.
+        # map is different enough that this needs its own routine.
         # Get merged info attributes shape, dtype, format, description.
         attrs = self.merge_cols_attributes(
             cols, metadata_conflicts, name, ("meta", "format", "description")
@@ -93,8 +161,10 @@ class MoonLocation(u.Quantity):
     This class uses the ME frame.
 
     Positions may be defined in Cartesian (x, y, z) coordinates with respect to the
-    center of mass of the Moon, or in ``selenodetic'' coordinates (longitude, latitude).
-    In selenodetic coordinates, positions are on the surface exactly.
+    center of mass of the Moon, or in ``selenodetic'' coordinates (longitude, latitude, height).
+
+    Selenodetic coordinates are defined with respect to a reference ellipsoid. The default is a
+    sphere of radius 1731.1e3 km, but other ellipsoids are available. See lunarsky.moon.SELENOIDS.
 
     See:
         "A Standardized Lunar Coordinate System for the Lunar Reconnaissance
@@ -110,11 +180,10 @@ class MoonLocation(u.Quantity):
     property.
     """
 
+    _ellipsoid = "SPHERE"
     _location_dtype = np.dtype({"names": ["x", "y", "z"], "formats": [np.float64] * 3})
     _array_dtype = np.dtype((np.float64, (3,)))
 
-    _lunar_radius = LUNAR_RADIUS
-
     # Manage the set of defined ephemerides.
     # Class attributes only
     _inuse_stat_ids = []
@@ -170,7 +239,9 @@ def _set_site_id(cls, inst):
             raise ValueError("Too many unique MoonLocation objects open at once.")
 
         for llh in llh_arr:
-            lonlatheight = "_".join(["{:.4f}".format(ll) for ll in llh])
+            lonlatheight = "_".join(
+                ["{:.4f}".format(ll) for ll in llh] + [inst._ellipsoid]
+            )
             if lonlatheight not in cls._existing_locs:
                 new_stat_id = cls._avail_stat_ids.pop()
                 cls._existing_locs.append(lonlatheight)
@@ -248,7 +319,7 @@ def from_selenocentric(cls, x, y, z, unit=None):
         return inst
 
     @classmethod
-    def from_selenodetic(cls, lon, lat, height=0.0):
+    def from_selenodetic(cls, lon, lat, height=0.0, ellipsoid=None):
         """
         Location on the Moon, from latitude and longitude.
 
@@ -264,6 +335,10 @@ def from_selenodetic(cls, lon, lat, height=0.0):
             Height above reference sphere (if float, in meters; default: 0).
             The reference sphere is a sphere of radius 1737.1 kilometers,
             from the center of mass of the Moon.
+        ellipsoid : str, optional
+            Name of the reference ellipsoid to use (default: 'SPHERE').
+            Available ellipsoids are:  'SPHERE', 'GRAIL23', 'CE-1-LAM-GEO'.
+            See docstrings for classes in ELLIPSOIDS dictionary for references.
 
         Raises
         ------
@@ -281,7 +356,11 @@ def from_selenodetic(cls, lon, lat, height=0.0):
         relative to a prime meridian, which is itself given by
         the mean position of the "sub-Earth" point on the lunar surface.
 
+        For the conversion to selenocentric coordinates, the ERFA routine
+        ``gd2gce`` is used.  See https://github.com/liberfa/erfa
+
         """
+        ellipsoid = _check_ellipsoid(ellipsoid, default=cls._ellipsoid)
         lon = Longitude(lon, u.degree, copy=False).wrap_at(180 * u.degree)
         lat = Latitude(lat, u.degree, copy=False)
         # don't convert to m by default, so we can use the height unit below.
@@ -294,23 +373,14 @@ def from_selenodetic(cls, lon, lat, height=0.0):
                     str(lon.shape), str(lat.shape)
                 )
             )
-        # get selenocentric coordinates. Have to give one-dimensional array.
-
-        lunar_radius = u.Quantity(cls._lunar_radius, u.m, copy=False)
-
-        Npts = lon.size
-        xyz = np.zeros((Npts, 3))
-        xyz[:, 0] = (lunar_radius + height) * np.cos(lat) * np.cos(lon)
-        xyz[:, 1] = (lunar_radius + height) * np.cos(lat) * np.sin(lon)
-        xyz[:, 2] = (lunar_radius + height) * np.sin(lat)
 
-        xyz = np.squeeze(xyz)
-
-        self = xyz.ravel().view(cls._location_dtype, cls).reshape(xyz.shape[:-1])
-        self._unit = u.meter
-        inst = self.to(height.unit)
+        # get selenocentric coordinates. Have to give one-dimensional array.
 
-        return inst
+        selenodetic = SELENOIDS[ellipsoid](lon, lat, height, copy=False)
+        xyz = selenodetic.to_cartesian().get_xyz(xyz_axis=-1) << height.unit
+        self = xyz.view(cls._location_dtype, cls).reshape(selenodetic.shape)
+        self.ellipsoid = ellipsoid
+        return self
 
     def __str__(self):
         return self.__repr__()
@@ -351,11 +421,19 @@ def selenodetic(self):
         """Convert to selenodetic coordinates."""
         return self.to_selenodetic()
 
-    def to_selenodetic(self):
+    @property
+    def ellipsoid(self):
+        """The default ellipsoid used to convert to selenodetic coordinates."""
+        return self._ellipsoid
+
+    @ellipsoid.setter
+    def ellipsoid(self, ellipsoid):
+        self._ellipsoid = _check_ellipsoid(ellipsoid)
+
+    def to_selenodetic(self, ellipsoid=None):
         """Convert to selenodetic coordinates (lat, lon, height).
 
-        Height is in reference to a sphere with radius `_lunar_radius`,
-        centered at the center of mass.
+        Height is in reference to the ellipsoid.
 
         Returns
         -------
@@ -364,16 +442,15 @@ def to_selenodetic(self):
             `~astropy.coordinates.Latitude`, and `~astropy.units.Quantity`
 
         """
+        ellipsoid = _check_ellipsoid(ellipsoid, default=self.ellipsoid)
         xyz = self.view(self._array_dtype, u.Quantity)
-        lld = CartesianRepresentation(xyz, xyz_axis=-1, copy=False).represent_as(
-            SphericalRepresentation
+        llh = CartesianRepresentation(xyz, xyz_axis=-1, copy=False).represent_as(
+            SELENOIDS[ellipsoid],
         )
-        return GeodeticLocation(
-            Longitude(lld.lon, u.degree, wrap_angle=180.0 * u.degree, copy=False),
-            Latitude(lld.lat, u.degree, copy=False),
-            u.Quantity(
-                lld.distance - (self._lunar_radius * u.meter), self.unit, copy=False
-            ),
+        return SelenodeticLocation(
+            Longitude(llh.lon, u.degree, wrap_angle=180.0 * u.degree, copy=False),
+            Latitude(llh.lat, u.degree, copy=False),
+            u.Quantity(llh.height, self.unit, copy=False),
         )
 
     @property
@@ -476,6 +553,8 @@ def __getitem__(self, item):
 
     def __array_finalize__(self, obj):
         super().__array_finalize__(obj)
+        if hasattr(obj, "_ellipsoid"):
+            self._ellipsoid = obj._ellipsoid
 
     def __len__(self):
         if self.shape == ():
diff --git a/lunarsky/spice_utils.py b/lunarsky/spice_utils.py
index 5563028..06611ce 100644
--- a/lunarsky/spice_utils.py
+++ b/lunarsky/spice_utils.py
@@ -12,13 +12,9 @@
 from .data import DATA_PATH
 
 _J2000 = Time("J2000")
-LUNAR_RADIUS = 1737.1e3  # m
 
 TEMPORARY_KERNEL_DIR = tempfile.TemporaryDirectory()
 
-# TODO --> Look for updated kernels and coordinate definitions.
-#          Use https://naif.jpl.nasa.gov/pub/naif/MER/misc/pck/pck00008.tpc for lunar radii
-
 
 def check_is_loaded(search):
     """
@@ -73,7 +69,7 @@ def furnish_kernels():
     return kernel_paths
 
 
-def lunar_surface_ephem(latitude, longitude, station_num=98):
+def lunar_surface_ephem(pos_x, pos_y, pos_z, station_num=98):
     """
     Make an SPK for the point on the lunar surface
 
@@ -81,10 +77,8 @@ def lunar_surface_ephem(latitude, longitude, station_num=98):
 
     Parameters
     ----------
-    latitude: float
-        Mean-Earth frame selenodetic latitude in degrees.
-    longitude: float
-        Mean-Earth frame selenodetic longitude in degrees.
+    pos_x, pos_y, pos_z: float
+        MCMF frame cartesian position in km
     station_num: int
         Station number
 
@@ -95,15 +89,10 @@ def lunar_surface_ephem(latitude, longitude, station_num=98):
     """
     point_id = 301000 + station_num
 
-    lat = np.radians(latitude)
-    lon = np.radians(longitude)
     ets = np.array([spice.str2et("1950-01-01"), spice.str2et("2150-01-01")])
-    pos_mcmf = spice.latrec(
-        LUNAR_RADIUS / 1e3, lon, lat
-    )  # TODO Use MoonLocation instead?
 
     states = np.zeros((len(ets), 6))
-    states[:, :3] = np.repeat(pos_mcmf[None, :], len(ets), axis=0)
+    states[:, :3] = np.repeat([[pos_x], [pos_y], [pos_z]], len(ets), axis=1).T
 
     center = 301
     frame = "MOON_ME"
diff --git a/lunarsky/tests/test_moon.py b/lunarsky/tests/test_moon.py
index 71bf1b9..8decb46 100644
--- a/lunarsky/tests/test_moon.py
+++ b/lunarsky/tests/test_moon.py
@@ -2,8 +2,13 @@
 import copy
 import numpy as np
 import astropy.units as unit
-from astropy.coordinates import Longitude, Latitude
+from astropy.coordinates import (
+    Longitude,
+    Latitude,
+)
+from astropy.tests.helper import quantity_allclose
 from lunarsky import MoonLocation, MoonLocationAttribute, MCMF
+from lunarsky.moon import SELENOIDS
 
 
 class TestsWithObject:
@@ -19,7 +24,7 @@ def setup_method(self):
         self.lat = Latitude([+0.0, 30.0, 60.0, +90.0, -90.0, -60.0, -30.0, 0.0], unit.deg)
         self.h = unit.Quantity([0.1, 0.5, 1.0, -0.5, -1.0, +4.2, -11.0, -0.1], unit.m)
         self.location = MoonLocation.from_selenodetic(self.lon, self.lat, self.h)
-        self.x, self.y, self.z = self.location.to_selenocentric()
+        self.x, self.y, self.z = self.location.selenocentric
 
     def test_input(self):
         cartesian = MoonLocation(self.x, self.y, self.z)
@@ -27,7 +32,12 @@ def test_input(self):
         cartesian = MoonLocation(self.x.value, self.y.value, self.z.value, self.x.unit)
         assert np.all(cartesian == self.location)
         spherical = MoonLocation(self.lon.deg, self.lat.deg, self.h.to(unit.m))
-        assert np.all(spherical == self.location)
+        # TODO this was just np.all(... == ...). Seeing floating point diff since v0.2.2
+        assert quantity_allclose(
+            spherical.mcmf.cartesian.xyz,
+            self.location.mcmf.cartesian.xyz,
+            atol=1e-9 * unit.m,
+        )
 
     def test_invalid(self):
         # incomprehensible by either raises TypeError
@@ -61,6 +71,10 @@ def test_invalid(self):
         with pytest.raises(ValueError):
             MoonLocation.from_selenodetic(self.lon, self.lat, self.h[:5])
 
+        # invalid ellipsoid
+        with pytest.raises(ValueError):
+            MoonLocation.from_selenodetic(self.lon, self.lat, self.h, ellipsoid="INVALID")
+
     def test_attributes(self):
         assert np.allclose(self.location.height, self.h)
         assert np.allclose(self.location.lon, self.lon)
@@ -130,7 +144,8 @@ def test_moonlocation_delete():
     assert check1 == before
 
 
-def test_station_ids():
+@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.
@@ -186,7 +201,11 @@ def test_station_ids():
 
     for gp in lonlatheights:
         lons, lats, heights = np.array(gp).T
-        locs.append(MoonLocation.from_selenodetic(lat=lats, lon=lons, height=heights))
+        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
@@ -209,7 +228,7 @@ def test_station_ids():
         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]))
+            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
diff --git a/lunarsky/tests/test_spice.py b/lunarsky/tests/test_spice.py
index 9cdf724..fcdfc51 100644
--- a/lunarsky/tests/test_spice.py
+++ b/lunarsky/tests/test_spice.py
@@ -108,12 +108,12 @@ def test_topo_kernel_setup():
 
     for filepath in lunarsky.spice_utils.KERNEL_PATHS:
         spice.furnsh(filepath)
-    lat, lon = 30, 20
+    loc = lunarsky.MoonLocation.from_selenodetic(lat="30d", lon="20d")
     station_name, idnum, frame_specs, latlon = lunarsky.spice_utils.topo_frame_def(
-        lat, lon, moon=True
+        loc.lat.deg, loc.lon.deg, moon=True
     )
     statnum = idnum - 1301000
-    lunarsky.topo._spice_setup(lat, lon, [statnum])
+    lunarsky.topo._spice_setup(loc, [statnum])
     try:
         assert lunarsky.spice_utils.check_is_loaded("*{}*".format(idnum))
     finally:
diff --git a/lunarsky/tests/test_transforms.py b/lunarsky/tests/test_transforms.py
index 1ead075..9afa83d 100644
--- a/lunarsky/tests/test_transforms.py
+++ b/lunarsky/tests/test_transforms.py
@@ -6,14 +6,26 @@
 from astropy.utils import IncompatibleShapeError, exceptions
 from astropy.tests.helper import assert_quantity_allclose
 import lunarsky
+from lunarsky.moon import SELENOIDS
 import pytest
 
+try:
+    from astropy.coordinates.angles.utils import angular_separation
+except ImportError:
+    from astropy.coordinates.angle_utilities import angular_separation
+
 # Lunar station positions
-Nangs = 3
+Nangs = 7
 latitudes = np.linspace(-90, 90, Nangs)
 longitudes = np.linspace(0, 360, Nangs)
 latlons_grid = [(lat, lon) for lon in longitudes for lat in latitudes]
 
+# Avoiding poles:
+latlons_grid_nopole = [
+    (lat, lon)
+    for lon in np.linspace(0.1, 360, Nangs, endpoint=False)
+    for lat in np.linspace(-70.0, 70.0, Nangs, endpoint=False)
+]
 
 # Times
 t0 = Time("2020-10-28T15:30:00")
@@ -54,10 +66,11 @@ def test_icrs_to_topo_long_time(time, lat, lon, grcat):
         (jd_4mo[:2], [10.3, 11.2], [0.0, 1.4]),
     ],
 )
-def test_transform_loops(obj, path, time, lat, lon):
+@pytest.mark.parametrize("ell", SELENOIDS)
+def test_transform_loops(obj, path, time, lat, lon, ell):
     # Tests from ICRS -> path -> ICRS
     obj = lunarsky.SkyCoord(obj)  # Ensure we're working with lunarsky-compatible SkyCoord
-    loc = lunarsky.MoonLocation.from_selenodetic(lat, lon)
+    loc = lunarsky.MoonLocation.from_selenodetic(lat, lon, ellipsoid=ell)
     fdict = {
         "lunartopo": lunarsky.LunarTopo(location=loc, obstime=time),
         "mcmf": lunarsky.MCMF(obstime=time),
@@ -75,9 +88,10 @@ def test_transform_loops(obj, path, time, lat, lon):
     assert_quantity_allclose(obj.cartesian.xyz, orig_pos, atol=tol)
 
 
-def test_topo_to_topo():
+@pytest.mark.parametrize("ell", SELENOIDS)
+def test_topo_to_topo(ell):
     # Check that zenith source transforms properly
-    loc0 = lunarsky.MoonLocation.from_selenodetic(lat=0, lon=90)
+    loc0 = lunarsky.MoonLocation.from_selenodetic(lat=0, lon=90, ellipsoid=ell)
     loc1 = lunarsky.MoonLocation.from_selenodetic(lat=0, lon=0)
 
     src = lunarsky.SkyCoord(alt="90d", az="0d", frame="lunartopo", location=loc0)
@@ -97,14 +111,13 @@ def test_mcmf_to_mcmf():
     )
     sph0 = src.spherical
     sph1 = orig_pos.spherical
-    res = ac.angle_utilities.angular_separation(
-        sph0.lon, sph0.lat, sph1.lon, sph1.lat
-    ).to("deg")
+    res = angular_separation(sph0.lon, sph0.lat, sph1.lon, sph1.lat).to("deg")
     assert_quantity_allclose(res, 177 * un.deg, atol=5 * un.deg)
     # Tolerance to allow for lunar precession / nutation.
 
 
-def test_earth_from_moon():
+@pytest.mark.parametrize("ell", SELENOIDS.keys())
+def test_earth_from_moon(ell):
     # Look at the position of the Earth in lunar reference frames.
 
     # The lunar apo/perigee vary over time. These are
@@ -133,7 +146,7 @@ def test_earth_from_moon():
 
     # Now test LunarTopo frame positions
     lat, lon = 0, 0  # deg
-    loc = lunarsky.MoonLocation.from_selenodetic(lat, lon)
+    loc = lunarsky.MoonLocation.from_selenodetic(lat, lon, ellipsoid=ell)
     topo = mcmf.transform_to(lunarsky.LunarTopo(location=loc, obstime=jd_4mo))
 
     # The Earth should stay within 10 deg of zenith of lat=lon=0 position
@@ -238,7 +251,7 @@ def test_incompatible_shape_error(Nt, Nl, success):
     if success:
         lunarsky.LunarTopo(location=locs, obstime=times)
     else:
-        with pytest.raises(ValueError, match="non-scalar attributes"):
+        with pytest.raises(ValueError, match="non-scalar data and/or attributes"):
             lunarsky.LunarTopo(location=locs, obstime=times)
 
 
@@ -268,22 +281,75 @@ def test_incompatible_transform(fromframe):
         src.transform_to(ltop)
 
 
-def test_finite_vs_spherical():
-    # Transform MCMF coordinates with distance and without
+@pytest.mark.parametrize("ell", SELENOIDS)
+def test_finite_vs_spherical(ell):
+    # Transform MCMF coordinates with distance and without units
+    # Check consistency with ellipsoid equatorial radius
     # 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)
+
+    R0 = 404789  # km
+    xyz = np.array([[R0, -R0], [0, 0], [0, 0]])
+    with_units = lunarsky.SkyCoord(lunarsky.MCMF(*(xyz * un.km)))
+    sans_units = lunarsky.SkyCoord(lunarsky.MCMF(*(xyz)))
+
+    loc = lunarsky.MoonLocation.from_selenodetic(
+        lon=180 * un.deg, lat=0, height=0, ellipsoid=ell
+    )
     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)
+    dists = R0 * un.km + np.array([1, -1]) * SELENOIDS[ell]._equatorial_radius
+    assert np.all(altaz_with_units.distance == dists)
     assert_quantity_allclose(altaz_sans_units.distance, 1.0)
+
+
+@pytest.mark.parametrize("ell", SELENOIDS)
+@pytest.mark.parametrize("lat,lon", latlons_grid_nopole)
+def test_topo_zenith_shift(ell, lat, lon):
+    # Verify that a given source at zenith in one topo frame shifts
+    # predictably when viewed from a topo frame with the same lat/lon but
+    # different ellipsoid
+
+    # Checking that the ellipsoid is interpreted correctly
+
+    # This test is a little sketchy... will need to review this later.
+    #   Some discrepancies for GRAIL23 selenoid when the source distance is large.
+    #       Choosing 1000 km for now.
+    #   Also fails near poles.
+
+    # Comparing against the SPHERE ellipsoid. Test fails for this due to divide by zero
+    if ell == "SPHERE":
+        return
+
+    lat *= un.deg
+    lon *= un.deg
+
+    loc0 = lunarsky.MoonLocation.from_selenodetic(
+        lon, lat, ellipsoid="SPHERE"
+    )  # For reference.
+    loc1 = lunarsky.MoonLocation.from_selenodetic(
+        lon, lat, ellipsoid=ell
+    )  # Same lat/lon = different place for different ellipsoid
+
+    # Zenith source at finite distance over loc0.
+    src0 = lunarsky.SkyCoord(
+        alt="90d",
+        az="0d",
+        distance=1e3 * un.km,
+        frame=lunarsky.LunarTopo(location=loc0, obstime=Time.now()),
+    )
+    src1 = src0.transform_to(lunarsky.LunarTopo(location=loc1))
+
+    # Law of sines:
+    #   Geometry among zenith angle, mcmf station vector,
+    #   and selenocentric vs. selenodetic latitudes
+
+    R1 = loc1.mcmf.cartesian.norm()
+    lat_cen = loc1.mcmf.spherical.lat
+    lat_det = loc1.lat
+    assert_quantity_allclose(
+        src1.distance / np.sin((np.abs(lat_det - lat_cen)).rad),
+        R1 / np.sin(src1.zen.rad),
+        atol=1 * un.km,
+    )
diff --git a/lunarsky/topo.py b/lunarsky/topo.py
index d45def0..bc39ff4 100644
--- a/lunarsky/topo.py
+++ b/lunarsky/topo.py
@@ -88,30 +88,25 @@ def zen(self):
 # -----------------
 
 
-def _spice_setup(latitude, longitude, station_id):
-
-    latlonids = np.stack(
-        [np.atleast_1d(latitude), np.atleast_1d(longitude), station_id]
-    ).T
-    if latlonids.ndim == 1:
-        latlonids = latlonids[None, :]
-
-    for lat, lon, sid in latlonids:
-        sid = int(sid)  # Station IDs must be ints, but are converted to float above.
+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(
-                lat, lon, station_num=sid
+                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(
-                lat, lon, moon=True, station_num=sid
+                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 make_transform(coo, toframe):
-
     ap_to_spice = {"icrs": ("J2000", 0), "mcmf": ("MOON_ME", 301)}
     # Get target frame and source frame names
     if isinstance(coo, LunarTopo):
@@ -138,7 +133,7 @@ def make_transform(coo, toframe):
     shape_out = check_broadcast(coo.shape, ets.shape, location.shape)
 
     # Set up SPICE ephemerides and frame details
-    _spice_setup(location.lat.deg, location.lon.deg, stat_ids)
+    _spice_setup(location, stat_ids)
 
     ets_ids = np.atleast_2d(np.stack(np.broadcast_arrays(ets, stat_ids)).T)
 
@@ -205,13 +200,11 @@ def make_transform(coo, toframe):
 # -----------------
 @frame_transform_graph.transform(FunctionTransformWithFiniteDifference, ICRS, LunarTopo)
 def icrs_to_lunartopo(icrs_coo, topo_frame):
-
     return make_transform(icrs_coo, topo_frame)
 
 
 @frame_transform_graph.transform(FunctionTransformWithFiniteDifference, LunarTopo, ICRS)
 def lunartopo_to_icrs(topo_coo, icrs_frame):
-
     return make_transform(topo_coo, icrs_frame)
 
 
@@ -220,13 +213,11 @@ def mcmf_to_lunartopo(mcmf_coo, topo_frame):
     # TODO:
     #   > What if mcmf_coo and topo_frame have different obstimes?
     #   > What if location has obstime?
-
     return make_transform(mcmf_coo, topo_frame)
 
 
 @frame_transform_graph.transform(FunctionTransformWithFiniteDifference, LunarTopo, MCMF)
 def lunartopo_to_mcmf(topo_coo, mcmf_frame):
-
     return make_transform(topo_coo, mcmf_frame)
 
 
diff --git a/setup.py b/setup.py
index 4fc2b8f..60461b5 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>3.0", "spiceypy", "jplephem"],
+    "install_requires": ["numpy>=1.15", "astropy>=6.0.0", "spiceypy", "jplephem"],
     "classifiers": [
         "Development Status :: 3 - Alpha",
         "Intended Audience :: Science/Research",