From 26b0488b28ef13071e60527f455d5633d7d01371 Mon Sep 17 00:00:00 2001 From: Clare Saunders Date: Fri, 2 Aug 2024 13:53:42 -0700 Subject: [PATCH] Add utility to compute apparent motion --- python/lsst/drp/tasks/gbdesAstrometricFit.py | 57 +++++++++ tests/test_gbdesAstrometricFit.py | 117 +++++++++++++++++++ 2 files changed, 174 insertions(+) diff --git a/python/lsst/drp/tasks/gbdesAstrometricFit.py b/python/lsst/drp/tasks/gbdesAstrometricFit.py index d13fb3a3..fac745d2 100644 --- a/python/lsst/drp/tasks/gbdesAstrometricFit.py +++ b/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,62 @@ ] +def calculate_apparent_motion(catalog, 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 + ---------- + catalog : `astropy.table.Table` + Table containing position, proper motion, parallax, and epoch for each + source, labeled by columns 'ra', 'dec', 'pmRA', 'pmDec', 'parallax', + and 'MJD'. + refEpoch : `astropy.time.Time` + 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 = catalog["ra"].to(u.rad) + dec_rad = catalog["dec"].to(u.rad) + + 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 + # 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 = catalog["parallax"].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 ): diff --git a/tests/test_gbdesAstrometricFit.py b/tests/test_gbdesAstrometricFit.py index 66a78665..a69e354d 100644 --- a/tests/test_gbdesAstrometricFit.py +++ b/tests/test_gbdesAstrometricFit.py @@ -23,6 +23,8 @@ import unittest from copy import copy +import astropy.table +import astropy.time import astropy.units as u import lsst.afw.geom as afwgeom import lsst.afw.table as afwTable @@ -32,6 +34,7 @@ import pandas as pd import wcsfit import yaml +from astropy.coordinates import Distance, SkyCoord from lsst import sphgeom from lsst.daf.base import PropertyList from lsst.drp.tasks.gbdesAstrometricFit import ( @@ -39,6 +42,7 @@ GbdesAstrometricFitTask, GbdesGlobalAstrometricFitConfig, GbdesGlobalAstrometricFitTask, + calculate_apparent_motion, ) from lsst.meas.algorithms import ReferenceObjectLoader from lsst.meas.algorithms.testUtils import MockRefcatDataId @@ -1136,6 +1140,119 @@ def test_useColor(self): self.assertEqual(len(outputs.colorParams["visit"]), len(self.testVisits)) +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.""" + + def setUp(self): + np.random.seed(12345) + + self.ras = (np.random.rand(10) + 150.0) * u.degree + self.decs = (np.random.rand(10) + 2.5) * u.degree + + self.pmRAs = (np.random.random(10) * 10) * u.mas / u.yr + self.pmDecs = (np.random.random(10) * 10) * u.mas / u.yr + + self.parallaxes = (abs(np.random.random(10) * 5)) * u.mas + + self.mjds = astropy.time.Time(np.random.rand(10) * 20 + 57000, format="mjd", scale="tai") + + self.refEpoch = astropy.time.Time(57100, format="mjd", scale="tai") + + 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 = astropy.table.Table( + { + "ra": self.ras, + "dec": self.decs, + "pmRA": self.pmRAs, + "pmDec": self.pmDecs, + "parallax": parallaxes, + "MJD": self.mjds, + } + ) + raCorrection, decCorrection = calculate_apparent_motion(data, self.refEpoch) + + pmRACosDec = self.pmRAs * np.cos(self.decs) + coords = SkyCoord( + ra=self.ras, + dec=self.decs, + pm_ra_cosdec=pmRACosDec, + pm_dec=self.pmDecs, + obstime=self.refEpoch, + ) + + newCoords = coords.apply_space_motion(new_obstime=self.mjds) + astropyDRA = newCoords.ra.degree - coords.ra.degree + astropyDDec = newCoords.dec.degree - coords.dec.degree + + # The proper motion values are on the order of 0-3 mas/yr and the + # differences are on the scale of 1e-8 mas. + self.assertFloatsAlmostEqual(astropyDRA, raCorrection.value, rtol=1e-6) + self.assertFloatsAlmostEqual(astropyDDec, decCorrection.value, rtol=1e-6) + + 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 + 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. + """ + # Set proper motions to zero to test parallax correction by itself. + pms = np.zeros(len(self.ras)) * u.mas / u.yr + + data = astropy.table.Table( + { + "ra": self.ras, + "dec": self.decs, + "pmRA": pms, + "pmDec": pms, + "parallax": self.parallaxes, + "MJD": self.mjds, + } + ) + raCorrection, decCorrection = calculate_apparent_motion(data, self.refEpoch) + + coords = SkyCoord( + 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, + 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=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 + + # The parallax values are on the scale of 1-3 mas, and the difference + # between the two calculation methods is ~0.05 mas. + self.assertFloatsAlmostEqual(astropyDRA, raCorrection.value, rtol=2e-2) + self.assertFloatsAlmostEqual(astropyDDec, decCorrection.value, rtol=2e-2) + + def setup_module(module): lsst.utils.tests.init()