Repository navigation
Determine why the Moon's AltAz coordinates vary across codes #4269
Description
Activity
@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.datetimeshould in fact be:
obs.date = time.tt.datetimeThis 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!
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, noelevation_mkeyword 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.csaysPrecision < 0.05" from -1369 to +2950and if by"they mean arcseconds then, nope, that is not the source of the difference we are seeing here!
- When the Skyfield
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/weeklyBut 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.datetimeThis 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?
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'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 @
masterand 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!
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:
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