Add utility to compute apparent motion

python/lsst/drp/tasks/gbdesAstrometricFit.py

@@ -44,6 +44,7 @@
 from smatch.matcher import Matcher
 
 __all__ = [
+    "calculate_apparent_motion",
     "GbdesAstrometricFitConnections",
     "GbdesAstrometricFitConfig",
     "GbdesAstrometricFitTask",
@@ -57,6 +58,66 @@
 ]
 
 
+def calculate_apparent_motion(data, refEpoch):
+    """Calculate shift from reference epoch to the apparent observed position
+    at another date.
+
+    This function calculates the shift due to proper motion combined with the
+    apparent motion due to parallax. This is not used in the
+    `GbdesAstrometricFitTask` or related child tasks, but is useful for
+    assessing results.
+
+    Parameters
+    ----------
+    data : `pd.DataFrame`
+        Table containing position, proper motion, parallax, and epoch for each
+        source.
+    refEpoch : `float`
+        Epoch of the reference position.
+
+    Returns
+    -------
+    apparentMotionRA : `np.ndarray` [`astropy.units.Quantity`]
+        RA shift in degrees.
+    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)
+
+    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)
+    sunLongitudes = sun.transform_to(frame).lon.radian
+
+    # These equations for parallax come from Equations 5.2 in Van de Kamp's
+    # book Stellar Paths. They differ from the parallax calculated in gbdes by
+    # ~0.01 mas, which is acceptable for QA and plotting purposes.
+    parallaxFactorRA = np.cos(wcsfit.EclipticInclination) * np.cos(ra_rad) * np.sin(sunLongitudes) - np.sin(
+        ra_rad
+    ) * np.cos(sunLongitudes)
+    parallaxFactorDec = (
+        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)
+    parallaxCorrectionRA = parallaxDegrees * parallaxFactorRA
+    parallaxCorrectionDec = parallaxDegrees * parallaxFactorDec
+
+    apparentMotionRA = properMotionRA + parallaxCorrectionRA
+    apparentMotionDec = properMotionDec + parallaxCorrectionDec
+
+    return apparentMotionRA, apparentMotionDec
+
+
 def _make_ref_covariance_matrix(
     refCat, inputUnit=u.radian, outputCoordUnit=u.marcsec, outputPMUnit=u.marcsec, version=1
 ):

tests/test_gbdesAstrometricFit.py

@@ -32,13 +32,16 @@
 import pandas as pd
 import wcsfit
 import yaml
+from astropy.coordinates import Distance, SkyCoord
+from astropy.time import Time
 from lsst import sphgeom
 from lsst.daf.base import PropertyList
 from lsst.drp.tasks.gbdesAstrometricFit import (
     GbdesAstrometricFitConfig,
     GbdesAstrometricFitTask,
     GbdesGlobalAstrometricFitConfig,
     GbdesGlobalAstrometricFitTask,
+    calculate_apparent_motion,
 )
 from lsst.meas.algorithms import ReferenceObjectLoader
 from lsst.meas.algorithms.testUtils import MockRefcatDataId
@@ -1058,6 +1061,114 @@ def test_inputCameraModel(self):
                     self.assertAlmostEqual(np.std(dDec), 0)
 
 
+class TestApparentMotion(lsst.utils.tests.TestCase):
+    """This class simulates an array of objects with random proper motions and
+    parallaxes and compares the results of
+    `lsst.drp.tasks.gbdesAstrometricFit.calculate_apparent_motion` against the
+    results calculated by astropy."""
+
+    @classmethod
+    def setUpClass(cls):
+        np.random.seed(12345)
+
+        cls.RAs = np.random.rand(10) + 150.0
+        cls.Decs = np.random.rand(10) + 2.5
+
+        cls.pmRAs = np.random.random(10) * 10
+        cls.pmDecs = np.random.random(10) * 10
+
+        cls.parallaxes = abs(np.random.random(10) * 5)
+
+        cls.mjds = np.random.rand(10) * 20 + 57000
+
+    def testProperMotion(self):
+        """Calculate the change in position due to proper motion for a given
+        time change, and compare the results against astropy.
+        """
+
+        parallaxes = np.zeros(len(self.RAs))
+        refEpoch = 57100
+
+        data = pd.DataFrame(
+            {
+                "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)
+
+        pmRA_cosDec = self.pmRAs * np.cos((self.Decs * u.degree).to(u.radian))
+        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"),
+        )
+
+        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
+
+        np.testing.assert_allclose(astropy_dRA, ra_correction.value, rtol=1e-6)
+        np.testing.assert_allclose(astropy_dDec, dec_correction.value, rtol=1e-6)
+
+    def testParallax(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
+        annual aberration, but it can be obtained by comparing the calculated
+        positions of an object with a given parallax and one with a much
+        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
+
+        data = pd.DataFrame(
+            {
+                "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)
+
+        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"),
+        )
+        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"),
+        )
+
+        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
+
+        np.testing.assert_allclose(astropy_dRA, ra_correction.value, rtol=1e-1)
+        np.testing.assert_allclose(astropy_dDec, dec_correction.value, rtol=1e-1)
+
+
 def setup_module(module):
     lsst.utils.tests.init()