Merge ebed159 into 5561f94

.gitignore

@@ -20,3 +20,4 @@ hisparc_sapphire.egg-info/
 ehthumbs.db
 Icon?
 Thumbs.db
+

.travis.yml

@@ -20,7 +20,7 @@ before_install:
   - conda update --yes conda
 
 install:
-  - conda create --yes -n build_env python=$TRAVIS_PYTHON_VERSION numpy scipy pytables numba
+  - conda create --yes -n build_env python=$TRAVIS_PYTHON_VERSION numpy scipy pytables numba astropy
   - source activate build_env
   - pip install coverage coveralls codecov
   - pip install -e .[dev]

sapphire/tests/transformations/test_celestial.py

@@ -1,12 +1,21 @@
 import calendar
 import datetime
 import unittest
+
 from math import pi
 
 import numpy as np
 
 from sapphire.transformations import base, celestial, clock
 
+# This is to switch off tests in case astropy is not present
+# Noqa used to silence flake8
+try:
+    import astropy  # noqa : F401
+    has_astropy = True
+except ImportError:
+    has_astropy = False
+
 
 class ZenithAzimuthHorizontalTests(unittest.TestCase):
 
@@ -19,12 +28,16 @@ def setUp(self):
 
     def test_zenithazimuth_to_horizontal(self):
         for zenith, altitude in zip(self.zenith, self.altitude):
-            self.assertEqual(celestial.zenithazimuth_to_horizontal(zenith, 0)[0], altitude)
-            self.assertEqual(celestial.horizontal_to_zenithazimuth(altitude, 0)[0], zenith)
+            self.assertEqual(celestial.zenithazimuth_to_horizontal(
+                zenith, 0)[0], altitude)
+            self.assertEqual(celestial.horizontal_to_zenithazimuth(
+                altitude, 0)[0], zenith)
 
         for azimuth, Azimuth in zip(self.azimuth, self.Azimuth):
-            self.assertEqual(celestial.zenithazimuth_to_horizontal(0, azimuth)[1], Azimuth)
-            self.assertEqual(celestial.horizontal_to_zenithazimuth(0, Azimuth)[1], azimuth)
+            self.assertEqual(celestial.zenithazimuth_to_horizontal(
+                0, azimuth)[1], Azimuth)
+            self.assertEqual(celestial.horizontal_to_zenithazimuth(
+                0, Azimuth)[1], azimuth)
 
 
 class EquatorialTests(unittest.TestCase):
@@ -34,7 +47,8 @@ class EquatorialTests(unittest.TestCase):
     Use references as tests also used by astropy.
 
     Source:
-    https://github.com/astropy/astropy/blob/master/astropy/coordinates/tests/accuracy/test_altaz_icrs.py
+    https://github.com/astropy/astropy/blob/master/
+    astropy/coordinates/tests/accuracy/test_altaz_icrs.py
 
     Licensed under a 3-clause BSD style license - see:
     https://github.com/astropy/astropy/blob/master/licenses/LICENSE.rst
@@ -82,11 +96,18 @@ def test_against_hor2eq(self):
         ra, dec = celestial.zenithazimuth_to_equatorial(latitude, longitude,
                                                         gps, zenith, azimuth)
 
+        # Test eq_to_zenaz merely against IDL
+        zencalc, azcalc = celestial.equatorial_to_zenithazimuth(
+            latitude, longitude, gps, ra_expected, dec_expected)
+
         self.assertAlmostEqual(ra, ra_expected, 1)
         self.assertAlmostEqual(ra, ra_astropy, 1)
         self.assertAlmostEqual(dec, dec_expected, 2)
         self.assertAlmostEqual(dec, dec_astropy, 2)
 
+        self.assertAlmostEqual(zencalc, zenith, 1)
+        self.assertAlmostEqual(azcalc, azimuth, 2)
+
     def test_against_pyephem(self):
         """Check for consistent results with one PyEphem example.
 
@@ -115,14 +136,20 @@ def test_against_pyephem(self):
         # SAPPHiRE
         gps = clock.utc_to_gps(calendar.timegm(utc.utctimetuple()))
         zenith, azimuth = celestial.horizontal_to_zenithazimuth(altitude, azi)
-        ra, dec = celestial.zenithazimuth_to_equatorial(latitude, longitude,
-                                                        gps, zenith, azimuth)
+        ra, dec = celestial.zenithazimuth_to_equatorial(
+            latitude, longitude, gps, zenith, azimuth)
+
+        zencalc, azcalc = celestial.equatorial_to_zenithazimuth(
+            latitude, longitude, gps, ra_expected, dec_expected)
 
         self.assertAlmostEqual(ra, ra_expected, 2)
         self.assertAlmostEqual(ra, ra_astropy, 2)
         self.assertAlmostEqual(dec, dec_expected, 2)
         self.assertAlmostEqual(dec, dec_astropy, 2)
 
+        self.assertAlmostEqual(zencalc, zenith, 2)
+        self.assertAlmostEqual(azcalc, azimuth, 2)
+
     def test_against_jpl_horizons(self):
         """Check for consistent results with the JPL Horizons example.
 
@@ -153,11 +180,145 @@ def test_against_jpl_horizons(self):
         ra, dec = celestial.zenithazimuth_to_equatorial(latitude, longitude,
                                                         gps, zenith, azimuth)
 
+        zencalc, azcalc = celestial.equatorial_to_zenithazimuth(latitude,
+                                                                longitude, gps,
+                                                                ra_expected,
+                                                                dec_expected)
+
         self.assertAlmostEqual(ra, ra_expected, 3)
         self.assertAlmostEqual(ra, ra_astropy, 3)
         self.assertAlmostEqual(dec, dec_expected, 2)
         self.assertAlmostEqual(dec, dec_astropy, 2)
 
+        self.assertAlmostEqual(zencalc, zenith, 2)
+        self.assertAlmostEqual(azcalc, azimuth, 2)
+
+
+class AstropyEquatorialTests(unittest.TestCase):
+    """
+    This tests the 4 new astropy functions. They should be very close to
+    Pyephem results and in this test they are compared to 10 different
+    coordinates from astropy.
+
+    """
+    @unittest.skipUnless(has_astropy, "astropy required.")
+    def test_pyephem_htoea(self):
+        """ Check celestial.horizontal_to_equatorial_astropy """
+
+        # This is the transform inputs
+        eq = [(-39.34633914878846, -112.2277168069694, 1295503840,
+               3.8662384455822716, -0.31222454326513827),
+              (53.13143508448587, -49.24074935964933, 985619982,
+               3.901575896592809, -0.3926720112815971),
+              (48.02031016860923, -157.4023812557098, 1126251396,
+               3.366278312183976, -1.3610394240813288)]
+
+        # result of pyephem hor->eq/zenaz-> eq
+        efemeq = [(5.620508199785029, -0.3651173667585858),
+                  (5.244630787139936, -0.7866376569183651),
+                  (2.276751381056623, -1.0406498066785745)]
+
+        htoea_test = []
+
+        # Produce horizontal_to_equatorial_astropy results
+        for latitude, longitude, gps, az, alt in eq:
+            result = celestial.horizontal_to_equatorial_astropy(
+                latitude, longitude, gps, [(az, alt)])
+            htoea_test.extend(result)
+
+        # Check if all inputs are correct
+        # Test horizontal_to_equatorial_astropy
+        np.testing.assert_almost_equal(efemeq, htoea_test, 4)
+
+    @unittest.skipUnless(has_astropy, "astropy required.")
+    def test_pyephem_etoha(self):
+        """Check celestial.equatorial_to_horizontal_astropy"""
+
+        # This is the transform inputs
+        eq = [(-39.34633914878846, -112.2277168069694, 1295503840,
+               3.8662384455822716, -0.31222454326513827),
+              (53.13143508448587, -49.24074935964933, 985619982,
+               3.901575896592809, -0.3926720112815971),
+              (48.02031016860923, -157.4023812557098, 1126251396,
+               3.366278312183976, -1.3610394240813288)]
+        # result of pyephem eq->hor
+        altaz = [(2.175107479095459, -0.19537943601608276),
+                 (5.25273323059082, -0.8308737874031067),
+                 (3.4536221027374268, -0.894329845905304)]
+
+        etoha_test = []
+
+        # Produce equatorial_to_horizontal_astropy results
+        for latitude, longitude, gps, ra, dec in eq:
+            result = celestial.equatorial_to_horizontal_astropy(
+                latitude, longitude, gps, [(ra, dec)])
+
+            etoha_test.extend(result)
+
+        # Check if all inputs are correct
+        np.testing.assert_almost_equal(altaz, etoha_test, 4)
+
+    @unittest.skipUnless(has_astropy, "astropy required.")
+    def test_pyephem_eqtozenaz(self):
+        """
+        celestial.equatorial_to_zenithazimuth_astropy
+        """
+
+        # This is the transform inputs
+        eq = [(-39.34633914878846, -112.2277168069694, 1295503840,
+               3.8662384455822716, -0.31222454326513827),
+              (53.13143508448587, -49.24074935964933, 985619982,
+               3.901575896592809, -0.3926720112815971),
+              (48.02031016860923, -157.4023812557098, 1126251396,
+               3.366278312183976, -1.3610394240813288)]
+        # result converted for eq->zenaz
+        zenaz = [(1.7661757628109793, -0.6043111523005624),
+                 (2.4016701141980032, 2.6012484033836625),
+                 (2.4651261727002005, -1.8828257759425302)]
+
+        eqtozenaz_test = []
+
+        # Produce equatorial_to_zenithazimuth_astropy results
+        for latitude, longitude, gps, ra, dec in eq:
+            result = celestial.equatorial_to_zenithazimuth_astropy(
+                latitude, longitude, gps, [(ra, dec)])
+            eqtozenaz_test.extend(result)
+
+        # Check if all inputs are correct, cast to numpy array for certainty
+        # Test equatorial_to_zenithazimuth_astropy
+        np.testing.assert_almost_equal(zenaz, eqtozenaz_test, 4)
+
+    @unittest.skipUnless(has_astropy, "astropy required.")
+    def test_pyephem_zenaztoeq(self):
+        """Check celestial.zenithazimuth_to_equatorial_astropy"""
+
+        # equatorial inputs from test_pyephem_eqtozenaz transformed into
+        # the shape of zenithazimuth coordinates so that they may be used for
+        # reverse benchmarking purposes.
+        zeneq = [(-39.34633914878846, -112.2277168069694,
+                  1295503840, 1.8830208700600348, -2.295442118787375),
+                 (53.13143508448587, -49.24074935964933, 985619982,
+                  1.9634683380764937, -2.3307795697979126),
+                 (48.02031016860923, -157.4023812557098, 1126251396,
+                  2.9318357508762256, -1.7954819853890793)]
+
+        # result of pyephem hor->eq/zenaz-> eq
+        efemeq = [(5.620508199785029, -0.3651173667585858),
+                  (5.244630787139936, -0.7866376569183651),
+                  (2.276751381056623, -1.0406498066785745)]
+
+        zenaztoeq_test = []
+
+        # Produce zenithazimuth_to_equatorial_astropy results
+        for latitude, longitude, gps, zen, az in zeneq:
+            result = celestial.zenithazimuth_to_equatorial_astropy(
+                latitude, longitude, gps, [(zen, az)])
+            zenaztoeq_test.extend(result)
+
+        # Check if all inputs are correct, cast to numpy array for certainty
+        # Test zenithazimuth_to_equatorial_astropy
+        np.testing.assert_almost_equal(efemeq, zenaztoeq_test, 4)
+
 
 if __name__ == '__main__':
     unittest.main()

sapphire/transformations/__init__.py

@@ -9,17 +9,19 @@
     rotation matrices and conversion between coordinate systems
 
 :mod:`~sapphire.transformations.base`
-    conversion betweeen decimal and sexagesimal
+    conversion between decimal and sexagesimal
 
 :mod:`~sapphire.transformations.celestial`
-    conversion betweeen celestial coordinate systems
+    conversion between celestial coordinate systems, contains both
+    legacy and new astropy functions
 
 :mod:`~sapphire.transformations.clock`
     conversion between different time keeping systems
 
 :mod:`~sapphire.transformations.geographic`
     geographic coordinate transformations (e.g. WGS84 to ENU)
 
+
 """
 from . import angles
 from . import axes
@@ -28,7 +30,6 @@
 from . import clock
 from . import geographic
 
-
 __all__ = ['angles',
            'axes',
            'base',