Add bary/heliocentric radial velocity correction by eteq · Pull Request #5752 · astropy/astropy

eteq · 2017-01-30T06:24:17Z

This PR implements functions for computing heliocentric velocities. The main practically useful bit is the helio_corr function. That computes the correction required to convert an observed radial velocity to a heliocentric velocity. This is a bread-and-butter calculation for basically all astronomical spectroscopic analysis. The neat thing is that the implementation here is pretty trivial now that the SS ephemerides stuff is in, and can even use JPL ephemerides if needed.

Now the question is: is this how we want to deal with velocities? Instead of doing these as piecemeal functions, this could be done as an actual transformation framework. That needs to be a separate discussion (I think there's already an issue about it?), but either way this code is fundamental. So I'm creating this PR now so that the code is here, and if we can get this figured out by 2.0 we put it in that framework, otherwise we merge this for 2.0 to make sure something is in there to do this fundamental activity.

One bonus bit to add if we go this route

add Solar System barycentric option - that's actually simpler than this, involving just skipping the bary->sun conversion step (essentially a 1-line change)

cc @adrn @mhvk @astrofrog

eteq · 2017-01-30T20:52:43Z

Hmm, a curious failure mode to matrix_product is revealed in the tests. @mhvk, do you have any idea what's going on in https://travis-ci.org/astropy/astropy/jobs/196493982 ?

mhvk · 2017-01-30T22:02:50Z

Haven't looked in much detail yet, but it looks like the inputs to matrix_product are not, in fact, matrices!

mhvk · 2017-01-30T23:33:58Z

Just to add back links: issue asking for barycentric correction: #3544; @StuartLittlefair's sample code: https://gist.github.com/StuartLittlefair/5aaf476c5d7b52d20aa9544cfaa936a1

(I'll try to look at the actual code here soon...)

mhvk · 2017-01-31T03:04:53Z

astropy/coordinates/velocities.py

+    gtarg = target.transform_to(GCRS(obstime=t, obsgeoloc=gcrs_p))
+    targxyz = gtarg.represent_as(UnitSphericalRepresentation).to_cartesian().xyz
+
+    res = matrix_product(vsuntarg_cartrepr.xyz, targxyz)


OK, what you really want here is the regular dot product, i.e., np.inner. Annoyingly, that function is not quantity-proof.

BUT we just changed representations to do this for us, so just do

gtarg /= gtarg.norm() res = gtarg.dot(vsuntarg_cartrepr)

@mhvk is almost right. Since you've already normalised gtarg through the conversion to a UnitSphericalRepresentation and back, you want:

targcart = gtarg.represent_as(UnitSphericalRepresentation).to_cartesian() res = targcart.dot(vsuntarg_cartrepr).to(KPS)

In my suggestion I wrote gtarg on purpose, i.e., I assumed the whole line targxyz would be removed (as it is unnecessarily slow compared to just dividing by the norm).

I went ahead and profiled it, and, quite surprisngly, the UnitSphericalRepresentation way is faster:

>>> %%timeit >>> targcart = gtarg.represent_as(UnitSphericalRepresentation).to_cartesian() >>>res = targcart.dot(vcorr_cart) 497 µs ± 14.1 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

and

>>> %%timeit >>> gtargd = gtarg.data >>> gtargd /= gtargd.norm() >>> res2 = gtargd.dot(vcorr_cart) 891 µs ± 28.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

That surprises me, but there it is... so will go that route.

OK, real data beats everything! I also think the route is also clearer than the original - it still avoids the expensive .xyz. Note that you can remove the .to_cartesian() -- .dot will do that for you.

But of course need to understand: represent_as(USph) is fast - just a new class initialization (this can become faster still...). Then, to_cartesian is slightly faster than it would have been before, since no multiplication with a distance (of 1.). While for my original suggestion there is the division by the distance (quite fast, but needless), a similar reinitialization, and then a slightly more expensive conversion to cartesian. My timings are not as different, but still the route via unitspherical is clearly better:

g = SphericalRepresentation(np.arange(0,100)*u.deg, np.arange(-45., 55.)*u.deg, np.arange(100,200)*u.kpc) v = CartesianRepresentation(5., 5., 5., unit=u.km/u.s) %timeit (g/g.norm()).dot(v) # 1000 loops, best of 3: 1.69 ms per loop %timeit g.represent_as(UnitSphericalRepresentation).dot(v) # 1000 loops, best of 3: 1.49 ms per loop %timeit g.represent_as(UnitSphericalRepresentation).to_cartesian().dot(v) # 1000 loops, best of 3: 1.49 ms per loop

StuartLittlefair · 2017-01-31T11:20:33Z

Hi @eteq,

I think it could be a good idea to implement this code now and worry about integration with a framework later. I'll add some inline style comments to the code, but I have a larger concern about precision.

If you read this paper - and in particular pp 842, you'll see that the approach we're taking here of adding and subtracting velocities, cannot be more accurate than a few m/s. Should we be aiming for better than this? Without cm/s precision, the exoplanet folk won't be interested in this, whereas the time corrections are good enough for all but pulsar astronomers. Given the exoplanet community is possibly the largest user base for this code, I'd argue we should be implementing a more precise algorithm.

On the subject of testing precision - the algorithm in the paper above has an online calculator which could be used to generate barycentric corrections for testing - http://astroutils.astronomy.ohio-state.edu/exofast/barycorr.html

StuartLittlefair · 2017-01-31T11:23:33Z

astropy/coordinates/velocities.py

+    gcrs_p, gcrs_v = loc.get_gcrs_posvel(t)
+
+    vsuntargxyz = vsunearth.xyz + gcrs_v.xyz
+    return CartesianRepresentation(vsuntargxyz.to(KPS))


This isn't needed any more, you should be able to do:

return vsunearth + gcrs_v

unless you really want it in km/s at this point

StuartLittlefair · 2017-01-31T11:31:04Z

astropy/coordinates/velocities.py

+    vsuntarg_cartrepr = helio_vector(t, loc)
+    gcrs_p, _ = loc.get_gcrs_posvel(t)
+
+    gtarg = target.transform_to(GCRS(obstime=t, obsgeoloc=gcrs_p))


For some reason I find this confusing, but shouldn't you be using the ICRS co-ordinate of the target here, and finding the dot product with the velocity vector? It will only make a ~1 m/s difference but it will be important if we shoot for a more precise algorithm

Hmm... I don't think that's right. It's sort of in the definition of the correction, but the way I think about barycentric correction is: "the velocity vector of the earth relative to the barycenter, projected along the line-of-sight from the observer to the target". If I used the ICRS coordinate instead of the GCRS coordinate, that would be the line-of-sight from the barycenter to the target. For anything beyond the solar system this is a tiny effect... but it is still what I would say is the correction.

Should probably be explicit about this in the docs, though.

I think @eteq is right: consider a not too-distance star that is exactly at the ecliptic pole, and which has a proper motion. As the earth orbits the sun, this proper motion would induce a radial-velocity modulation, because as seen from the geocenter, the position of the star oscillates around the ecliptic pole. But in ICRS, the star doesn't move.

StuartLittlefair · 2017-01-31T11:32:50Z

astropy/coordinates/velocities.py

+    gtarg = target.transform_to(GCRS(obstime=t, obsgeoloc=gcrs_p))
+    targxyz = gtarg.represent_as(UnitSphericalRepresentation).to_cartesian().xyz
+
+    res = matrix_product(vsuntarg_cartrepr.xyz, targxyz)


@mhvk is almost right. Since you've already normalised gtarg through the conversion to a UnitSphericalRepresentation and back, you want:

targcart = gtarg.represent_as(UnitSphericalRepresentation).to_cartesian() res = targcart.dot(vsuntarg_cartrepr).to(KPS)

StuartLittlefair · 2017-01-31T11:34:14Z

astropy/coordinates/velocities.py

+    return CartesianRepresentation(vsuntargxyz.to(KPS))
+
+
+def helio_corr(t, loc, target):


In my gist I used the signature (Time, SkyCoord, EarthLocation=None), to be consistent with how we implemented time corrections. Again, for consistency, if the EarthLocation is None, we try and get the location from the Time object.

I'm actually really against any location information from the Time object. The whole point of having EarthLocation existing is to get around the awkwardness there, and its more natural anyway for something in coordinates. But the rest sounds fine!

We had a long discussion for the barycentric time corrections about where these things should live. For velocity corrections, to me SkyCoord is more logical, in which the signature would become (SkyCoord, Time, Location) - with defaults of sc.obstime and sc.location (and, if a time is passed in, a secondary default of time.location).

That said, I think it would be fine to have a short-cut on Time (and similarly a short-cut for time corrections on SkyCoord).

If we do add it to Time as well, I would stick to @StuartLittlefair's signature for consistency with light_travel_time - it seems logical to look for defaults first on self.

mhvk · 2017-01-31T14:15:43Z

On the implementation: I would suggest initially to follow the precedent of the barycentric time corrections, and make the public API be a method rather than a stand-alone version. The method is probably most logical on SkyCoord (and one may want the ability to calculate the time difference there too, even though it breaks the "only one way" rule). If so, the most logical signature might be

def radial_velocity_correction(self, obstime=None, kind='barycentric', location=None, ephemeris=None)

where obstime would by default get taken from the coordinate, as would location (from obsgeoloc, with obstime.location as a secondary default). Not completely sure whether it should be location or obsgeoloc -- I do not like the latter for its abbreviation, but it is what is stored in some SkyCoord.

mhvk · 2017-01-31T14:17:03Z

Longer term, I guess we'd want SkyCoord to be able to have both proper motion and radial velocity. It may get slightly tricky to deal with situations where some of those components (and/or distance) are missing, but we have to face this one way or the other...

mhvk · 2017-01-31T14:22:59Z

@StuartLittlefair - Looking at the reference you cited, it would seem what we need to provide is just the velocity correction, but we have to make clear in the documentation that one applies this as

v_corrected = c*[(1+v_measured/c) * (1 + v_correction/c) - 1.]

(and we'll need to implement that ourselves later...)

StuartLittlefair · 2017-01-31T14:56:04Z

Yes, that's true.

However, I think we are still missing important (i.e m/s) physics. There's a python interface to the online calculator I linked to earlier. If you install that you can see that neither my gist, nor @eteq 's code agrees with the results of BARYCORR to a few m/s.

This gist shows an example...

mhvk · 2017-06-12T16:40:13Z

@eteq - changed the milestone to 3.0 for now -- I think we have enough on our plate without this (and if you get the velocities right, this comes for free with the proper coordinate transform...).

eteq · 2017-06-16T08:00:25Z

Note that, while in principal this functionality is now supported (with a lot more) in #6226 and the PR's that it are derived from, the finite difference calculation in that PR has some numerical stability issues for some relevant use cases. So I think the best thing to do is include this PR (updated with the comments, of course), with a note in the docs that it may at some point be adjusted to use the full velocity framework if the numerical stability problems can be solved.

Will update this PR ASAP with that plan in mind.

StuartLittlefair · 2017-06-16T09:56:10Z

Great @eteq - perhaps add a note in the docs that this PR doesn't yet provide m/s precision and link to the barycorr Python module? Otherwise there's a risk we'll get lynched by some exoplanet people...

mhvk · 2017-06-16T11:45:23Z

I'm still rather uncomfortable with the API here, with a new module. Since this is rather simple, I'd much rather have a method on SkyCoord (or possibly Time).

More generally, I think with the velocities we have more than enough on our plate, and I prefer to think about this better -- another option is to simply document this as part of solar_system.

eteq · 2017-06-16T22:18:07Z

I'm about to try to re-work this with the implementation comments above in mind, but just a couple thoughts on the bigger picture from @StuartLittlefair:

I think reaching the cm/s precision will take a lot more care than we have done thus far, as you see. Physics aside, our default ephemerides via erfa aren't good enough for that right now. I think this would be a great project to sic someone on, though... once we have the start of the framework in hand.
Agreed that this means it should come with a clear warning of its precision level. But actually there are plenty lot of use cases for this that are not exoplanet folks - pretty much anyone who wants to do Galactic or extragalactic work with spectra has to do this, and km/s level precision is usually more than sufficient for those cases. So hopefully we can get them something now, and build up to the higher precision case.

eteq · 2017-06-16T22:22:28Z

@mhvk - I'm receptive to the idea of it being in a method, as you suggest... the main concern, though, is that it might lead to a lot more confusion with the velocities work @adrn and I are trying to finish because there'd be two ways to do basically the same thing. A function, meanwhile, is much easier to cast as sort of a "temporary" workaround or an alternate approach (even if we support it indefinitely by having it eventually be based on the properly-done coordinates stuff). People will be much more likely to see it as a distinct implementation.

…iews

eteq · 2017-06-20T05:07:49Z

With the commits I just pushed up I think I addressed all the comments. I did not add an ephemeris keyword, but did add an example of how to use the existing machinery to use an alternate ephemeris.

I think this is very close now (assuming tests pass), but two remaining unresolved items:

In the current version I am disallowing any time.location, on the assumption that that's just too many layers of indirection/ambiguity. I could do the chain @mhvk suggested (location arg, if not: SkyCoord.location, if not: obstime.location, if not: SkyCoord.obstime.location). That seems likely to cause a lot of unforeseen oddities, though...
In the latest commit, obstime/location cannot be passed in as arguments if they're on the SkyCoord. But this means that if a SkyCoord is in a frame that requires obstime or location (most important of these is AltAz which needs both), there' s no way to change obstime other than changing it on the SkyCoord). Are we all ok with that? (I am, because I think it demonstrates the original reason for not having it on the method... but I just wanted to make sure we all understood the challenge here)

I think @adrn and @mhvk might want to respond to those, so won't merge now. But this is very close so a clear candidate for backporting later this week, and hence I won't push the milestone.

StuartLittlefair · 2017-06-20T08:48:39Z

On disallowing time.location: my problem with that is that it is totally inconsistent with how we have implemented light_travel_time. Indeed the docs instruct the user to put the location onto the Time object before calculating light_travel_time corrections. I can understand you wanting to enforce some rigour here, but I think the user is entitled to feel frustrated if their next step is to calculate RV corrections and they have make a new Time object with the location removed, and add it to their SkyCoord.

Other than elegance, is there any real issue with having a clear hierarchy, e.g SkyCoord.location > Time.location > location, as long as it is clearly documented?

On the second restriction, I can only think of very contrived use cases where you would want to use a different obstime or location - at least for the AltAz case. Again, my issue here is one of convenience to the user.

mhvk · 2017-06-20T13:11:17Z

@StuartLittlefair's point I think is valid: typical use in my "ICRS coordinate case" would be to calculate both barycentric corrections for the same Time -- I think we might as well support the location inside Time, although I think we're OK with insisting that SkyCoord does not have an obstime, and ignoring any location inside SkyCoord.obstime (since it should, if used, be already copied to SkyCoord.location).

…correction

eteq · 2017-06-21T06:10:48Z

OK, I updated this so that it will use obstime.location if it's given and neither of the other are present. Is this now what you were thinking, @StuartLittlefair? I was not entirely clear if you wanted to allow there to be multiple obstime's and just pick the one higher in "precedence", or if your > symbols meant a recommendation?

I'm agains the idea of allowing any combination (e.g. obstime.location and SkyCoord.location are both not None), because there's just no way to tell what the user intended - did they mean to override, or just not realize that there's a location on that time? (similar for SkyCoord). It's easy enough to get such objects from other contexts without specifically attaching the location, so I think it's better to be conservative, After all, "In the face of ambiguity, refuse the temptation to guess." ™️... But I am receptive to the convenience argument, so if there is a common use case where you definitely expect both to be present I can live with that.

mhvk

This now looks all OK to me. Missing though, is documentation.

eteq · 2017-06-21T18:23:30Z

@mhvk - good point re: docs, but my plan was to include this in the section that is to be added #6219 - so that would be easiest post-merge, and possibly in between the beta and final release. We should create an issue accounting for that once this is merged and make sure its milestoned 2.0.

(Although there's a real test failure in the docstring which I'll fix shortly)

adrn · 2017-06-21T19:54:40Z

I'm fine with adding this to the (coming) velocity documentation, but yes, an issue milestoned 2.0 is required!

mhvk · 2017-06-21T20:54:29Z

Same, fine to add later, as long as there is a 2.0 issue for it.

[docs only]

eteq · 2017-06-22T05:37:37Z

I added a note about adding the documentation for this in #6260. Assuming the tests pass, I think this is now ready!

StuartLittlefair · 2017-06-22T14:54:47Z

Apologies. Thought I'd replied - looks fine, merge away!

eteq · 2017-06-22T14:55:36Z

I'm going to go ahead and merge this, but @StuartLittlefair, if I was misunderstanding your intent above, feel free to raise it as a separate issue and/or a PR of your own. If we can get to that by next Monday we can still pull it into 2.0 as that would be a "bug" in my interpretation of your comment 😉

astrobot · 2017-06-22T14:56:02Z

@eteq - 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!

If you believe the above to be incorrect (which I - @astrobot - very much doubt) you can ping @astrofrog

eteq · 2017-06-22T14:57:20Z

Whoops! Forgot the changelog... I think I'll just backport this to v2.0.0 now then and fix it directly in master to avoid confusion...

Add bary/heliocentric radial velocity correction

eteq · 2017-06-23T04:57:03Z

(commits listed above are for the changelog in master, the backport of this PR to v2.0.x, and the backport of the changelog commit, respectively)

eteq added coordinates Experimental labels Jan 30, 2017

eteq added this to the v2.0.0 milestone Jan 30, 2017

eteq requested a review from adrn January 30, 2017 06:24

eteq changed the title ~~Add heliocentric correction~~ Add heliocentric radial velocity correction Jan 30, 2017

mhvk reviewed Jan 31, 2017

View reviewed changes

StuartLittlefair reviewed Jan 31, 2017

View reviewed changes

mhvk modified the milestones: v3.0.0, v2.0.0 Jun 12, 2017

eteq mentioned this pull request Jun 16, 2017

Support Velocities in coordinates by having transforms that use finite differences #6226

Merged

eteq added 6 commits June 16, 2017 18:27

add helio velocities correction

7ffd296

add tests

7e5b1ef

fixes for code based on tests

34d70fe

work around apparent matrix_product bug

f37e52e

update to radial_velocity_correction and give bary option

7c5f914

update to new implementation based on @mhvk and @StuartLittlefair rev…

c45114a

…iews

eteq added 3 commits June 20, 2017 00:22

clean up if/else structures and be more consistent

7814282

updates to radial_velocity_correction from review

18a1c14

fix some typos

1541318

Update pathway for choosing which location to use in radial_velocity_…

48f2689

…correction

mhvk approved these changes Jun 21, 2017

View reviewed changes

adrn approved these changes Jun 22, 2017

View reviewed changes

Fix sphinx warnings

2bdc92b

[docs only]

eteq mentioned this pull request Jun 22, 2017

Add documentation for velocity support in coordinates #6260

Merged

eteq added the zzz 💤 merge-when-ci-passes Do not use: We have auto-merge option now. label Jun 22, 2017

eteq merged commit 3d53b73 into astropy:master Jun 22, 2017

eteq deleted the add-vhelio branch June 22, 2017 14:56

eteq modified the milestones: v2.0.0, v3.0.0 Jun 22, 2017

eteq added a commit that referenced this pull request Jun 22, 2017

Added missing changelog entry for #5752

cb93baa

eteq added a commit that referenced this pull request Jun 22, 2017

Merge pull request #5752 from eteq/add-vhelio

8e0ecb8

Add bary/heliocentric radial velocity correction

eteq added a commit that referenced this pull request Jun 22, 2017

Added missing changelog entry for #5752

132db4d

StuartLittlefair mentioned this pull request Nov 16, 2017

Precise barycentric radial velocity corrections. #6861

Merged

		return CartesianRepresentation(vsuntargxyz.to(KPS))


		def helio_corr(t, loc, target):

Uh oh!

Conversation

eteq commented Jan 30, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

eteq commented Jan 30, 2017

Uh oh!

mhvk commented Jan 30, 2017

Uh oh!

mhvk commented Jan 30, 2017

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

StuartLittlefair commented Jan 31, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

mhvk commented Jan 31, 2017

Uh oh!

mhvk commented Jan 31, 2017

Uh oh!

mhvk commented Jan 31, 2017

Uh oh!

StuartLittlefair commented Jan 31, 2017

Uh oh!

mhvk commented Jun 12, 2017

Uh oh!

eteq commented Jun 16, 2017

Uh oh!

StuartLittlefair commented Jun 16, 2017

Uh oh!

mhvk commented Jun 16, 2017

Uh oh!

eteq commented Jun 16, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

eteq commented Jun 16, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

eteq commented Jun 20, 2017

Uh oh!

StuartLittlefair commented Jun 20, 2017

Uh oh!

mhvk commented Jun 20, 2017

Uh oh!

eteq commented Jun 21, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

mhvk left a comment

Choose a reason for hiding this comment

Uh oh!

eteq commented Jun 21, 2017 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

adrn commented Jun 21, 2017

eteq commented Jan 30, 2017 •

edited

Loading

StuartLittlefair commented Jan 31, 2017 •

edited

Loading

eteq commented Jun 16, 2017 •

edited

Loading

eteq commented Jun 16, 2017 •

edited

Loading

eteq commented Jun 21, 2017 •

edited

Loading

eteq commented Jun 21, 2017 •

edited

Loading