Allow use of JPL ephemerides in coordinate transforms and light travel time corrections #5273

StuartLittlefair · 2016-08-28T00:51:08Z

This PR is based upon the changes in #5231, and perhaps should include changes in #5253 and #5254 as well; I will rebase once #5231 is merged and also if #5253 and #5254 are merged before this.

This PR fixes #4942 by allowing the coordinate transformations, and the light travel time corrections which use them to respect the state of the solar_system_ephemeris context manager. light_travel_time takes an explicit ephemeris argument whilst the co-ordinate transformations don't. To set the ephemeris for coordinate transformations, one uses solar_system_ephemeris.set(). The reason for this decision was that adding an ephemeris argument to the co-ordinate transforms meant adding a kwargs argument to all frame transforms when this was not really needed for most of them.

With this PR, one can use:

from astropy.time import Time
from astropy.coordinates import GCRS, solar_system_ephemeris, SkyCoord, EarthLocation

ippeg = SkyCoord.from_name('IP Peg')
now = Time.now()
greenwich = EarthLocation.of_site('greenwich')
now.location = greenwich
gcrs_frame = GCRS(obstime=now)

dt = now.light_travel_time(ippeg)  # use builtin ephemeris
dt = now.light_travel_time(ippeg, ephemeris='jpl') # use jpl ephemeris
new_coo = ippeg.transform_to(gcrs_frame) # use builtin

solar_system_ephemeris.set('jpl')
dt = now.light_travel_time(ippeg)  # use jpl ephemeris
new_coo = ippeg.transform_to(gcrs_frame) # use jpl ephemeris

StuartLittlefair · 2016-08-28T00:57:02Z

It's worth commenting on performance changes; this makes marginal difference to the performance of the transformations using the builtin ephemeris. There's a circa 5 percent slowdown due to additional overhead in preparing astropy objects for erfa calls.

The JPL ephemerides are roughly 2x slower than the builtin equivalents.

StuartLittlefair · 2016-09-30T11:47:25Z

This is now rebased following the changes in #5231. @mhvk perhaps you could take a look?

mhvk · 2016-10-03T03:53:25Z

Now at the right PR: @StuartLittlefair - I've been travelling quite a bit, so may be a bit slow, but remain very interested! If I still don't get to this by next week, do feel free to ping again.

StuartLittlefair · 2016-10-17T11:05:33Z

Pinging @mhvk as requested...

mhvk · 2016-10-17T18:36:02Z

astropy/coordinates/builtin_frames/ecliptic_transforms.py

-    # first offset to heliocentric
-    heliocart = (from_coo.cartesian.xyz).T + delta_bary_to_helio * u.au
+    # offset to heliocentric
+    heliocart = CartesianRepresentation(from_coo.cartesian.x + bary_sun_pos.x, from_coo.cartesian.y + bary_sun_pos.y,


move "y" to its own line so the first line is shorter and it is easier to see what happens
(of course, with #5301, this would just be from_coo.cartesian + bary_sun_pos...)

good plug :)

mhvk · 2016-10-17T18:38:19Z

astropy/coordinates/builtin_frames/ecliptic_transforms.py

+    # get barycentric sun coordinate
+    # this goes here to avoid circular import errors
+    from ..solar_system import get_body_barycentric, solar_system_ephemeris
+    ephemeris = solar_system_ephemeris.get()


Getting the ephemeris is not needed; the get_body_barycentric function does this by default.

mhvk · 2016-10-17T18:38:34Z

astropy/coordinates/builtin_frames/ecliptic_transforms.py

+    from ..solar_system import get_body_barycentric, solar_system_ephemeris
+
+    # get barycentric sun coordinate
+    ephemeris = solar_system_ephemeris.get()


As above: this should not be needed.

mhvk · 2016-10-17T18:38:50Z

astropy/coordinates/builtin_frames/ecliptic_transforms.py

+    ephemeris = solar_system_ephemeris.get()
+    bary_sun_pos = get_body_barycentric('sun', from_coo.equinox, ephemeris=ephemeris)
+
+    newrepr = CartesianRepresentation(intermed_repr.x - bary_sun_pos.x, intermed_repr.y - bary_sun_pos.y,


As above: split "y" to its own line

mhvk · 2016-10-17T18:44:16Z

astropy/coordinates/builtin_frames/utils.py

+    from ..solar_system import (get_body_barycentric, get_body_barycentric_posvel,
+                                solar_system_ephemeris)
+    # get barycentric position and velocity of earth
+    ephemeris = solar_system_ephemeris.get()


Again, if you're just using the default ephemeris, this is not necessary.

mhvk · 2016-10-17T19:13:25Z

astropy/time/core.py

+                                   HCRS, ICRS, GCRS, solar_system_ephemeris)
+
+        # get the ephemeris for later (we don't want to set it to None)
+        if ephemeris is None:


This slightly begs the question whether we're not doing the wrong thing for solar_system_ephemeris.set(). Looking at the ScienceState code, we could None stand for "don't change anything" by using

@classmethod def validate(cls, value): if value is None: return cls._value # Set up Kernel; if the file is not in cache, this will download it. cls.get_kernel(value) return value

mhvk · 2016-10-17T19:15:22Z

astropy/time/tests/test_corrs.py

+
+    @remote_data
+    @pytest.mark.skipif('not HAS_JPLEPHEM')
+    def test_ephemerides(self):


Again, ideally we'd check against a known-to-be-correct number...

As above. This may mean finally biting the bullet and installing and figuring out TEMPO2.

mhvk · 2016-10-17T19:17:55Z

docs/coordinates/transforming.rst

+For locations at large distances from the Solar system, using the JPL ephemerides will make a
+negligible difference, on the order of micro-arcseconds. For nearby objects, such as the Moon,
+the difference can be of the order of milli-arcseconds. For more details about what ephemerides
+are available, see :ref:`astropy-coordinates-solarsystem`.


I would write here "..., including the requirements for using JPL ephemerides, see ..."

With that, you could remove the note below (which I feel is better kept in just one place).

mhvk · 2016-10-17T19:18:35Z

docs/time/index.rst

+The difference between the builtin ephemerides and the JPL ephemerides is normally
+of the order of 1/100th of a millisecond, so the builtin ephemerides should be suitable
+for most purposes. For more details about what ephemerides are available,
+see :ref:`astropy-coordinates-solarsystem`.


Same as above.

mhvk · 2016-10-17T19:20:09Z

docs/time/index.rst

+
+    >>> ltt_bary_jpl = times.light_travel_time(ip_peg, ephemeris='jpl') # doctest: +REMOTE_DATA +IGNORE_OUTPUT
+
+The difference between the builtin ephemerides and the JPL ephemerides is normally


Would it be an idea to actually show the various values that have been gotten here and above?

to the precision that's output by default, they look identical to the barycentric corrections with the default ephemeris, so I wonder if it will just look confusing.

Hmm, I see your point. I guess in principle that is interesting on its own, and you could make the fact that the difference is small explicit by displaying

(ltt_bary_jpl - ltt_bary).to(u.us)

(simularly, in the other cases, it may be useful to display ...to(u.min) after just once showing the normal result (and thus that normally one gets a TimeDelta out).

Obviously, none of this is very important, but I think it will help users get a sense of when they would need to use the jpl ephemeris.

mhvk · 2016-10-17T19:24:31Z

@StuartLittlefair - now had a good look and like this very much. My comments are mostly cosmetical, so should take little effort to address.

@eteq, @taldcroft, @astrofrog - with this PR, solar_system_ephemeris will start affecting coordinate transformations. I think this is proper but thought you should also consider whether this is indeed the right way forward.

mhvk · 2016-10-18T14:10:02Z

docs/time/index.rst

@@ -901,10 +901,14 @@ calculate light travel times to the barycentre as follows::
    >>> times = time.Time([56325.95833333, 56325.978254], format='mjd',
    ...                   scale='utc', location=greenwich)
    >>> ltt_bary = times.light_travel_time(ip_peg)
+    >>> ltt_bary


Here and in the next output, add # doctest: +FLOAT_CMP, just so we don't get hit by rounding errors.

mhvk · 2016-10-18T14:10:33Z

Looks good! A few small comments remaining...

mhvk · 2016-10-19T13:40:10Z

astropy/coordinates/builtin_frames/ecliptic_transforms.py

-    # first offset to heliocentric
-    heliocart = (from_coo.cartesian.xyz).T + delta_bary_to_helio * u.au
+    # offset to heliocentric
+    heliocart = CartesianRepresentation(from_coo.cartesian.x + bary_sun_pos.x,


Hmm, now I wonder whether to optimize when really we just want to have #5301 in... The way this is written, cartesian gets called three times, which each time will calculate all of x, y, and z, so this can be sped up nearly a factor of three by doing

from_coo_cartesian = from_coo.cartesian heliocart = ... from_coo_cartesian.x + ...

mhvk · 2016-10-19T13:41:38Z

astropy/coordinates/builtin_frames/icrs_cirs_transforms.py

-    astrom = erfa.apcs13(jd1, jd2, pv)
+    earth_pv, earth_heliocentric = prepare_earth_position_vel(gcrs_frame.obstime)
+
+    # get astrometry context


This left-over comment should be deleted, I think.

The "get astrometry context"? If so, no, that's a reference to the astrom object. I've made the comments in this section clearer.

mhvk · 2016-10-19T13:45:34Z

docs/time/index.rst

+
+    >>> ltt_bary_jpl = times.light_travel_time(ip_peg, ephemeris='jpl') # doctest: +REMOTE_DATA +IGNORE_OUTPUT
+    >>> ltt_bary_jpl # doctest: +REMOTE_DATA +FLOAT_CMP
+    <TimeDelta object: scale='tdb' format='jd' value=[-0.0037715  -0.00377286]> # doctest: +REMOTE_DATA +FLOAT_CMP


You don't need to put the # doctest on the output lines (also the line below).

mhvk · 2016-10-19T13:46:48Z

@StuartLittlefair - yet more comments, all trivial.

@taldcroft, @eteq, @astrofrog - unless I hear that you would like to have a look, I'll merge this tomorrow.

mhvk · 2016-10-20T20:54:22Z

This all looks OK to me, so I'll merge. Thanks very much, @StuartLittlefair!!

astrobot · 2016-10-20T20:54:29Z

@mhvk - thanks for merging this! However, I noticed the following issue with this pull request:

Changelog entry not present (or pull request number missing) and neither the Affects-dev nor the no-changelog-entry-needed label are set

Would it be possible to fix this? Thanks!

This is an experimental bot being written by @astrofrog - let me know if the message above is incorrect!

eteq · 2016-10-26T19:06:50Z

@astrofrog - any idea why the bot failed here? there is a changelog entry

pllim · 2016-10-26T19:16:30Z

No PR number attached to the change log entry.

MSeifert04 · 2016-10-26T19:24:29Z

Just out of interest: When should a PR be squashed?

I ask this here because I've noticed that it contains 28 commits with ~300 line changes, several of them about formatting, fixing tests, etc. On the other hand I've seen PRs (today one of mine, #5214 with 8 commits and ~400 line changes) with fewer commits and/or more changed lines where someone requested to squash them.

Not that I mind, but it would help me decide when to ask someone to squash their PR (or not).

As an afterthought: Should I have asked this on astropy-dev?

eteq · 2016-10-26T20:02:48Z

@mhvk @StuartLittlefair - this all looks OK to me except for a few minor tweaks - I went ahead and made a separate PR on this, though (#5436), given that I never found the time to review this when it was actually submitted...

eteq · 2016-10-26T20:05:00Z

@MSeifert04 - that's a matter of opinion. Personally I never ask for squashing unless there's a substantial file size problem (i.e. adding and then removing of a large data file). I find it useful at times to see the logic in a series of commits. But I know others prefer more coherent commits. So far we haven't had a clear policy other than "whatever the primary reviewer thinks"...

mhvk · 2016-10-26T21:06:00Z

@MSeifert04 - just to add to @eteq's comment: I think for astropy we normally don't ask for squashing (though I also do if the implementation changed completely), but numpy does always ask it. They're also different on insisting on more than a one-liner for the commit, as they don't want to rely on the github history to be around forever (probably smart, but then, they did move system, so know what it involves...).

pllim · 2016-10-27T14:32:09Z

I am guilty for requesting squashing. I like neat commit history. 😅

Minor adjustments to #5273

pllim added the coordinates label Aug 29, 2016

StuartLittlefair added 12 commits September 30, 2016 10:00

transforms now use current solar_system_ephemeris

2961e46

light_travel_time gains an ephemeris argument

22c6795

add docstring to get_cip

5049b4e

raise exception on ICRS-HCRS trans with no dist

c189aaa

move prepare_earth_position_vel to utilities

a08a76f

ecliptic transforms use solar_system_ephemeris

04de028

remove unnecessary rollaxis use

7e63ec5

added documentation of changes

e77ba84

added tests for ephemerides

8669dd8

fix bugs in docs

fc3fd3d

bug in changelog

3919f6d

accomodate name changes in PR#5231

de00e2e

mhvk added time Enhancement labels Oct 17, 2016

mhvk added this to the v1.3.0 milestone Oct 17, 2016

mhvk assigned eteq and mhvk and unassigned eteq Oct 17, 2016

mhvk reviewed Oct 17, 2016

View reviewed changes

StuartLittlefair added 5 commits October 17, 2016 22:07

remove redundant accessing of default ephemeris

cfc2443

neater formatting

645c2a5

move prepare_earth_position_vel

fb8ee08

drop empty line and document rollaxis operations

f9cb4f8

make solar_system_ephemeris do nada if set to None

ebbb1e8

mhvk reviewed Oct 18, 2016

View reviewed changes

StuartLittlefair added 6 commits October 18, 2016 21:58

more redundant ephemeris calls removed

feea830

adding FLOAT_CMP to doctests

51feb43

more FLOAT_CMPs

249e71d

make sure we have a regular array for heliocentre

ed8d604

removed duplicate erfa import

51cf310

got rid of rollaxis in favour of newaxis broadcast

a439f7d

mhvk reviewed Oct 19, 2016

View reviewed changes

StuartLittlefair added 3 commits October 19, 2016 22:12

optimisation of cartesian math

0cd4303

clarifying astrometry setup

768e3a1

remove doctest on output

1a20b34

mhvk merged commit 1809afc into astropy:master Oct 20, 2016

mhvk added the Affects-release label Oct 20, 2016

eteq added a commit to eteq/astropy that referenced this pull request Oct 26, 2016

fixes to changelog from astropy#5273

84a5eb8

eteq mentioned this pull request Oct 26, 2016

Tweaks from 5273 #5436

Merged

eteq added a commit that referenced this pull request Oct 27, 2016

Merge pull request #5436 from eteq/tweaks-from-5273

3256f45

Minor adjustments to #5273

StuartLittlefair mentioned this pull request May 15, 2017

WIP: improve speed of coordinate transformations for arrays of obstime #6068

Closed


		>>> ltt_bary_jpl = times.light_travel_time(ip_peg, ephemeris='jpl') # doctest: +REMOTE_DATA +IGNORE_OUTPUT

		The difference between the builtin ephemerides and the JPL ephemerides is normally

Allow use of JPL ephemerides in coordinate transforms and light travel time corrections #5273

Allow use of JPL ephemerides in coordinate transforms and light travel time corrections #5273

Conversation

StuartLittlefair commented Aug 28, 2016

StuartLittlefair commented Aug 28, 2016

StuartLittlefair commented Sep 30, 2016

mhvk commented Oct 3, 2016

StuartLittlefair commented Oct 17, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mhvk commented Oct 17, 2016

Choose a reason for hiding this comment

mhvk commented Oct 18, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mhvk commented Oct 19, 2016

mhvk commented Oct 20, 2016

astrobot commented Oct 20, 2016

eteq commented Oct 26, 2016

pllim commented Oct 26, 2016

MSeifert04 commented Oct 26, 2016 • edited

eteq commented Oct 26, 2016

eteq commented Oct 26, 2016

mhvk commented Oct 26, 2016

pllim commented Oct 27, 2016

MSeifert04 commented Oct 26, 2016 •

edited