Are round-trip conversions between geocentric and geodetic EarthLocation expected to be stable?

Astropy
,
1

Hello, I’ve run into an issue with Astropy causing test failures for us. I’m not sure if this is expected behavior / something I’m doing wrong, but when converting between geocentic and geodetic EarthLocation coordinates, I don’t always get the same answer when converting from one reference frame to the other and back again.

Here’s a small self-contained example:


import random
from astropy.coordinates import EarthLocation


def try_roundtrip():
    rand = lambda: random.random() * 10e6
    before = EarthLocation.from_geocentric(
        x=rand(), y=rand(), z=rand(), unit="m"
    )
    geodetic = before.to_geodetic()
    after = EarthLocation.from_geodetic(
        lat=geodetic.lat,
        lon=geodetic.lon,
        height=geodetic.height,
    )
    assert before == after, (
        f"initial value {before} doesn't match conversion round-trip value {after}"
    )


try_roundtrip()

The assertion fails for me often, but not every single time. This is on an Ubuntu 24.04 workstation in Python 3.12 using Astropy 7.1.1.

The numbers I get back are close, but not exactly the same. I’m assuming there’s some floating point rounding happening internally somewhere. Is there any way to guarantee stable output in tests, or do I need to add pytest.approx() and live with it?

2

You should generally never use == for floating point comparisons, which will readily fail already (on the face) trivial tests like math.sqrt(2)**2 == 2. np.allclose or, for quantities like here, units.allclose provide better control of the expected precision, e.g.

from astropy import units as u
u.allclose([after.x, after.y, after.z], [before.x, before.y, before.z], rtol=1e-10)

this should still round-trip for rtol well below the default of 1e-5.