Implement ICRS<->AltAz transformation stack #3217
Conversation
|
ping @dpshelio |
|
@cdeil - if you have a chance to work on the coordinates-benchmark, can you also test out this branch against the other codes? |
There was a problem hiding this comment.
@eteq - Can you extract this generation of test coordinates into a utility function?
Below you are using the same code.
In Gammapy I have one which randomly samples the sphere ... feel free to copy it if you think it's useful:
https://gammapy.readthedocs.org/en/latest/_modules/gammapy/utils/random.html#sample_sphere
Also I don't see the point why grids of tests points on the sphere should be random, even if it's probably generating the same coordinates on almost all machines.
Instead a function that generates a grid of values that includes critical points like pole and lon=0 line where numerical problems can occur would be better.
But the first quick step is to extract this into a function (generate_coordinates_test_grid or something that has a few options, including N_TO_TEST) to have it in a central place.
There was a problem hiding this comment.
@cdeil - I've implemented abstracting out the coordinate-generation and making it uniform on the sphere instead of in ra/dec - will push up the commits when I finish other work on this PR.
But I don't think it makes sense to use a grid: I did that for one of the earlier coordinate tests and ran into situations where chance coincidences of numbers (basically, a form of aliasing) would occasionally cause tests to succeed where they shouldn't have. Random points, though, don't have those sorts of failure modes.
On repeatability, my interpretation is that NumpyRNGContext(12345) should always give the same answer, regardless of architecture. So it's still 100% repeatable on every platform. (If that's not true, presumably its a bug!)
|
@eteq This is huge, thanks! @astrofrog - I'll try to find time this afternoon and tomorrow to develop precision tests for this in My guess would be that it will take weeks and discussions with coordinate experts to set up precision tests for the large grid of transforms and cases and to understand all the comparisons against other packages ... if there's bugs in this implementation they can be fixed in later PRs. |
CHANGES.rst
Outdated
There was a problem hiding this comment.
SkyCoord.get_sun seems a bit weird to me ... will you add SkyCoord.get_mars and SkyCoord.get_moon as well?
If this has been discussed, OK, no strong objections, but if it's still up for discussion I'd prefer to keep SkyCoord simple and put the function to get the solar position somewhere else.
There was a problem hiding this comment.
I have a feeling that this will be very useful for us, however, it does seem a little odd to have it attached to SkyCoord why not a separate function?
ping @dpshelio
There was a problem hiding this comment.
Agreed that this is great to have but that it might be best to have it as a separate function, at least for now. Eventually we'll have access to the JPL ephemerides and therewith also to planets, as well as even more accurate solar positions.
There was a problem hiding this comment.
@mhvk What's the plan for JPL ephemerides? These would be unspeakably useful!
There was a problem hiding this comment.
Yes, and pulsar timing needs it too. I've been hoping for quite a while to port https://github.com/brandon-rhodes/python-jplephem over to astropy (as discussed with @brandon-rhodes quite a while back). My own sense right now is to first get things shaken out with erfa, including barycentric and heliocentric corrections (see #2244), and then add a high-precision option using the JPL ephemerides.
... actually let me open an issue so we don't forget... #3219..
There was a problem hiding this comment.
My reasoning for putting it in SkyCoord was that the sun is a somewhat special case because it's location is tied in with the definitions of coordinate frames in a way that other ephemerides are not. Also, for optical observing it has special significance, of course. So my thinking was that in v1.1 we'd implement a more general ephemerides approach (shall we move discussion of that to #3219 ?).
But it seems pretty clear that everyone here would rather have it an external function, so I'll just move it to funcs.py for now.
There was a problem hiding this comment.
Sounds good. It may be that we end up reconsidering (your reasoning is certainly logical), but earlier to add a method later than to remove it!
|
@eteq - The example calculation and plot you give in this PR is great! Can you copy it into the Sphinx docs? |
|
Fantastic! 👍 on the vectorisation for |
|
@eteq - very nice work! I've had a chance to try this out while traveling and I've outlined a few issues below. Bug in GCRS to AltAzI have tested this with a few cases and ran into a problem due to degeneracies in the transform graph. Basically if one does: c3 = SkyCoord.get_sun(times).transform_to(AltAz(...))The original and:: These are the functions in the composite transform obtained with: frame_transform_graph.get_transform(GCRS, AltAz)Now if I run the following example: import numpy as np
from astropy.coordinates import SkyCoord, EarthLocation, AltAz
from astropy import units as u
from astropy.time import Time
times = Time(np.linspace(2456293.25, 2456657.25, 8575) * u.day,
format='jd', scale='utc')
loc = EarthLocation(lon=10 * u.deg, lat=80. * u.deg)
aaframe = AltAz(obstime=times, location=loc)
c3 = SkyCoord.get_sun(times).transform_to(aaframe)I sometimes get a failure, sometimes I don't, and this is due to the degeneracy in the transform graph. The exception only happens when going through ITRS, and is: This is just an issue with the order of the dimensions in the dot product. It took me a while to figure out what was causing this because the error was stochastic, but once I realized that the transform graph is where the stochasticity comes from it helped. So if you can force the transformation to go this way, it should be easy to debug. Degeneracies in transform graphA secondary recommendation is that we adjust the weights in the transform graph so that the fastest route is always picked. In fact, maybe we could actually do some speed benchmarks of each branch in the transform graph and adjust the weights accordingly? This will incidentally also ensure that there are no degeneracies. Of course the speed will depend on the machine, but if we use a 'typical' machine it should be fine for most cases. This is probably too much effort for now, but I would at least advocate that we adjust weights as needed to avoid degenerate paths. Generalize
|
There was a problem hiding this comment.
altitudes below 5 degrees -> do you mean latitudes, or is the unit wrong?
There was a problem hiding this comment.
This surely is correct: how high the source is over the horizon (too low, and refraction is more complicated -- no equivalent of the green flash!).
There was a problem hiding this comment.
altitude is right. But to avoid confusion, I could edit this to say something like "altitudes below 5/zenith angles above 85". This is why ERFA uses zenith angles, apparently... too bad astronomers are slaves to tradition 😉
There was a problem hiding this comment.
Ahh, sorry, of course, I was thinking altitude as in height ;)
|
@astrofrog - just quickly on |
astropy/utils/iers/iers.py
Outdated
There was a problem hiding this comment.
This will never get called, but still probably better to return np.zeros_like(i).
astropy/utils/iers/iers.py
Outdated
There was a problem hiding this comment.
Ah, this is neat, automatically taking the unit from the Column. Should do the same for ut1_utc, though obviously not urgent (#3223).. For a bit of speed-up, you may want to add copy=False in the Quantity initialisation.
|
@eteq - just looked over the |
|
On |
|
To start with: I confess that I am shocked at this, but the tests passed on the first try! @mhvk thanks for the review! re:
|
|
I've worked my way through @astrofrog's comments (except that I need to add the example), and pushed up the work so far in case @cdeil or @Cadair want to examine some of it. More coming soon, though. |
|
I don't have much to add about the code, but I just want to note that this is extremely cool and I think we should definite write a tutorial on "preparing for an observing run" once this is merged. I'm happy to spearhead the tutorial, but drop any ideas about what would go in such a tutorial here: |
eteq
force-pushed
the
coords-icrs2altaz
branch
from
702e9fd to
7d8128b
Compare
|
I think all the comments have been addressed, so this is good to go if the tests pass. @cdeil - I've xfailed some of the tests you have in there because they're incomplete as you said. We can give that a more thorough shakeout in beta (as at that point it would be bug-fixing) |
|
I had a final look at this and it is good to go. Once this is merged we need to spend some time verifying the accuracy (which is where the coordinates-benchmark works comes in). |
Implement ICRS<->AltAz transformation stack
|
@cdeil - just for info, this is merged, so it should make the coordinates-benchmark work easier to do. |
|
Opened an issue to remind us that accuracy tests are critical before 1.0 is released: #3244 |
|
It looks like there's no example using @eteq Do you have time to add one? |
|
@cdeil - an expanded version of the examples is in http://astropy.readthedocs.org/en/latest/coordinates/observing-example.html But its weird you don't see that with the search... I just pushed a commit up that makes this document slightly easier to find (it's now mentioned in the "getting started" section). Maybe that will make the search work...? we'll see when the docs finish building. |
|
@eteq - Thanks, this is a very nice example! The reason it doesn't show up in the Sphinx search is because apparently only normal text is included, not code or links. E.g. for So It's still pretty hard to find and unless there's a way to configure the Sphinx search to include references, only mentioning keywords like "horizontal", "altitude", "attitude" in the main text would help. One minor issue I noticed is that if you run the example you get this warning that goes on for pages ... can you maybe not print the Time object fully here? |
7 participants
As has long been promised, this PR implements the full stack necessary to do ICRS <-> AltAz transformations. So it allows stuff like:
As well as the following plot:
which I made with:
While it's in principal a lot of code, a fair amount is testing, and most of the "guts" of all of it is from ERFA, so it's already well tested. I few items of note that reviewers might want to focus on:
TimeFrameAttributes have to be scalar. The logic here is that for celestial->altaz, you sometimes want the same celestial thing at different times. Doing it this way lets us use the vectorized ERFA machinery, which is much faster than a bunch of individualSkyCoord's. However, it may not work for all of the non-ERFA transforms. Right now I have it giving a warning that says basically "you might encounter problems with some transforms if you do this." We can fix that up in a future release, but it's a huge step forward to allow this for e.g. observation planning (see 2 below)QuantityFrameAttributeandEarthLocationFrameAttribute. The former might be useful for other frames (@adrn's Implement galactocentric frame #2761, perhaps?), and the latter is basically a way to allow earth locations to either use the existingEarthLocationor a coordinate frame (as long as it transforms to ITRS). This is a neat trick because it's sort of "nested coordinates" but shows the power of the flexible transformation scheme.cc @astrofrog @adrn @mhvk @taldcroft @Cadair