Solar system and lunar ephemerides via jplephem - #4746

@bmorris3 - great! Even if it does mean one GSOC project is now partially done...

One general comment: Right now, you hard code the specific ephemeris file that should be used; I think it would be best to leave an option for selecting it. Indeed, even better if we could make it a sufficiently abstracted interface that one could also use non-JPL ephemerides. Specifically, it turns out that within erfa there are pretty good approximations to the planet positions, which one can access via plan9 (see #4736, where it is used to calculate the position of the moon). From a quick look, it seems one could fairly easily substitute these.

Conversely, in the new Time.light_travel_time method (#4314), erfa is used to calculate an approximate position of the Earth; clearly, that could use an extension to use the precise one (but this is trivially added later).

I am glad that jplephem might prove useful to AstroPy!

When I see a difference with HORIZONS that is small for the distant planets but large for the Moon, it usually involves time, since a small error in time causes a bigger problem for the fast-moving Moon than for the slower planets. I quickly go check things like whether my leap-second table is up to date and (if Earth rotation is involved) my Delta-T estimate. Try combing through the HORIZONS output, looking at every variable and model in use for time and for the bodies involved — any difference could be the source of the difference.

I have been able to get Skyfield to 0.001 arcsecond agreement with the Naval Observatory, and am now honing in on HORIZONS myself. I have been thinking of drawing up a list of errors that can happen, and maybe even writing a little routine that, given a date and position for something like the Moon, goes through all combinations of time scales and frames of reference, trying to determine which switches need to be flipped how to generate the number that a user is seeing from a web site.

It's a huge issue — users always want to check astronomy libraries against web sites (often not even good ones like HORIZONS), and trying to explain why an astronomy library is giving an answer that's an arcminute off is difficult when the web site is a black box and you can't see its methodology. I've often wished that the HORIZONS source code were available.

@bmorris3 to try and answer some of your questions:

• GCRS coordinates without an EarthLocation currently default to "true" GCRS coordinates - i.e Geocentric. I think this is a sensible default in your case as well.
• I don't think inconsistencies between the geocentre are responsible for the differences between your code and JPL HORIZONS. As written the GCRS frame you are making will transparently pass longitude (ra) and latitude (dec) to the SkyCoord; any discrepancies are already present in the earth_to_planet_vector.

As a comment on your PR, get_absolute_planet_position currently implements an ICRS->GCRS transformation by hand. It might be an idea to simply create an ICRS object directly from the barycentric position returned by jplephem for the specified body and use the coordinate transformation machinery from then on. This at least has the merit of making sure that all coordinate transforms are kept in one place.

I also think as written your coordinates are not properly taking into account the location of the observatory; merely setting the obsgeoloc attribute of the GCRS coordinate doesn't actually give a different position. Creating an ICRS object and transforming to a GCRS coordinate with a specified obsgeoloc would address this.

Also - I think you should be supplying times in TDB to the jplephem kernel.

The tests will fail in the ci chain because I've put jplephem in the optional dependencies, but is not currently accessible via conda (you can grab it with pip install jplephem). I'm not sure if this is a showstopper for us – would it be possible to push a copy of jplephem to conda?

@bmorris3 - Put in under PIP_DEPENDENCIES (but only for the builds where the other optionals are tested, eg where CONDA_DEPENDENCIES=$CONDA_ALL_DEPENDENCIES. It's perfectly fine to have optional dependencies that are not (yet) available on a conda channel, and install them via pip.

On the other hand, I see that there is a conda recipe in the jplephem repo. Would you @brandon-rhodes consider it adding to the openastronomy conda channel? https://github.com/OpenAstronomy/conda-channel

Thanks @bsipocz, I hadn't noticed that line in the .travis.yml file. Maybe that last commit will produce more interesting test results.

@bsipocz What are the steps for adding it to the conda repository? I've never put a package in one before! :)

@brandon-rhodes - Most probably adding and entry to https://github.com/OpenAstronomy/conda-channel/blob/master/requirements.yml will do the job. If it doesn't work with the repo's default recipe then you can also add the one from jplephem.
I only added packages that worked out of the box, thus never dealt with custom recipes, but I'm sure @mwcraig will know the solution for most problems.

@bmorris3 - what is the status of this? Is there any chance we can finish it up by Friday?

@astrofrog - agree it would be great to get this in, so let me try to summarize what I think are the hurdles. The most important one -- which should really be resolved before we merge -- is the question whether the positions match properly. It may be that much of this is in the transformation from ICRS coordinates relative to the solar system barycentre to coordinates relative to some observer on Earth. @StuartLittlefair suggested using the coordinate machinery for this transformation, rather than do it by hand, but that means one starts using approximate positions for the Earth, when one of the beauties of jplephem is that it gives us much more accurate ones.

My suggestion would be to proceed by first writing tests for the absolute (ICRS) positions -- as I noted above, these are very useful in themselves, so it is also good to have a public method for those (maybe just a different value of location; e.g., if we allow strings which by default get passed to EarthLocation.of_site, which could include items for barycentre and geocentre).

Apart from the above more serious hurdle, I think the main modification to be made is to make strings the default for planet names (I think it is safest not to allow integers at all for now).

Another smaller question was whether to allow a choice of ephemeris used; this could either be another DE... file, or plan9 in erfa. But I'd guess this could be done later too.

@astrofrog – I won't be able to work on this again this week.

If someone wants to push this forward, a place to start would be to review @StuartLittlefair's PR to my branch bmorris3#1, which does some of the more rigorous comparisons that @mhvk, @brandon-rhodes and others have encouraged.

In the pull request @bmorris3 mentions above I've included the ability to request planet positions by index or string. I've also partially resolved the discrepancies with JPL Horizons.

Requested positions are in agreement(better than 1 arcsec) for geocentric observers, but there's a problem for observers with an EarthLocation specified that I think is to do with the GCRS coordinate transform machinery, but which I have not managed to pin down yet. The differences can be up to 1 degree for the Moon, but is less than an arc minute for other objects.

Arrggh, should have commented here rather than in @bmorris3's fork:

@StuartLittlefair - I looked at your improvements and generally they look good! Since @bmorris3 doesn't have time, might it be an idea to open a new PR with his and your commits combined, so one can comment there?

One question on principle: I noticed that you correct the time for light travel time. However, you calculate positions of both the body and the earth at this slightly earlier time; I think you should change the time only for the body in question.

On the final conversion to apparent position for a non-geocentric location: I wonder if the GCRS conversion doesn't somewhere implicitly assume that distances are large..

Finally, I think the routine just getting the barycentric position of the body should be public, as it is very useful.

Great suggestions as always @mhvk.

If I made a single, new, pull request, how would I make sure that @bmorris got proper credit?

Just don't squash all the commits! You presumably already have his branch locally, so what you would need to do is to rebase it on master (maybe git rebase -i upstream/master so you can squash e.g. the last 3 commits here into one, but keep yours separate). And then push it back to your branch on github and do a pull request from there.

As requested by @mhvk I have created PR #4890 to include @bmorris3 's original code and my modifications therein. That PR should be considered as superseding this one.

Closing this since development moved to #4890.