astropy/astropy — pull request #3217 by eteq (2014-12-17 08:50:47+00:00)
Title: Implement ICRS<->AltAz transformation stack

As has long been promised, this PR implements the full stack necessary to do ICRS -> AltAz transformations. So it allows stuff like:
```
>>> m33 = SkyCoord.from_name('M31')
>>> newhaven = EarthLocation(lat=41.3*u.deg, lon=-73*u.deg,height=50*u.m)
>>> time = Time('2008-8-1 2:30:00')
>>> aaframe = AltAz(obstime=time,location=newhaven)
>>> 'alt={0.alt}, az={0.az}'.format(m33.transform_to(aaframe))
'alt=21.9927170464 deg, az=53.6342622223 deg'
```

As well as the following plot:
![out](https://cloud.githubusercontent.com/assets/346587/5468635/bcec1752-859f-11e4-8b17-a330bfb883b9.png)
which I made with:
```
>>> t2 = Time(time.jd + np.arange(1000)/1000.+.6,format='jd')
>>> dthr = (t2.jd - 2454680.5)*24
>>> aaframe2 = AltAz(obstime=t2,location=newhaven)
>>> sunazs = SkyCoord.get_sun(t2).transform_to(aaframe2)
>>> m33azs = m33.transform_to(aaframe2)
>>> plt.plot(dthr, sunazs.alt,c='y', label='Sun')
>>> plt.scatter(dthr, m31azs.alt,c = m31azs.az, lw=0, s=5, label='M33')
>>> plt.colorbar().set_label('Az [deg]')
>>> plt.legend(loc=0)
>>> plt.xlabel('hrs from 2 Aug 2008 UTC midnight')
>>> plt.ylabel('Alt [deg]')
>>> plt.xlim(dthr[0], dthr[-1])
>>> plt.ylim(0, 90)
```

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:

1. I relaxed the requirement that `TimeFrameAttribute`s 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 individual `SkyCoord`'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)
2. One thing that comes for "free" is the location of the sun (it's needed for various other steps). So I added `SkyCoord.get_sun`` to do just that. Combined with the point above, this makes the plot I've shown above possible. Note that that would have taken ~1 min to run with just a hundred points without the change from 1 above.
3. I added a `QuantityFrameAttribute` and` EarthLocationFrameAttribute`. The former might be useful for other frames (@adrn's #2761, perhaps?), and the latter is basically a way to allow earth locations to either use the existing `EarthLocation` *or* 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.
4. The ERFA transformations from ICRS to GCRS and from ITRS to AltAz throw away the distance. So to keep distances, this just uses regular cartesian distances. However, that makes it (strictly speaking) incorrect because it ignores relativistic effects. These are of course quite important for the direction of a celestial body (aberration and light deflection are crucial here), but less so for the distance. More importantly, we'd have to have a discussion about what we even *mean* by distance in a fully GR-consistent way. My feeling is that this is fine as-is, and we can always refine it in future versions if anyone ever cares to this level of precision.
5. There are not explicit accuracy tests. There are some spot-checks and a bunch of consistency tests (including comparing with raw ERFA), but no "I know what the alt/az and ra/dec are for these objects". Part of the problem is that there are a lot of subtlties in the various steps in the transformations, so you don't expect other people to agree exactly. But this is what astropy/coordinates-benchmark is for, so I'm deferring to that.

cc @astrofrog @adrn @mhvk @taldcroft