Determine why the Moon's AltAz coordinates vary across codes

@mhvk

This comes out of the discussion in #4254, as well as some offline discussion about astroplan development.

If you execute the gist at https://gist.github.com/eteq/08dba433082e57b7f2af (with astropy master, pyephem, and skyfield installed), you see this:

PyEphem (Alt, Az, Dist): (55.4920153458 deg, 239.031895666 deg, 356424.848107 km)
Skyfield (Alt, Az, Dist): (55.32651159 deg, 239.298837015 deg, 356413.96585 km)
Skyfield(de430) (Alt, Az, Dist): (55.3265117716 deg, 239.298837058 deg, 356413.965628 km)
jplephem/Astropy (Alt, Az, Dist): (56.0578669005 deg, 239.043478848 deg, 356423.493806 km)
jplephem/Astropy(ICRS) (Alt, Az, Dist): (56.0596259457 deg, 239.03677597 deg, 356421.587215 km)

That demonstrates that there are disagreements about the as-observed location of the moon at the ~0.5-1 degree level (and ~10 km distance), but between astropy and the other two, as well as between skyfield and pyephem.

It's quite possible this level of difference is driven by subtle differences in how the various frames are defined rather than any sort of implementation bug... But we should try to get to the bottom of the differences, if possible.

And once we've figured it out, even if there are no bugs to be fixed, hopefully this investigation can provide a thorough battery of tests to be added to the testing suite.

cc @mhvk @bmorris3 @brandon-rhodes

@eteq

Also @eteq – I see that 1.0.6 was released, but I can't grab it with conda update astropy. Should it be on conda soon?

@bmorris3

@bmorris3 - sometimes it takes them a little while to update it... This has gotten to be a while though, so I should probably ping them.

In the meantime, you can always

conda create -n astropy1.0.6 astropy
source activate astropy1.0.6

pip install --upgrade astropy==1.0.6
<or>
cd /path/to/astropy
git checkout v1.0.6
python setup.py install

(or just do a manual install

Skyfield should always agree with the highest precision results from HORIZONS. While PyEphem does not offer that precision, I would have expected it to agree within an arcsecond or less for an object whose position is known with as high a precision as the Moon’s.

But in this case — to focus on the azimuth, since it is not affected by refraction and so involves fewer variables — the difference is an appalling 961 arcseconds. What is wrong?

A first problem is that PyEphem does not understand time scales at all. It only knows Terrestrial Time. So the statement:

        obs.date = time.datetime

should in fact be:

        obs.date = time.tt.datetime

This will bring PyEphem’s prediction to within about 53 arcseconds of Skyfield — which, while it is still quite appallingly bad, is at least within one arcminute. I will look for smaller effects that might reduce the difference even more, and report back!

@brandon-rhodes

Thanks for that clarification @brandon-rhodes, I wouldn't have thought of that.

So that I do not run in circles and re-investigate topics over again in a few days when I have time to return to this problem, here are two effects which are not large enough to explain the difference:

When the Skyfield topos() method is called, no elevation_m keyword is provided to let Skyfield know the observer’s elevation. I would have thought that a difference in position of 2,798 m would make more of a difference, but adding in the keyword argument makes a difference of only 6.27e-05 arcseconds (!) to the azimuth.
The C library inside PyEphem thinks that the Earth’s flatness is 298.257 whereas Skyfield knows the 2010 estimate of 298.25642. It makes a difference of 3.31e-05 arcseconds.
I wondered whether the Moon model used in PyEphem might be less accurate than DE421, but a comment in moon.c says Precision < 0.05" from -1369 to +2950 and if by " they mean arcseconds then, nope, that is not the source of the difference we are seeing here!

Instead of simply comparing both libraries to each other, it is best to ask HORIZONS for these values and then figure out why each library diverges from that best-available number. So here goes a big cut-and-paste:

*******************************************************************************
Ephemeris / WWW_USER Mon Nov  2 06:14:49 2015 Pasadena, USA      / Horizons    
*******************************************************************************
Target body name: Moon (301)                      {source: DE-0431LE-0431}
Center body name: Earth (399)                     {source: DE-0431LE-0431}
Center-site name: (user defined site below)
*******************************************************************************
Start time      : A.D. 1984-Jan-21 12:00:00.0000 UT      
Stop  time      : A.D. 1984-Jan-22 12:00:00.0000 UT      
Step-size       : 1440 minutes
*******************************************************************************
Target pole/equ : IAU_MOON                        {East-longitude +}
Target radii    : 1737.4 x 1737.4 x 1737.4 km     {Equator, meridian, pole}    
Center geodetic : 254.180000,32.7800000,2.7980000 {E-lon(deg),Lat(deg),Alt(km)}
Center cylindric: 254.180000,5370.07631,3434.9857 {E-lon(deg),Dxy(km),Dz(km)}
Center pole/equ : High-precision EOP model        {East-longitude +}
Center radii    : 6378.1 x 6378.1 x 6356.8 km     {Equator, meridian, pole}    
Target primary  : Earth
Vis. interferer : MOON (R_eq= 1737.400) km        {source: DE-0431LE-0431}
Rel. light bend : Sun, EARTH                      {source: DE-0431LE-0431}
Rel. lght bnd GM: 1.3271E+11, 3.9860E+05 km^3/s^2                              
Atmos refraction: NO (AIRLESS)
RA format       : DEG
Time format     : CAL 
EOP file        : eop.151030.p160121                                           
EOP coverage    : DATA-BASED 1962-JAN-20 TO 2015-OCT-30. PREDICTS-> 2016-JAN-20
Units conversion: 1 au= 149597870.700 km, c= 299792.458 km/s, 1 day= 86400.0 s 
Table cut-offs 1: Elevation (-90.0deg=NO ),Airmass (>38.000=NO), Daylight (NO )
Table cut-offs 2: Solar Elongation (  0.0,180.0=NO ),Local Hour Angle( 0.0=NO )
*****************************************************************************************************
 Date__(UT)__HR:MN     R.A.___(J2000.0)___DEC. R.A._(airls-appar)_DEC. Azi_(a-appr)_Elev        CT-UT
*****************************************************************************************************
$$SOE
 1984-Jan-21 12:00  m  164.7068409  11.5136952 164.4980406  11.5991097 239.0335  55.4910    54.184520
 1984-Jan-22 12:00  m  178.5589980   5.5460418 178.3537413   5.6347851 213.6816  58.6372    54.184548
$$EOE
*****************************************************************************************************
Column meaning:

TIME

  Prior to 1962, times are UT1. Dates thereafter are UTC. Any 'b' symbol in
the 1st-column denotes a B.C. date. First-column blank (" ") denotes an A.D.
date. Calendar dates prior to 1582-Oct-15 are in the Julian calendar system.
Later calendar dates are in the Gregorian system.

  Time tags refer to the same instant throughout the universe, regardless of
where the observer is located.

  The dynamical Coordinate Time scale is used internally. It is equivalent to
the current IAU definition of "TDB". Conversion between CT and the selected
non-uniform UT output scale has not been determined for UTC times after the
next July or January 1st.  The last known leap-second is used over any future
interval.

  NOTE: "n.a." in output means quantity "not available" at the print-time.

SOLAR PRESENCE (OBSERVING SITE)
  Time tag is followed by a blank, then a solar-presence symbol:

        '*'  Daylight (refracted solar upper-limb on or above apparent horizon)
        'C'  Civil twilight/dawn
        'N'  Nautical twilight/dawn
        'A'  Astronomical twilight/dawn
        ' '  Night OR geocentric ephemeris

LUNAR PRESENCE (OBSERVING SITE)
  The solar-presence symbol is immediately followed by a lunar-presence symbol:

        'm'  Refracted upper-limb of Moon on or above apparent horizon
        ' '  Refracted upper-limb of Moon below apparent horizon OR geocentric
             ephemeris

 R.A.___(J2000.0)___DEC. =
   J2000.0 astrometric right ascension and declination of target center.
Adjusted for light-time. Units: DEGREES

 R.A._(airls-appar)_DEC. =
   Airless apparent right ascension and declination of the target center with
respect to the Earth's true-equator and the meridian containing the Earth's
true equinox-of-date. Adjusted for light-time, the gravitational deflection of
light, stellar aberration, precession and nutation.  Units: DEGREES

 Azi_(a-appr)_Elev =
   Airless apparent azimuth and elevation of target center. Adjusted for
light-time, the gravitational deflection of light, stellar aberration,
precession and nutation. Azimuth measured North(0) -> East(90) -> South(180) ->
West(270) -> North (360). Elevation is with respect to plane perpendicular
to local zenith direction.  TOPOCENTRIC ONLY. Units: DEGREES

 CT-UT =
   Difference between uniform Coordinate Time scale and Earth-rotation
dependent Universal Time. Prior to 1962, the difference is with respect
to UT1 (CT-UT1).  For 1962 and later, the delta is with respect to UTC
(CT-UTC).  Values beyond the next July or January 1st may change if a
leap-second is introduced. Units: SECONDS


 Computations by ...
     Solar System Dynamics Group, Horizons On-Line Ephemeris System
     4800 Oak Grove Drive, Jet Propulsion Laboratory
     Pasadena, CA  91109   USA
     Information: http://ssd.jpl.nasa.gov/
     Connect    : telnet://ssd.jpl.nasa.gov:6775  (via browser)
                  telnet ssd.jpl.nasa.gov 6775    (via command-line)
     Author     : [email protected]

Studying the difference between these numbers and the ones for Skyfield, the RA and dec seem to be within an arcsecond — and while Skyfield’s unit tests insist on a far closer match to HORIZONS, in this case two different ephemerides are in play — DE-421 for the Skyfield example code, versus DE-0431 for HORIZONS — so I am going to consider sub-arcsecond agreement to be fine for the moment.

But the alt and az values are farther out. My attention is drawn to the fact that Skyfield does not yet know how to download and use high-precision Earth orientation data of the sort that HORIZONS uses as a matter of course: the delta-t offset of UT1, and the x- and y-offsets of the pole. The latter are sub-arcsecond effects:

http://maia.usno.navy.mil/bullapm_xy.png
http://www.usno.navy.mil/USNO/earth-orientation/eo-products/weekly

But the UT1 value is important for an alt-az computation because it gives Skyfield a much higher precision idea of where the Earth, and thus an observatory on the Earth, is located. I will see about getting an automated delta-t solution added to Skyfield at the sprint on Monday, but for the moment, we can bring it into agreement with HORIZONS by manually interpolating from the tables at:

http://www.usno.navy.mil/USNO/earth-orientation/eo-products/long-term

Thus:

    delta_t = (53.7882 * 20 / 31) + (53.8367 * 11 / 31)
    jd = api.JulianDate(tt=time.tt.jd, delta_t=delta_t)
    moon_obs = loc.at(jd).observe(planets['moon'])

This should bring the Skyfield results within an arcsecond of HORIZONS. Now I just need to figure out why PyEphem is so far off!

Oh — I see now! This shows the importance of using something like HORIZONS as a base reference: my first fix was not a fix at all, because I was actually making PyEphem less accurate by moving it towards Skyfield!

Now that we are giving Skyfield the crucial Earth-orientation information that it needs, its estimate is actually quite close to the original one that PyEphem was giving. So that explains why the PyEphem documentation does not constantly tell its users to supply TT instead of UTC — because, if I am reading its old C code correctly:

PyEphem treats the time you supply as TT for the purposes of Solar System positions.
PyEphem treats the time you supply as UTC for the purposes of Earth orientation.
Planets and the Moon move a little bit in 53 seconds.
But the Earth turns your observatory a lot in 53 seconds!
So it is more important for the time given to PyEphem to approximate UTC than TT.

So we should back out the change I first suggested, and return the PyEphem caller to doing:

            obs.date = time.datetime

This brings PyEphem and Skyfield into within 10 arcseconds of each other. If I am correct here about what PyEphem is doing with dates, then it is possible that sub-arcsecond agreement with HORIZONS is not going to be possible for the old PyEphem code. So for now I will let the matter rest, but if I get curious at some later date then I will share whatever else I might learn!

Maybe this can be closed, too?

@eteq

I'll let @eteq decide on this one: do we now have convergence?

@StuartLittlefair

Prompted by #10994 by @StuartLittlefair I updated the original @eteq gist https://gist.github.com/astrojuanlu/c3ab83fcdd293d60fcd95a7afe1bda5b and these are the results with recent versions of all the dependencies:

$ pip list | grep -E 'astropy|astroplan|skyfield|pyephem'
astroplan                0.7
astropy                  4.3.dev234+g2c6fd69b4 /home/juanlu/Personal/Astropy/astropy
pyephem                  3.7.7.0
pytest-astropy           0.8.0
pytest-astropy-header    0.1.2
skyfield                 1.33
$ python test_fix_gcrs.py 
With astropy v4.3.dev234+g2c6fd69b4
PyEphem (Alt, Az, Dist): (55.492015345752876 deg, 239.0318956655817 deg, 356424.8481068663 km)
Skyfield (Alt, Az, Dist): (55.49124334373045 deg, 239.03351681046558 deg, 356403.6954445954 km)
Skyfield(de430) (Alt, Az, Dist): (55.49124352738164 deg, 239.03351684329135 deg, 356403.69525128463 km)
jplephem/Astropy (Alt, Az, Dist): (56.48500486073724 deg, 238.7107927029978 deg, 357213.21025216894 km)
jplephem/Astropy(ICRS) (Alt, Az, Dist): (55.41730784517011 deg, 239.37137207735177 deg, 355352.7959745067 km)

jplephem/Astropy(ICRS) got closer to what PyEphem and Skyfield do, but I don't have enough context at the moment to judge if this issue can be closed already.

@eteq

@eteq's gist is substantially out of date. The Skyfield code should correct for polar motion (not yet available in a released version of Skyfield, but is OK as installed from GitHub. The astropy code should be updated to:

return astropy.coordinates.get_body('moon', time, location, 'de430').transform_to(AltAz(obstime=time, location=location))

Below is a comparison of the output of Skyfield master and Astropy master for the Moon:

Indeed, with Skyfield @ master and the updated Astropy code, these are the results of the gist:

$ python test_fix_gcrs.py 
With astropy v4.3.dev234+g2c6fd69b4
PyEphem (Alt, Az, Dist): (55.492015345752876 deg, 239.0318956655817 deg, 356424.8481068663 km)
Skyfield (Alt, Az, Dist): (55.49124334373045 deg, 239.03351681046558 deg, 356403.6954445954 km)
Skyfield(de430) (Alt, Az, Dist): (55.49124352738164 deg, 239.03351684329135 deg, 356403.69525128463 km)
jplephem/Astropy (Alt, Az, Dist): (55.49118251346562 deg, 239.03406534146313 deg, 356399.58247342776 km)

Let's close this one then!

Navigation Menu

Repository navigation

Determine why the Moon's AltAz coordinates vary across codes #4269

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Determine why the Moon's AltAz coordinates vary across codes #4269

Description

Activity

bmorris3 commented on Oct 26, 2015

eteq commented on Oct 30, 2015

brandon-rhodes commented on Nov 2, 2015

bmorris3 commented on Nov 2, 2015

brandon-rhodes commented on Nov 2, 2015

brandon-rhodes commented on Nov 2, 2015

brandon-rhodes commented on Nov 2, 2015

bsipocz commented on May 14, 2016

mhvk commented on May 14, 2016

astrojuanlu commented on Nov 23, 2020

StuartLittlefair commented on Nov 23, 2020

astrojuanlu commented on Nov 23, 2020

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions