Implement alt,az to ra,dec transform

[dkirkby]

Fixes #14. The following round trip from ra,dec to alt,az and back again is now included in the unit tests:

    # Specify the observing conditions
    where = observatories['APO']
    when = Time(56383, format='mjd')
    wlen = 5400 * u.Angstrom
    temperature = 5 * u.deg_C
    pressure = 800 * u.kPa
    obs_model = create_observing_model(where=where, when=when,
        wavelength=wlen, temperature=temperature, pressure=pressure)

    # Perform the round trip.
    sky_in = SkyCoord(ra=0.5*u.rad, dec=1.5*u.rad, frame='icrs')
    altaz_out = sky_to_altaz(sky_in, obs_model)
    sky_out = altaz_to_sky(altaz_out.alt, altaz_out.az, obs_model, frame='icrs')

    # Check that we get back the initial coordinates.
    assert np.allclose(sky_in.ra, sky_out.ra)
    assert np.allclose(sky_in.dec, sky_out.dec)

[dkirkby]

@sbailey Does the example above have the pieces you need for desisim (together with the alt,az <-> x,y transforms that are already implemented) ?

[sbailey]

To 0th order, it provides what I need, though the np.allclose calls fail for me with:

In [45]: assert np.allclose(sky_in.ra, sky_out.ra)
---------------------------------------------------------------------------
UnitsError                                Traceback (most recent call last)
<ipython-input-45-20dac39209c9> in <module>()
----> 1 assert np.allclose(sky_in.ra, sky_out.ra)

/Users/sbailey/anaconda/lib/python2.7/site-packages/numpy/core/numeric.pyc in allclose(a, b, rtol, atol, equal_nan)
   2272 
   2273     """
-> 2274     return all(isclose(a, b, rtol=rtol, atol=atol, equal_nan=equal_nan))
   2275 
   2276 def isclose(a, b, rtol=1.e-5, atol=1.e-8, equal_nan=False):

/Users/sbailey/anaconda/lib/python2.7/site-packages/numpy/core/numeric.pyc in isclose(a, b, rtol, atol, equal_nan)
   2352     yfin = isfinite(y)
   2353     if all(xfin) and all(yfin):
-> 2354         return within_tol(x, y, atol, rtol)
   2355     else:
   2356         finite = xfin & yfin

/Users/sbailey/anaconda/lib/python2.7/site-packages/numpy/core/numeric.pyc in within_tol(x, y, atol, rtol)
   2335     def within_tol(x, y, atol, rtol):
   2336         with errstate(invalid='ignore'):
-> 2337             result = less_equal(abs(x-y), atol + rtol * abs(y))
   2338         if isscalar(a) and isscalar(b):
   2339             result = bool(result)

/Users/sbailey/anaconda/lib/python2.7/site-packages/astropy/units/quantity.pyc in __array_prepare__(self, obj, context)
    315                                      "argument is not a quantity (unless the "
    316                                      "latter is all zero/infinity/nan)"
--> 317                                      .format(function.__name__))
    318             except TypeError:
    319                 # _can_have_arbitrary_unit failed: arg could not be compared

UnitsError: Can only apply 'add' function to dimensionless quantities when other argument is not a quantity (unless the latter is all zero/infinity/nan)

Casting to an array np.allclose(np.array(sky_in.ra), np.array(sky_out.ra)) works.

However, when testing with other values I find large errors in the round trip:

sky_in = SkyCoord(ra=(1+np.arange(10))*u.deg, dec=(1+np.arange(10))*u.deg, frame='icrs')
altaz_out = sky_to_altaz(sky_in, obs_model)
sky_out = altaz_to_sky(altaz_out.alt, altaz_out.az, obs_model, frame='icrs')
sky_in.ra - sky_out.ra
<Angle [ 0.0141732 , 0.00643555, 0.00321101, 0.00172966, 0.00099198,
         0.00059921, 0.00037803, 0.00024744, 0.00016715, 0.00011604] deg>

[dkirkby]

@sbailey Thanks for these checks!

The unit test I implemented does not have the np.allclose problem because it uses a single sky coordinate, not an array.

The largish (0.01 deg is not that large) errors you found are possibly because you are using an older version of astropy. With 1.0.4, I find:

>>> sky_in.ra - sky_out.ra
<Angle [ -3.38149508e-11, -3.39119843e-11, -3.40016904e-11,
         -3.41779938e-11, -3.44932971e-11, -3.49729135e-11,
         -3.56967789e-11, -3.66462416e-11, -3.79500875e-11,
         -3.95896649e-11] deg>

What version are you using? How does your altaz_out compare with mine:

>>> altaz_out
<SkyCoord (AltAz: obstime=2001-01-01T00:00:00.000, location=(-1994502.6043061386, -5037538.54232911, 3358104.9969029757) m, pressure=788.060134183 hPa, temperature=15.0 deg_C, relative_humidity=0.0, obswl=0.54 micron):00:00.000, location=(-1994502.6043061386, -5037538.54232911, 3358104.9969029757) m, pressure=788.060134183 hPa, temperature=15.0 deg_C, relative_humidity=0.0, obswl=0.54 micron): (az, alt) in deg
    [(157.7199411, 57.07493532), (155.3480744, 57.67822132),
     (152.90073054, 58.23672114), (150.38099892, 58.74858406),
     (147.79291974, 59.21207576), (145.14150998, 59.62561518),
     (142.43275748, 59.98781218), (139.673578, 60.29750454),
     (136.87173264, 60.55379242), (134.03570553, 60.75606855)]>

Also, did you use the same time and location as in my example?

[dkirkby]

@sbailey I just realized that I did the test above with zero pressure, with effectively turns off the refraction model. Using the correct observing conditions:

    where = observatories['APO']
    when = Time(56383, format='mjd')
    wlen = 5400 * u.Angstrom
    temperature = 5 * u.deg_C
    pressure = 800 * u.kPa
    obs_model = create_observing_model(where=where, when=when,
        wavelength=wlen, temperature=temperature, pressure=pressure)

I now get the same results as you for sky_in.ra - sky_out.ra. The reason for these large-ish errors is that your input sky coordinates are low on the horizon, where the refraction model is less physically reliable and less numerically stable:

>>> altaz_out.alt
<Latitude [  7.78756766,  9.04673426, 10.32969921, 11.63048343,
            12.94410003, 14.26670592, 15.59536474, 16.92779336,
            18.26216841, 19.59699526] deg>

I added a test_altaz_array_roundtrip unit test with a looser allclose absolute tolerance of 0.015 deg to make sure this does not get any worse, but I think we have to live with it.

[sbailey]

@dkirkby 0.015 deg = 54 arcsec which is quite a lot, though if this is only problematic at the horizon then it doesn't really matter for our purposes. I'd suggest a unit test with much tighter tolerance, while scanning around +/- 45 degrees of zenith, e.g.

from specsim.transform import sky_to_altaz, altaz_to_sky
from specsim.transform import observatories, create_observing_model
from astropy.coordinates import SkyCoord, ICRS
from astropy.time import Time
import astropy.units as u

where = observatories['APO']
when = Time(56383, format='mjd')
wlen = 5400 * u.Angstrom
temperature = 5 * u.deg_C
pressure = 800 * u.kPa
obs_model = create_observing_model(where=where, when=when,
    wavelength=wlen, temperature=temperature, pressure=pressure)

zenith = altaz_to_sky(90*u.deg, 0*u.deg, obs_model, frame='icrs')
ra = np.linspace(zenith.ra.value-45, zenith.ra.value+45, 11)*u.deg
dec = np.linspace(zenith.dec.value-45, zenith.dec.value+45, 11)*u.deg
altaz = sky_to_altaz(SkyCoord(ra, dec, frame='icrs'), obs_model)
radec = altaz_to_sky(altaz.alt, altaz.az, obs_model)
np.all( radec.separation(ICRS(ra, dec)) < 0.1*u.arcsec )

(which passes)

[dkirkby]

I took a more careful look at these round trip errors and added a notebook (under docs/nb/ but excluded in the packaging MANIFEST) to document this:

I added a unit test to verify that the round-trip altitude error is < 0.5 arcsec for altitude > 20 deg and updated the code to issue a UserWarning when an altitude < 20 deg is used (but this threshold can be changed via a model-level attribute).

@sbailey I think this is ready to merge now - do you agree?

[sbailey]

@dkirkby Plot looks good. Consider also adding a check that the error < 0.1 arcsec for altitude > 30 deg since that reflects the level of tolerance and range of telescope pointings that we will likely need. Otherwise looks good to merge. Thanks.

[dkirkby]

Ok, I added an additional check in the unit test for < 0.1 arcsec error above 30 deg altitude.

Keep in mind that the round trip error is measuring the numerical invertibility of the transformation, not the physical accuracy, which might be better or worse. The underlying refraction model comes from astropy, whose docs have the following note (ERFA has its own github repo):

The refraction model is based on that implemented in ERFA, which is fast but becomes inaccurate for altitudes below about 5 degrees. Near and below altitudes of 0, it can even give meaningless answers, and in this case transforming to AltAz and back to another frame can give highly discrepent results. For much better numerical stability, leaving the pressure at 0 (the default), disabling the refraction correction (yielding “topocentric” horizontal coordinates).

Also note that typical optical distortions across the focal plane are 1-10 arcsec, so we are only aiming for reasonably fast and straightforward sky transforms that are good enough for simulations here, not do-or-die operations.