Improve get_body by using location of observer in light travel time correction - #10292

Just letting you know. The Sun isn't a very good test for this. It's much better to use the Moon. I'm currently seeing +/-200 mas error in RA/Dec and +/-500 meters in distance over the course of a day using an older version of astropy.

@mkbrewer - you're right, it would probably be good to have a test for the Moon. Might you be able to make a PR?

Independently, thanks again for pointing out these problems!

Since I haven't looked into your test procedures, I'm not sure what you would want. What I did was simply compare the Moon's position in RA/DEC, distance and Alt/Az in one minute intervals for 24 hours with my own solar system ephemeris program based on SOFA. I tested the output of my program for this period on the Moon against both Skyfield and JPL Horizons and got agreement to 0.5 mas in Declination, 1.5 mas in Right Ascension and 0.7 meters in distance against Horizons and similar results against Skyfield so long as I disabled polar motion. Skyfield seems to have a problem with that, but it's only about 12 mas in Declination and 4 meters in distance. As does SOFA, btw, but it only affects Alt/Az. That error is a rather huge at +/-1 arcsecond! (It affects astropy as well. If you are interested, please see the issue I opened on the Skyfield repository). Horizons also has an issue with apparent Right Ascension, being 53 mas off from SOFA. I pointed this out to them last June and, of course, they pointed their fingers at SOFA. At least they are now acknowledging the issue in their print out though. Of course, these errors are all very small compared to the errors I got testing the older version of astropy.

I could certainly provide you with my test data and the script I used to generate the astropy data to compare against it if that would be helpful.

@mkbrewer that would be fantastic, thank you!

Perhaps open a new issue describing the Moon problems and provide the data there?

astropy_tests.zip

I was just informing you of the issues you might have comparing astropy Moon data to Skyfield or Horizons. There may be a problem with polar motion in SOFA that affects astropy, but I'm still looking into that. I will certainly raise an issue on that if I determine that there is, in fact, a problem, but I'm not there yet.

The file I attached contains a script for printing topocentric astrometric, apparent postions wrt the true equator and equinox, distance and alt/az that can be directly compared to equivalent Skyfield and Horizons output if you so choose. I also included the output from my program that you can compare to output from the script for astropy both before and after the modification in this merge to see how much of an improvement resulted from the change you made. I used the same EOPs (dut1,xp,yp) that astropy used on the test date in the script for the input to my program, so my data should be completely compatible with astropy in that respect.

@StuartLittlefair please let me know if this works for you.

Thanks @mkbrewer. This helped me a lot dig into what is happening. The short story is that with this PR, the astropy positions for the Moon are pretty damn good! For the long story, see below.

I used @mkbrewer's scripts together with my additions in this PR to get the following:

from astropy.time import Time
from astropy.coordinates import solar_system_ephemeris, EarthLocation, get_body, get_body_barycentric, SkyCoord
from astropy.coordinates.solar_system import _apparent_position_in_true_coordinates

# For comparison with JPL Horizons as this is what they use
solar_system_ephemeris.set('de430')

class_lat = -22.959748*u.deg
class_lon = -67.787260*u.deg
class_elev = 5186*u.m
class_loc = EarthLocation.from_geodetic(class_lon, class_lat, class_elev)

t = Time('2017-04-06 00:00', scale='utc')
moon_barycentric = _get_apparent_body_position('moon', t, 'de430', 
                                    class_loc.get_gcrs_posvel(t)[0]) 
obs = get_body_barycentric('earth', t, 'de430') + class_loc.get_gcrs_posvel(t)[0]

# ASTROMETRIC ICRF
SkyCoord(coord.GCRS(moon_barycentric - obs, obstime=t)).to_string(style='hmsdms')
'09h05m23.649s +15d36m28.9186s'

c.f JPL Horizons "09 05 23.65 +15 36 28.9" and SkyField "09 05 23.64 +15 36 28.5"

What about JPL's True Equator, True Equinox of date ("apparent") coordinates? Let's look at what the result of calling astropy's get_body is with the changes in this PR:

In [165]: _apparent_position_in_true_coordinates(get_body('moon', t, location=class_loc)).to_string(style='hmsdms')                                                
Out[165]: '09h06m21.1377s +15d32m12.2391s'

c.f JPL "09 06 21.13 +15 32 12.2" and SkyField "09 05 24.31 +15 36 25.3".

You can see there is now excellent agreement with astropy and JPL Horizons for the moon! In fact, we're quite a bit more accurate than Skyfield for this one example.

This also gives us confidence that implementing the transform between GCRS and TETE requested in #10230 will make it easier for people to use JPL Horizon query results within astropy and keep precision.

Re: #10292 (comment)

@StuartLittlefair , so are we expecting a follow-up PR to this? Or is this like a astropy-tutorial material?

Very good news! And particularly nice that it works so well with what really is a small bug fix! Would be nice to have the Moon example as a test case and/or as an example for get_body in the documentation.

I think what I will do is work on #10230 when I get the time. Including lunar test cases and a tutorial about the different frames, how they relate and when to use each might fit well in a PR for that issue.

Here are the comparisons I'm seeing from the coordinates that you posted with my own program based on SOFA, JPL Horizons and Skyfield.

                       RA (deg)      Delta RA (mas)     Dec (deg)  Delta Dec(mas)
mkb       astrometric 136.348537482     0.0648        15.608032955    -0.0380
          apparent    136.588073175     2.0700        15.536733269    -0.6684

jpl       astrometric 136.348537637    -0.4932        15.608032907     0.1348
           apparent   136.588058954    53.2656        15.536733166    -0.2976

skyfield  astrometric 136.348537209     1.0476        15.608034242    -4.6712
             apparent 136.588072886     3.1105        15.536734528    -5.2008

For my program, the astrometric positions are good to the uas level. The discrepancy in apparent positions is a bit perplexing. I am using the EOPs used by astropy 3.2.2 here. Could they be slightly different?

The JPL results are slightly different due to the fact that they have their own set of EOPs, plus the 53 mas offset in apparent RA. Hint: When comparing with Horizons, it's best to choose decimal degrees for the output with "extra precision" checked.

The Skyfield results do not include polar motion, so the discrepancy there is not surprising.

@StuartLittlefair could you please run the script In the astropy_tests.zip file that I posted above and upload the output as a .txt file? I'd like to compare your new positions over a 24 hour period that are generated by that script.

@pllim I was responding to @mhvk's request for a PR for testing on the Moon.

@mkbrewer - you're right, it would probably be good to have a test for the Moon. Might you be able to make a PR?

Independently, thanks again for pointing out these problems!

Since I don't know how to submit a PR for tests, that file astropy_tests.zip was my response. It contains a script to output Moon data from astropy over a 24 hour period and output from my own program based on SOFA as comparison data..

@mkbrewer , to be credited properly for your contribution, it is still best for you to submit the PR yourself. Perhaps you might find https://docs.astropy.org/en/latest/development/workflow/development_workflow.html useful. Thanks!

One possibility for the discrepancy in apparent position between astropy and my program is the issue with gravitational light deflection #10291 which I see is still open.

I doubt it. For the example I looked at it made a negligible difference.

Is there a reason you don’t want to upgrade astropy btw? Not saying you should but there may other bug fixes you are missing.

I’ll post your example run on astropy dev as soon as I can. Lots to do and it’s early evening my time.

Hi @mkbrewer

I've attached the full output from your code to this post. Here's a snippet. First the file of positions from your SOFA program that you attached above:

UT,LST,ICRS_RA,ICRS_DEC,RA,DEC,AZ,EL,DIST,DUT1,XP,YP 
2017-04-06 00:00:00.000000,8.445790414,136.348537301,15.608033005,136.588072995,15.536733318,15.032679644,50.303109080,376252883.251,0.462380,0.009587,0.384603
2017-04-06 00:01:00.000000,8.462502712,136.354077845,15.606277514,136.593609981,15.534971532,14.674616403,50.362682839,376251255.709,0.462379,0.009588,0.384604

And what we get on astropy master:

UTC,LST,ICRS_RA,ICRS_DEC,RA,DEC,AZ,EL,DIST,DUT1,XP,YP
2017-04-06 00:00,8.445791832,136.348537485,15.608032954,136.588073934,15.536733085,15.032729808,50.303166721,376252883.195,0.463626,0.009587,0.384603
2017-04-06 00:01,8.462504130,136.354078025,15.606277464,136.593610917,15.534971300,14.674666644,50.362740603,376251255.697,0.463625,0.009588,0.384604

Note that amongst the differences between my version and yours are discrepancies in LST and UTC-UT1.

Here's the output from the JPL web form for comparison:

*******************************************************************************************************************************************
 Date__(UT)__HR:MN     R.A._______(ICRF)_______DEC R.A.__(airless-appar)__DEC. Azi_(a-appr)_Elev            delta      deldot 1-way_down_LT
*******************************************************************************************************************************************
$$SOE
 2017-Apr-06 00:00  m  136.348537637  15.608032907 136.588058954  15.536733166  15.0327  50.3030 3.7625288339E+05  -0.0279606    0.02091741
 2017-Apr-06 01:00  m  136.678112032  15.494425098 136.917408421  15.422751330 352.5663  51.3299 3.7633487126E+05   0.0737328    0.02092197
$$EOE

The astropy coordinates I get are 0.5 mas (astrometric) and 52 mas (apparent) different from JPL.
moon_topo_astropy_master.txt

Thank You. That is very helpful. Here my results:

The average difference in DUT1 is 250 ns, so that is negligible. The difference in LST is due to the polar motion correction to longitude from equation (6.20) in USNO Circular 179. I am deriving my LST from equation (6.19) there. The hour angles are supposed to be the same and SOFA doesn't appear to add that correction to the local earth rotation angle, so I added it to LST. This is still under investigation on my part, so I'm not sure it is correct yet. In any case, it has no effect on anything else as I'm not using LST in any of my calculations.

The discrepancy in apparent position is quite evident. @StuartLittlefair, you say that you looked at an example after correcting the light deflection problem then? By "negligible", do you mean much less than the 2.5 mas I'm seeing? If so, this could be due to the less than rigorous way that SOFA handles diurnal aberration when transforming from CIRS to observed. I am still investigating this and will post further comments about that under issue #10230.

There is a bit of noise in your distance calculation which could very well be due to the 1.0e-8 seconds (3 meters at the speed of light) that you are iterating down to.

There is a fairly large discrepancy in Alt/Az, I believe that this is due to some inconsistency regarding polar motion, but I'm still investigating that. This data will be very helpful in that regard.

As you can see, this test against the Moon is very good for finding small problems. I will look into converting it into a pass/fail test that I could submit as a PR if I have the time. I'm more interested in actually solving the problems I see, rather than getting credit for a PR containing this test though, so it won't be a high priority for me. In the meantime, if one of you experts want to go ahead and implement this as a test, please feel free to use my script and data to do so.

On second thought, the discrepancy in apparent position is unlikely to be due to diurnal aberration as it's not diurnal. Instead there is gradient that appears to be related to the Moon's position wrt the Earth. This makes light deflection a better candidate for this discrepancy.

Yes! The discrepancy is due to light deflection. I reverted to using the incorrect deflection that SOFA uses in Ldsun() in my program and the discrepancy is reduced by 90%.

Thanks for the sleuthing! And thanks for cross-posting at #10291 - indeed, let's move further discussion there since people will not easily find it here.

@StuartLittlefair — Could you share the Skyfield code that you used to generate the coordinates that do not quite match the other libraries?

I just wrote a script to see the results myself and try fixing Skyfield:

from skyfield.api import Topos, load
ts = load.timescale(builtin=True)
t = ts.utc(2017, 4, 6)
planets = load('de430t.bsp')
topos = Topos('22.959748 S', '67.787260 W', elevation_m=5186.0)
pos = (planets['earth'] + topos).at(t).observe(planets['moon'])
ra, dec, distance = pos.radec()
print('Astrometric  Right ascension {}  Declination {}'.format(ra, dec))
ra, dec, distance = pos.apparent().radec(epoch='date')
print('Apparent     Right ascension {}  Declination {}'.format(ra, dec))

But, alas, the output — at least here on my laptop — agrees perfectly with the HORIZONS output and Astropy output that you quote above:

Astrometric  Right ascension 09h 05m 23.65s  Declination +15deg 36' 28.9"
Apparent     Right ascension 09h 06m 21.14s  Declination +15deg 32' 12.2"

To be more exact: Skyfield here agrees most closely with AstroPy, while we both dissent very slightly from the HORIZONS result with respect to the arcseconds of the apparent RA, where HORIZONS has 21.13 but both Skyfield and AstroPy prefer 21.14.

If you’ll share your code, I'll be happy to improve Skyfield's docs to better help programmers avoid whatever snag your code ran into that hurt its accuracy. Thanks!

@brandon-rhodes do you have an alert for mentions of Skyfield ;) Thanks very much for following up...

You code matches mine exactly apart from the following differences:

1. I did not specify builtin=True for the timescale, but I suspect this can't hurt.
2. I used the de421 ephemeris (didn't know about de430t and de430 did not work).
3. I did not specify epoch='date' for the apparent position.

I suspect the last one explains most of the difference, since the astrometric coordinates are very close. The use of radec is very clearly documented in SkyField, but I was rushing and didn't read the docs! My apologies if I've caused you work as a result...

@brandon-rhodes do you have an alert for mentions of Skyfield ;)

Alas, I do not! That’s not a bad idea, though. In this case it was an issue opened by @mkbrewer that mentioned a similar effort to reach agreement between his library and Astropy that led me to cross-reference with his issues over here.

I’m glad you found the docs clear! All I will do, then, is try to make sure that my example code is rigorous in using of-date RA and dec every time it calls apparent(), so that whichever example of that call a programmer’s eyes might fall on is likely to include the pattern. Thanks for your quick response!

As you’ll see in the commit just mentioned, upon review I think some blame might be placed on Skyfield’s radec() docs.

Instead of showing how to use radec() with real coordinates, the example — in a misguided effort to make the example “simple” and save CPU — instantiated a meaningless ICRS([1, 2, 3]) vector and then asked its RA and declination. First, as I was reminded recently when responding to another user question, users are far better served by seeing examples where methods are called in realistic situations; many users never instantiate an ICRS manually, and so are confused about when they would be using one. Second, if the coordinate is a fake one like [1, 2, 3], then there's no reason to prefer ICRS or epoch-of-date, so the example can’t illustrate the difference.

So instead I’ve pivoted the radec() docstring to a full realistic example of asking Skyfield for coordinates, and per best practices generate ICRS for my astrometric but epoch-of-date for my apparent coordinates. Hopefully that will make it a bit more likely that user code makes the same choice except in cases where the user really deliberately wants the idiosyncratic combination of apparent with ICRS.