Position of the Moon (topocentric) using Astrometry.jl (SOFA/ERFA) routines

[tamasgal]

OK I finally ask for help here because I am losing my head

I want to calculate the Moon coordinates in an equatorial coordinate system for an observer on Earth at a given time using Astrometry.jl which uses SOFA/ERFA routines..

What I got so far is calculating the Moon coordinates in a geocentric coordinate system using this tiny function

function moon(dt::DateTime)
    ht = 0.0  # sitting on Earht's surface
    jd = (datetime2julian(dt), 0.0)
    tt = SOFA.taitt(SOFA.utctai(jd...)...)
    pm, vm = SOFA.moon98(tt...)
    rc, dc = SOFA.c2s(pm)  # cartesian to spherical
    ra, dec = mod2pi(rc), dc
    return (ra=ra, dec=dec)
end

and it matches AstroPy’s get_body(...) when called without a location.

Fine, the problem I have now is to calculate the coordinates for a given location, which means I need to calculate it in topocentric system.

Using AstroPy (btw. I get the same results using SLALIB’s RDPLAN routine) I obtain the following values for a given location and time:

from astropy.coordinates import EarthLocation, AltAz, get_body
from astropy.time import Time
import astropy.units as u

location = EarthLocation(
    lat=0.7470155743891863*u.rad,
    lon=0.10510862459055528*u.rad,
    height=0.0*u.m
)

d = Time(["2028-11-08T23:34:54.00"], scale="utc")

print(get_body("moon", d, location))
WARNING: Tried to get polar motions for times after IERS data is valid. Defaulting to polar motion from the 50-yr mean for those. This may affect precision at the arcsec level. Please check your astropy.utils.iers.conf.iers_auto_url and point it to a newer version if necessary. [astropy.coordinates.builtin_frames.utils]

<SkyCoord (GCRS: obstime=['2028-11-08T23:34:54.000'], obsgeoloc=[(3144039.81787578, 3484619.75485466, 4300693.72117526)] m, obsgeovel=[(-254.1002625, 228.37771983, 0.71891401)] m / s): (ra, dec, distance) in (deg, deg, AU)
    [(127.74528269, 16.71982861, 0.00250404)]>

So we get

• RA = 127.74528269
• DEC = 16.71982861

Note that our moon function above yields coordinates roughly 1deg off due to the parallax of the Moon (it’s in the geocentric system and not taking account any location on Earth)

julia> dt = DateTime("2028-11-08T23:34:54.00")
2028-11-08T23:34:54

julia> pos = moon(dt)
(ra = 2.2167352184921754, dec = 0.30209193883117474, azimuth = 1.4716666159254528, elevation = 0.34136871844980865)

julia> rad2deg.([pos.ra, pos.dec])
2-element Vector{Float64}:
 127.00957231761205
  17.30859311995054

Now comes my function (see below) which tries to do everything right, but doesn’t, and I can’t figure out where the problem is. I tried to comment everything so that it’s maybe a bit easier to spot the error.

This is what I get

julia> pos = moon_topocentric(lon, lat, dt)
(ra = 2.21935735847888, dec = 0.299922564561618, distance_au = 0.0025181965734329967)

julia> rad2deg.([pos.ra, pos.dec])
2-element Vector{Float64}:
 127.15980987214272
  17.18429713012066

Which is still off by around half a degree or so!

Does anybody know what’s wrong here? I got great help two years ago from @DaveMacMahon here Calculate positions of astronomical objects from a terrestrial observatory where I calculated the Sun coordinates but that did not involve any parallax corrections.

I currently think that the nutation/precession correction is the last missing piece but cannot get it right.

Here is the function I am sitting at since yesterday

using Astrometry.SOFA
using Dates

function moon_topocentric(lon, lat, dt::DateTime; height=0.0)
    # UTC as JD
    utc1, utc2 = datetime2julian(dt), 0.0
    
    # Earth orientation (use IERS values for high precision)
    dut1 = 0.0
    xp = 0.0
    yp = 0.0
    
    # Disable refraction, we are observing neutrinos so we don't care right now ;)
    phpa = 0.0   # pressure (hPa)
    tc   = 0.0   # temperature (C)
    rh   = 0.0   # relative humidity
    wl   = 0.0   # wavelength (micron)
    
    # Build astrometry context
    astrom, eo = SOFA.apco13(
        utc1, utc2,
        dut1,
        lon, lat, height,
        xp, yp,
        phpa, tc, rh, wl
    )
    
    # TT for Moon ephemeris
    tai1, tai2 = SOFA.utctai(utc1, utc2)
    tt1, tt2   = SOFA.taitt(tai1, tai2)
    
    # Moon geocentric GCRS position (AU and AU/day)
    pm, vm = SOFA.moon98(tt1, tt2)
    
    # Calculate distance in AU
    distance_au = sqrt(pm[1]^2 + pm[2]^2 + pm[3]^2)
    
    # Parallax in arcseconds: π = 1 AU / distance
    parallax_arcsec = 1.0 / distance_au
    
    # Cartesian -> spherical (GCRS)
    ra_gc, dec_gc = SOFA.c2s(pm)
    
    # Calculate radial velocity (AU/day)
    rv_au_day = (pm[1]*vm[1] + pm[2]*vm[2] + pm[3]*vm[3]) / distance_au
    # Convert to km/s (1 AU/day ≈ 1731.46 km/s)
    rv_km_s = rv_au_day * 1731.46
    
    # GCRS -> topocentric apparent place with parallax
    # Note: For Moon, proper motion during the observation is negligible
    # compared to parallax effect, so we can use 0.0 for pmra/pmdec
    ra_topo, dec_topo = SOFA.atciq(
        ra_gc, dec_gc,
        0.0, 0.0,           # proper motion (negligible for Moon)
        parallax_arcsec, 
        rv_km_s,            # radial velocity (km/s)
        astrom
    )
    
    return (
        ra  = mod2pi(ra_topo),
        dec = dec_topo,
        distance_au = distance_au
    )
end

[tamasgal]

I added the alt/az calculation and compared to AstroPy by using AstroPy to generate a dataset of 24 hours.

What you can see is the moon function (geocentric) which ignores the Earth location and the slightly extended moon_topocentric function (topocentric). The differences are periodic and it seems to me that I am not at the right “location” on Earth.

The azimuth and right ascension are correct twice a day

Here is the current code

function moon_topocentric2(lon, lat, dt::DateTime; height=0.0)
    # UTC as JD
    utc1, utc2 = datetime2julian(dt), 0.0
    
    # Earth orientation (use IERS values for high precision)
    dut1 = 0.0
    xp = 0.0
    yp = 0.0
    
    # Disable refraction
    phpa = 0.0   # pressure (hPa)
    tc   = 0.0   # temperature (C)
    rh   = 0.0   # relative humidity
    wl   = 0.0   # wavelength (micron)
    
    # Build astrometry context
    astrom, eo = SOFA.apco13(
        utc1, utc2,
        dut1,
        lon, lat, height,
        xp, yp,
        phpa, tc, rh, wl
    )
    
    # TT for Moon ephemeris
    tai1, tai2 = SOFA.utctai(utc1, utc2)
    tt1, tt2   = SOFA.taitt(tai1, tai2)
    
    # Moon geocentric GCRS position (AU and AU/day)
    pm, vm = SOFA.moon98(tt1, tt2)
    
    # Calculate distance in AU
    distance_au = sqrt(pm[1]^2 + pm[2]^2 + pm[3]^2)
    
    # Parallax in arcseconds: π = 1 AU / distance
    parallax_arcsec = 1.0 / distance_au
    
    # Cartesian → spherical (GCRS)
    ra_gc, dec_gc = SOFA.c2s(pm)
    
    # Calculate radial velocity (AU/day)
    rv_au_day = (pm[1]*vm[1] + pm[2]*vm[2] + pm[3]*vm[3]) / distance_au
    # Convert to km/s (1 AU/day ≈ 1731.46 km/s)
    rv_km_s = rv_au_day * 1731.46
    
    # GCRS → topocentric apparent place with parallax
    # Note: For Moon, proper motion during the observation is negligible
    # compared to parallax effect, so we can use 0.0 for pmra/pmdec
    ra_topo, dec_topo = SOFA.atciq(
        ra_gc, dec_gc,
        0.0, 0.0,           # proper motion (negligible for Moon)
        parallax_arcsec,    # parallax (arcsec) - THIS IS THE KEY!
        rv_km_s,            # radial velocity (km/s)
        astrom
    )
    
    moon_dist_au = sqrt(sum(pm.^2))
    # 180 × 3600 / π, for reference ;)
    parallax = 1.0 / (moon_dist_au * 206264.806247)  # parallax in arcseconds
    moonobs = SOFA.atco13(ra_topo, dec_topo, 0.0, 0.0, parallax, 0.0, utc1, utc2, 0.0, lon, lat, height, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

    az, el = moonobs.azi, π/2 - moonobs.zen

    return (ra=mod2pi(ra_topo), dec=dec_topo, distance_au = distance_au, azimuth=az, elevation=el)
end

[cormullion]

Can NASA software help?

Horizons System

[mihalybaci]

Piggybacking off this, I would say that Python’s skyfield package is much better designed for ephemeris calculations than even Astropy is. It might be good to take inspiration from that package, it says it is “pure python” so the calculations may be exposed in a way you could translate to Julia.

[tamasgal]

Thanks @mihalybaci I will have a look but as a side note; we also have a Slalib based package which agrees with the dataset generated by AstroPy within 0.0001 degree. The fact that I cannot get the Moon right within 0.7deg is very likely a user mistake on my side

[tamasgal]

I will check this out as well, as a fourth reference

[mihalybaci]

Yeah, it’s tricky and I am certainly not well-versed in planetary motion. If I had to guess, and I don’t, I would take a long look at the lunar parallax again as you already mentioned. In your code, the parallax is simply a function of the distance, but the parallax will be a function of time and location on the Earth as well. Your code may already take that into account and I missed it.

Edit: It looks like SPICE.jl may be the current best way to implement this kind of calculation, unless you do want/need to roll your own.

[tamasgal]

I am now going through my code and the SOFA docs again; I must have made mistakes, as you also pointed out, probably with the parallax calculation.

I would like to stick to Astrometry.jl due to the pure Julia and basically zero-dependency nature of it. It has to work somehow

[DaveMacMahon]

moon98 returns the position relative to GCRS. I think you can use the c2t06a function to get a rotation matrix that you can left multiply the position vector by to convert it to the International Terrestrial Reference System. Then I think it’s just longitude/latitude rotations to get to local topocentric vector, which can be converted to azimuth and elevation.

I’ll see if I can work up a more specific example this evening.

[DaveMacMahon]

Here is a script showing how to compute the topocentric azimuth and elevation of the Moon for a given observer location and date/time. This script uses the Green Bank Telescope as the observer location and 2028-11-08T23:34:54 as the time, but it should be pretty easy to change as desired. The script also includes some hard coded verification values derived from JPL HORIZONS. This script could be adapted to use HORIZONS.jl to dynamically generate the verification values from HORIZONS for whichever location and time the user specifies.

Here is the output from the script:

┌ Info: moon geocentric RA,dec at 2028-11-08T23:34:54
│   geora = 127.00957231761205
│   jplra = 127.007799642
│   geodec = 17.308593119950537
└   jpldec = 17.308994763
┌ Info: moon az,el at 2028-11-08T23:34:54 for observer location
│   obsaz = 17.218596681255377
│   jplaz = 17.219373403
│   obsel = -33.37961826856053
└   jplel = -33.379180688

The observer in Green Bank might be disappointed that the moon is still below the horizon.

Here is the script:

using Dates
using Astrometry.SOFA: utctai, taitt, moon98, AU, c2t06a, gd2gc, c2s, hd2ae

#---
# Setup observer location (Green Bank Telescope according to JPL HORIZONS)

lon = deg2rad(280.1601614)
lat = deg2rad(38.433129)
htm = 823.664

#---
# Setup time of observation

dt = DateTime("2028-11-08T23:34:54.000")
jdutc = datetime2julian(dt)
jdtai = utctai(jdutc, 0.0)
jdtt = taitt(jdtai...)

#---
# Get geocentric moon position and velocity at time of observation (in m and m/d)

moonp, moonv = moon98(jdtt...) .* AU

#---
# Get RA/dec in degrees for moon geocentric position at time of observation
# This is only used to compare with other references (e.g. JPL HORIZONS).
#
# Calculated: ra=127.00957231761205 dec=17.308593119950537
# HORIZONS:   ra=127.007799642      dec=17.308994763

geora, geodec = rad2deg.(c2s(moonp))
jplra, jpldec = 127.007799642, 17.308994763
@info "moon geocentric RA,dec at $dt" geora jplra geodec jpldec

#---
# Get rotation matrix to convert from GCRS to ITRS.  By using UTC we are
# implicitly ignoring `dut1`.  We explicilty ignore polar motion.

rc2t = c2t06a(jdtt..., jdutc, 0.0, 0.0, 0.0)

#---
# Rotate geocentric moon position vector into ITRS frame

gcmoonitrs = rc2t * moonp

#---
# Get geocentric ITRS position vector from geodetic observer location

gcobsitrs = gd2gc(:WGS84, lon, lat, htm)

#---
# Subtract geocentric observer ITRS vecor from geocentrim moon ITRS vector to
# get observer-to-moon ITRS vector.

obsmoonitrs = gcmoonitrs - gcobsitrs

#---
# Get hour angle and declination of moon from observer-to-moon ITRS vector

ha, dec = c2s(obsmoonitrs)

#---
# Convert (hour angle, dec) to (az, el)
#
# There are two tricky bits here.  One is that the "hour angle" from `c2s` has
# the opposite sign convention from hour angle.  The other is that it is
# relative to the prime meridian, so we have to adjust for the observer's
# longitude to make it local.  This is why we pass `lon-ha` as the hour angle
# argument to `hd2ae`.
#
# For comparison to JPl HORIZONS:
#
# Calculated: az=17.218596681255377 el=-33.37961826856053
# HORIZONS:   az=17.219373403       el=-33.379180688

obsaz, obsel = rad2deg.(Tuple(hd2ae(lon-ha, dec, lat)))
jplaz, jplel = 17.219373403, -33.379180688

@info "moon az,el at $dt for observer location" obsaz jplaz obsel jplel

[tamasgal]

Awesome @DaveMacMahon azimuth and elevation are spot on, many many thanks!

I only have one problem left: how do I get the topocentric right ascension and declination? I thought (given your comments regarding the tricky part of hd2ae) I could use lon - ha and dec but the values are off.

The declination is off by a seemingly constant offset of roughly 0.5 deg (seems to vary, depending on the date, but it’s round 0.5 deg) in my example and the right ascension seems to rotate away

That’s my return statement using your function (I made sure to return rad-values)

    return (ra=mod2pi(lon-ha), dec=dec, azimuth=obsaz, elevation=obsel)

This is the input of my dataset

 longitude: 6.022280579463798 deg (0.10510862459055528)
 latitude: 42.800839643041364 deg (0.7470155743891863)
height: 0.0
datetime: 2022-03-02T18:20:35.300000   # 1min steps, over 24 hours

here the plot, comparing to AstroPy (the “geocentric” is simply my calculation, which is obvously wrong, but I left it there for comparison)

[tamasgal]

Ah, I got it!

Here is the final function. I am using rc2t to rotate back to GCRS and then that one to calculate ra and dec from the resulting obsmoon_gcrs (which is in the inertial frame)

I left and put some comments, I hope I got everything right. Thanks again @DaveMacMahon !

function moon_dave(lon, lat, dt::DateTime; height=0.0)
    jdutc = datetime2julian(dt)
    jdtai = utctai(jdutc, 0.0)
    jdtt = taitt(jdtai...)
    
    # geocentric moon position and velocity
    moonp, moonv = moon98(jdtt...) .* AU
    
    # rot matrix to convert from GCRS to ITRS
    rc2t = c2t06a(jdtt..., jdutc, 0.0, 0.0, 0.0)
    
    # rotate geocentric moon position vector into ITRS frame
    gcmoonitrs = rc2t * moonp
    
    # geocentric ITRS position vector from geodetic observer location
    gcobsitrs = gd2gc(:WGS84, lon, lat, height)
    
    # observer-to-moon ITRS vector
    obsmoonitrs = gcmoonitrs - gcobsitrs
    
    # rotate topocentric vector back to GCRS to get RA/Dec**
    obsmoon_gcrs = rc2t' * obsmoonitrs
    
    # ra and dec from the GCRS topocentric vector
    ra, dec = c2s(obsmoon_gcrs)
    
    # hour angle and declination from ITRS vector for az/el calculation
    ha, dec_itrs = c2s(obsmoonitrs)
    
    # convert (hour angle, dec) to (az, el) using the "tricky" hd2ae function
    obsaz, obsel = hd2ae(lon-ha, dec_itrs, lat)
    
    return (ra=mod2pi(ra), dec=dec, azimuth=obsaz, elevation=obsel)
end

Comparision (the final version called “Dave alt”) to AstroPy, an also the initial “geocentric” version:

Final results, comparing only AstroPy

[tamasgal]

Btw. wondering where the final wobble difference is coming from…

[DaveMacMahon]

Using c2s to get RA/dec from obsmoon_gcrs results in ra and dec being in the GCRS frame. When comparing these values with other RA/dec values from other sources, it’s important to ensure that both RA/dec pairs are in the same frame. Do you know which frame your other (e.g. Python derived) RA/dec pair is in?

Continuing from my earlier example, the following code computes the moon RA/dec in the ICRS frame as seen from the observer location and compares with the value from JPL HORIZONS:

#---
# Get obsmoon vector in GCRS frame

obsmoongcrs = rc2t' * obsmoonitrs

#---
# Get obsmoon RA/dec in GCRS frame

obsmoongcrsra, obsmoongcrsdec = c2s(obsmoongcrs)

#---
# Convert obsmoongcrsra, obsmoongcrsdec to ICRS

astrom = SOFA.apcg13(jdtt...)
obsmoonicrsra, obsmoonicrsdec = SOFA.aticq(obsmoongcrsra, obsmoongcrsdec, astrom)

#---
# Convert to degrees and display alongside HORIZONS values

obsmoongcrsradeg, obsmoobgcrsdecdeg = rad2deg.((obsmoongcrsra, obsmoongcrsdec))
obsmoonicrsradeg, obsmoobicrsdecdeg = rad2deg.((obsmoonicrsra, obsmoonicrsdec))
jplmoonicrsradeg, jplmoonicrsdecdeg = 127.215173462, 16.526197595

@info("moon RA,dec at $dt from observer location",
    obsmoongcrsradeg, obsmoonicrsradeg, jplmoonicrsradeg,
    obsmoobgcrsdecdeg, obsmoobicrsdecdeg, jplmoonicrsdecdeg
)

which produces:

┌ Info: moon RA,dec at 2028-11-08T23:34:54 from observer location
│   obsmoongcrsradeg = 127.21691359389354
│   obsmoonicrsradeg = 127.21564743505404
│   jplmoonicrsradeg = 127.215173462
│   obsmoobgcrsdecdeg = 16.525870840471995
│   obsmoobicrsdecdeg = 16.525924620360865
└   jplmoonicrsdecdeg = 16.526197595

[mihalybaci]

Awesome that you’ve got this working! Do you intend to make this a registered package? It would be pretty awesome to have this alongside SolarPosition.jl by @langestefan.