Switch to astropy table

python/lsst/drp/tasks/gbdesAstrometricFit.py

@@ -58,7 +58,7 @@
 ]
 
 
-def calculate_apparent_motion(data, refEpoch):
+def calculate_apparent_motion(catalog, refEpoch):
     """Calculate shift from reference epoch to the apparent observed position
     at another date.
 
@@ -69,10 +69,11 @@ def calculate_apparent_motion(data, refEpoch):
 
     Parameters
     ----------
-    data : `pd.DataFrame`
+    catalog : `astropy.table.Table`
         Table containing position, proper motion, parallax, and epoch for each
-        source.
-    refEpoch : `float`
+        source, labeled by columns 'ra', 'dec', 'pmRA', 'pmDec', 'parallax',
+        and 'MJD'.
+    refEpoch : `astropy.time.Time`
         Epoch of the reference position.
 
     Returns
@@ -82,20 +83,15 @@ def calculate_apparent_motion(data, refEpoch):
     apparentMotionDec : `np.ndarray` [`astropy.units.Quantity`]
         Dec shift in degrees.
     """
-    ra_rad = (data["ra"].values * u.deg).to(u.rad)
-    dec_rad = (data["dec"].values * u.deg).to(u.rad)
+    ra_rad = catalog["ra"].to(u.rad)
+    dec_rad = catalog["dec"].to(u.rad)
 
-    dt = (astropy.time.Time(data["MJD"].values, format="mjd") - astropy.time.Time(refEpoch, format="mjd")).to(
-        u.yr
-    )
-    pmRA_deg = (data["pmRA"].values * u.mas).to(u.deg) / u.yr
-    pmDec_deg = (data["pmDec"].values * u.mas).to(u.deg) / u.yr
-    properMotionRA = pmRA_deg * dt
-    properMotionDec = pmDec_deg * dt
-
-    obsTimes = astropy.time.Time(data["MJD"], format="mjd")
-    sun = astropy.coordinates.get_body("sun", time=obsTimes)
-    frame = astropy.coordinates.GeocentricTrueEcliptic(equinox=obsTimes)
+    dt = (catalog["MJD"] - refEpoch).to(u.yr)
+    properMotionRA = catalog["pmRA"].to(u.deg / u.yr) * dt
+    properMotionDec = catalog["pmDec"].to(u.deg / u.yr) * dt
+
+    sun = astropy.coordinates.get_body("sun", time=catalog["MJD"])
+    frame = astropy.coordinates.GeocentricTrueEcliptic(equinox=catalog["MJD"])
     sunLongitudes = sun.transform_to(frame).lon.radian
 
     # These equations for parallax come from Equations 5.2 in Van de Kamp's
@@ -108,7 +104,7 @@ def calculate_apparent_motion(data, refEpoch):
         np.sin(wcsfit.EclipticInclination) * np.cos(dec_rad)
         - np.cos(wcsfit.EclipticInclination) * np.sin(ra_rad) * np.sin(dec_rad)
     ) * np.sin(sunLongitudes) - np.cos(ra_rad) * np.sin(dec_rad) * np.cos(sunLongitudes)
-    parallaxDegrees = (data["parallax"].to_numpy() * u.mas).to(u.degree)
+    parallaxDegrees = catalog["parallax"].to(u.degree)
     parallaxCorrectionRA = parallaxDegrees * parallaxFactorRA
     parallaxCorrectionDec = parallaxDegrees * parallaxFactorDec
 

tests/test_gbdesAstrometricFit.py

@@ -33,7 +33,8 @@
 import wcsfit
 import yaml
 from astropy.coordinates import Distance, SkyCoord
-from astropy.time import Time
+from astropy.table import Table
+import astropy.time
 from lsst import sphgeom
 from lsst.daf.base import PropertyList
 from lsst.drp.tasks.gbdesAstrometricFit import (
@@ -1149,53 +1150,54 @@ class TestApparentMotion(lsst.utils.tests.TestCase):
     def setUpClass(cls):
         np.random.seed(12345)
 
-        cls.RAs = np.random.rand(10) + 150.0
-        cls.Decs = np.random.rand(10) + 2.5
+        cls.ras = (np.random.rand(10) + 150.0) * u.degree
+        cls.decs = (np.random.rand(10) + 2.5) * u.degree
 
-        cls.pmRAs = np.random.random(10) * 10
-        cls.pmDecs = np.random.random(10) * 10
+        cls.pmRAs = (np.random.random(10) * 10) * u.mas / u.yr
+        cls.pmDecs = (np.random.random(10) * 10) * u.mas / u.yr
 
-        cls.parallaxes = abs(np.random.random(10) * 5)
+        cls.parallaxes = (abs(np.random.random(10) * 5)) * u.mas
 
-        cls.mjds = np.random.rand(10) * 20 + 57000
+        cls.mjds = astropy.time.Time(np.random.rand(10) * 20 + 57000, format="mjd", scale="tai")
 
-    def testProperMotion(self):
-        """Calculate the change in position due to proper motion for a given
-        time change, and compare the results against astropy.
-        """
+        cls.refEpoch = astropy.time.Time(57100, format="mjd", scale="tai")
 
-        parallaxes = np.zeros(len(self.RAs))
-        refEpoch = 57100
+    def test_proper_motion(self):
+        """Calculate the change in position due to proper motion only for a
+        given time change, and compare the results against astropy.
+        """
+        # Set parallaxes to zero to test proper motion part by itself.
+        parallaxes = np.zeros(len(self.ras)) * u.mas
 
-        data = pd.DataFrame(
+        data = Table(
             {
-                "ra": self.RAs,
-                "dec": self.Decs,
+                "ra": self.ras,
+                "dec": self.decs,
                 "pmRA": self.pmRAs,
                 "pmDec": self.pmDecs,
                 "parallax": parallaxes,
                 "MJD": self.mjds,
             }
         )
-        ra_correction, dec_correction = calculate_apparent_motion(data, refEpoch)
+        raCorrection, decCorrection = calculate_apparent_motion(data, self.refEpoch)
 
-        pmRA_cosDec = self.pmRAs * np.cos((self.Decs * u.degree).to(u.radian))
+        pmRACosDec = self.pmRAs * np.cos(self.decs)
         coords = SkyCoord(
-            ra=self.RAs * u.degree,
-            dec=self.Decs * u.degree,
-            pm_ra_cosdec=pmRA_cosDec * u.mas / u.yr,
-            pm_dec=self.pmDecs * u.mas / u.yr,
-            obstime=Time(refEpoch, format="mjd"),
+            ra=self.ras,
+            dec=self.decs,
+            pm_ra_cosdec=pmRACosDec,
+            pm_dec=self.pmDecs,
+            obstime=self.refEpoch,
         )
 
-        newCoords = coords.apply_space_motion(new_obstime=Time(self.mjds, format="mjd"))
-        astropy_dRA = newCoords.ra.degree - coords.ra.degree
-        astropy_dDec = newCoords.dec.degree - coords.dec.degree
+        newCoords = coords.apply_space_motion(new_obstime=self.mjds)
+        astropyDRA = newCoords.ra.degree - coords.ra.degree
+        astropyDDec = newCoords.dec.degree - coords.dec.degree
 
-        np.testing.assert_allclose(astropy_dRA, ra_correction.value, rtol=1e-6)
-        np.testing.assert_allclose(astropy_dDec, dec_correction.value, rtol=1e-6)
+        self.assertFloatsAlmostEqual(astropyDRA, raCorrection.value, rtol=1e-6)
+        self.assertFloatsAlmostEqual(astropyDDec, decCorrection.value, rtol=1e-6)
 
-    def testParallax(self):
+    def test_parallax(self):
         """Calculate the change in position due to parallax for a given time
         change and compare the results against astropy. Astropy will not give
         the parallax part separate from other sources of space motion, such as
@@ -1204,47 +1206,48 @@ def testParallax(self):
         further distance. The result is still only approximately the same, but
         is close enough for accuracy needed here.
         """
-        pms = np.zeros(len(self.RAs))
-        refEpoch = 57100
+        # Set proper motions to zero to test parallax correction by itself.
+        pms = np.zeros(len(self.ras)) * u.mas / u.yr
 
-        data = pd.DataFrame(
+        data = Table(
             {
-                "ra": self.RAs,
-                "dec": self.Decs,
+                "ra": self.ras,
+                "dec": self.decs,
                 "pmRA": pms,
                 "pmDec": pms,
                 "parallax": self.parallaxes,
                 "MJD": self.mjds,
             }
         )
-        ra_correction, dec_correction = calculate_apparent_motion(data, refEpoch)
+        raCorrection, decCorrection = calculate_apparent_motion(data, self.refEpoch)
 
         coords = SkyCoord(
-            ra=self.RAs * u.degree,
-            dec=self.Decs * u.degree,
-            pm_ra_cosdec=pms * u.mas / u.yr,
-            pm_dec=pms * u.mas / u.yr,
-            distance=Distance(parallax=self.parallaxes * u.mas),
-            obstime=Time(refEpoch, format="mjd"),
+            ra=self.ras,
+            dec=self.decs,
+            pm_ra_cosdec=pms,
+            pm_dec=pms,
+            distance=Distance(parallax=self.parallaxes),
+            obstime=self.refEpoch,
         )
+        # Make another set of objects at the same coordinates but a much
+        # greater distant in order to isolate parallax when applying space
+        # motion.
         refCoords = SkyCoord(
-            ra=self.RAs * u.degree,
-            dec=self.Decs * u.degree,
-            pm_ra_cosdec=pms * u.mas / u.yr,
-            pm_dec=pms * u.mas / u.yr,
-            distance=1e10 * Distance(parallax=self.parallaxes * u.mas),
-            obstime=Time(refEpoch, format="mjd"),
+            ra=self.ras,
+            dec=self.decs,
+            pm_ra_cosdec=pms,
+            pm_dec=pms,
+            distance=1e10 * Distance(parallax=self.parallaxes),
+            obstime=self.refEpoch,
         )
 
-        newCoords = coords.apply_space_motion(new_obstime=Time(self.mjds, format="mjd")).transform_to("gcrs")
-        refNewCoords = refCoords.apply_space_motion(new_obstime=Time(self.mjds, format="mjd")).transform_to(
-            "gcrs"
-        )
-        astropy_dRA = newCoords.ra.degree - refNewCoords.ra.degree
-        astropy_dDec = newCoords.dec.degree - refNewCoords.dec.degree
+        newCoords = coords.apply_space_motion(new_obstime=self.mjds).transform_to("gcrs")
+        refNewCoords = refCoords.apply_space_motion(new_obstime=self.mjds).transform_to("gcrs")
+        astropyDRA = newCoords.ra.degree - refNewCoords.ra.degree
+        astropyDDec = newCoords.dec.degree - refNewCoords.dec.degree
 
-        np.testing.assert_allclose(astropy_dRA, ra_correction.value, rtol=1e-1)
-        np.testing.assert_allclose(astropy_dDec, dec_correction.value, rtol=1e-1)
+        self.assertFloatsAlmostEqual(astropyDRA, raCorrection.value, rtol=1e-1)
+        self.assertFloatsAlmostEqual(astropyDDec, decCorrection.value, rtol=1e-1)
 
 
 def setup_module(module):