Use astropy.time.Time for epoch

python/lsst/meas/algorithms/loadIndexedReferenceObjects.py

@@ -59,7 +59,7 @@ def loadSkyCircle(self, ctrCoord, radius, filterName=None, epoch=None):
         @param[in] radius  radius of search region (an lsst.geom.Angle)
         @param[in] filterName  name of filter, or None for the default filter;
             used for flux values in case we have flux limits (which are not yet implemented)
-        @param[in] epoch  Epoch for proper motion correction (MJD TAI), or None
+        @param[in] epoch  Epoch for proper motion correction (astropy.time.Time), or None
 
         @return an lsst.pipe.base.Struct containing:
         - refCat a catalog of reference objects with the

python/lsst/meas/algorithms/loadReferenceObjects.py

@@ -25,6 +25,8 @@
 
 import abc
 
+import astropy.time
+import astropy.units
 import numpy
 
 import lsst.geom
@@ -229,7 +231,7 @@ def loadPixelBox(self, bbox, wcs, filterName=None, calib=None, epoch=None):
         @param[in] wcs  WCS (an lsst.afw.geom.SkyWcs)
         @param[in] filterName  name of camera filter, or None or blank for the default filter
         @param[in] calib  calibration, or None if unknown
-        @param[in] epoch  Epoch for proper motion correction (MJD TAI), or None
+        @param[in] epoch  Epoch for proper motion correction (astropy.time.Time), or None
 
         @return an lsst.pipe.base.Struct containing:
         - refCat a catalog of reference objects with the
@@ -268,7 +270,7 @@ def loadSkyCircle(self, ctrCoord, radius, filterName=None, epoch=None):
         @param[in] radius  radius of search region (an lsst.geom.Angle)
         @param[in] filterName  name of filter, or None for the default filter;
             used for flux values in case we have flux limits (which are not yet implemented)
-        @param[in] epoch  Epoch for proper motion correction (MJD TAI), or None
+        @param[in] epoch  Epoch for proper motion correction (astropy.time.Time), or None
 
         Note that subclasses are responsible for performing the proper motion
         correction, since this is the lowest-level interface for retrieving
@@ -519,7 +521,7 @@ def getMetadataBox(self, bbox, wcs, filterName=None, calib=None, epoch=None):
         @param[in] wcs  WCS (an lsst.afw.geom.SkyWcs)
         @param[in] filterName  name of camera filter, or None or blank for the default filter
         @param[in] calib  calibration, or None if unknown
-        @param[in] epoch  Epoch for proper motion correction (MJD TAI), or None
+        @param[in] epoch  Epoch for proper motion correction (astropy.time.Time), or None
         @return metadata (lsst.daf.base.PropertyList)
         """
         circle = self._calculateCircle(bbox, wcs)
@@ -535,7 +537,7 @@ def getMetadataCircle(self, coord, radius, filterName, calib=None, epoch=None):
         @param[in] radius  radius of circle (lsst.geom.Angle)
         @param[in] filterName  name of camera filter, or None or blank for the default filter
         @param[in] calib  calibration, or None if unknown
-        @param[in] epoch  Epoch for proper motion correction (MJD TAI), or None
+        @param[in] epoch  Epoch for proper motion correction (astropy.time.Time), or None
         @return metadata (lsst.daf.base.PropertyList)
         """
         md = PropertyList()
@@ -603,9 +605,10 @@ def applyProperMotions(self, catalog, epoch):
             - ``pm_dec`` : Proper motion in Declination (rad/yr,
                 North positive)
             - ``pm_decErr`` : Error in ``pm_dec`` (rad/yr), optional.
+            - ``epoch`` : Mean epoch of object (an astropy.time.Time)
 
-        epoch : `float`
-            Epoch to which to move objects (TAI MJD).
+        epoch : `astropy.time.Time`
+            Epoch to which to move objects.
         """
         if ("epoch" not in catalog.schema or "pm_ra" not in catalog.schema or "pm_dec" not in catalog.schema):
             if self.config.requireProperMotion:
@@ -614,16 +617,25 @@ def applyProperMotions(self, catalog, epoch):
             return
         if not catalog.isContiguous():
             raise RuntimeError("Catalog must be contiguous")
-        self.log.debug("Correcting reference catalog for proper motion to %f", epoch)
-        timeDiffsYears = (epoch - catalog["epoch"])/365.25  # Time difference, years
+        catEpoch = astropy.time.Time(catalog["epoch"], scale="tai", format="mjd")
+        self.log.debug("Correcting reference catalog for proper motion to %r", epoch)
+        # Use `epoch.tai` to make sure the time difference is in TAI
+        timeDiffsYears = (epoch.tai - catEpoch).to(astropy.units.yr).value
         coordKey = catalog.table.getCoordKey()
         # Compute the offset of each object due to proper motion
         # as components of the arc of a great circle along RA and Dec
-        offsetsRaRad = catalog["pm_ra"]*timeDiffsYears
-        offsetsDecRad = catalog["pm_dec"]*timeDiffsYears
+        pmRaRad = catalog["pm_ra"]
+        pmDecRad = catalog["pm_dec"]
+        offsetsRaRad = pmRaRad*timeDiffsYears
+        offsetsDecRad = pmDecRad*timeDiffsYears
         # Compute the corresponding bearing and arc length of each offset
         # due to proper motion, and apply the offset
-        offsetBearingsRad = numpy.arctan2(offsetsDecRad, offsetsRaRad)
+        # The factor of 1e6 for computing bearing is intended as
+        # a reasonable scale for typical values of proper motion
+        # in order to avoid large errors for small values of proper motion;
+        # using the offsets is another option, but it can give
+        # needlessly large errors for short duration
+        offsetBearingsRad = numpy.arctan2(pmDecRad*1e6, pmRaRad*1e6)
         offsetAmountsRad = numpy.hypot(offsetsRaRad, offsetsDecRad)
         for record, bearingRad, amountRad in zip(catalog, offsetBearingsRad, offsetAmountsRad):
             record.set(coordKey,