DM-24277: Apply proper motion and parallax while loading refcats in Jointcal

python/lsst/jointcal/ccdImage.cc

@@ -56,6 +56,9 @@ void declareCcdImage(py::module &mod) {
     cls.def("getCcdId", &CcdImage::getCcdId);
     cls.def_property_readonly("ccdId", &CcdImage::getCcdId);
 
+    cls.def("getMjd", &CcdImage::getMjd);
+    cls.def_property_readonly("mjd", &CcdImage::getMjd);
+
     cls.def("getImageFrame", &CcdImage::getImageFrame, py::return_value_policy::reference_internal);
     cls.def_property_readonly("imageFrame", &CcdImage::getImageFrame,
                               py::return_value_policy::reference_internal);

python/lsst/jointcal/jointcal.py

@@ -22,6 +22,7 @@
 import collections
 import os
 
+import astropy.time
 import numpy as np
 import astropy.units as u
 
@@ -370,9 +371,10 @@ def setDefaults(self):
                     'base_PsfFlux_flag', 'base_PixelFlags_flag_suspectCenter']
         self.sourceSelector['science'].flags.bad = badFlags
 
-        # Default to Gaia-DR2 for astrometry and PS1-DR1 for photometry,
-        # with a reasonable initial filterMap.
+        # Default to Gaia-DR2 (with proper motions) for astrometry and
+        # PS1-DR1 for photometry, with a reasonable initial filterMap.
         self.astrometryRefObjLoader.ref_dataset_name = "gaia_dr2_20200414"
+        self.astrometryRefObjLoader.requireProperMotion = True
         self.astrometryRefObjLoader.filterMap = {'u': 'phot_g_mean',
                                                  'g': 'phot_g_mean',
                                                  'r': 'phot_g_mean',
@@ -654,6 +656,22 @@ def _get_refcat_coordinate_error_override(self, refCat, name):
         else:
             return self.config.astrometryReferenceErr
 
+    def _compute_proper_motion_epoch(self, ccdImageList):
+        """Return the proper motion correction epoch of the provided images.
+
+        Parameters
+        ----------
+        ccdImageList : `list` [`lsst.jointcal.CcdImage`]
+            The images to compute the appropriate epoch for.
+
+        Returns
+        -------
+        epoch : `astropy.time.Time`
+            The date to use for proper motion corrections.
+        """
+        mjds = [ccdImage.getMjd() for ccdImage in ccdImageList]
+        return astropy.time.Time(np.mean(mjds), format='mjd', scale="tai")
+
     def _do_load_refcat_and_fit(self, associations, defaultFilter, center, radius,
                                 filters=[],
                                 tract="", profile_jointcal=False, match_cut=3.0,
@@ -705,9 +723,11 @@ def _do_load_refcat_and_fit(self, associations, defaultFilter, center, radius,
                         associations.fittedStarListSize())
 
         applyColorterms = False if name.lower() == "astrometry" else self.config.applyColorTerms
+        epoch = self._compute_proper_motion_epoch(associations.getCcdImageList())
         refCat, fluxField = self._load_reference_catalog(refObjLoader, referenceSelector,
                                                          center, radius, defaultFilter,
-                                                         applyColorterms=applyColorterms)
+                                                         applyColorterms=applyColorterms,
+                                                         epoch=epoch)
         refCoordErr = self._get_refcat_coordinate_error_override(refCat, name)
 
         associations.collectRefStars(refCat,
@@ -742,7 +762,7 @@ def _do_load_refcat_and_fit(self, associations, defaultFilter, center, radius,
         return result
 
     def _load_reference_catalog(self, refObjLoader, referenceSelector, center, radius, filterName,
-                                applyColorterms=False):
+                                applyColorterms=False, epoch=None):
         """Load the necessary reference catalog sources, convert fluxes to
         correct units, and apply color term corrections if requested.
 
@@ -760,6 +780,9 @@ def _load_reference_catalog(self, refObjLoader, referenceSelector, center, radiu
             The name of the camera filter to load fluxes for.
         applyColorterms : `bool`
             Apply colorterm corrections to the refcat for ``filterName``?
+        epoch : `astropy.time.Time`, optional
+            Epoch to which to correct refcat proper motion and parallax,
+            or `None` to not apply such corrections.
 
         Returns
         -------
@@ -770,7 +793,8 @@ def _load_reference_catalog(self, refObjLoader, referenceSelector, center, radiu
         """
         skyCircle = refObjLoader.loadSkyCircle(center,
                                                radius,
-                                               filterName)
+                                               filterName,
+                                               epoch=epoch)
 
         selected = referenceSelector.run(skyCircle.refCat)
         # Need memory contiguity to get reference filters as a vector.
@@ -780,11 +804,7 @@ def _load_reference_catalog(self, refObjLoader, referenceSelector, center, radiu
             refCat = selected.sourceCat
 
         if applyColorterms:
-            try:
-                refCatName = refObjLoader.ref_dataset_name
-            except AttributeError:
-                # NOTE: we need this try:except: block in place until we've completely removed a.net support.
-                raise RuntimeError("Cannot perform colorterm corrections with a.net refcats.")
+            refCatName = refObjLoader.ref_dataset_name
             self.log.info("Applying color terms for filterName=%r reference catalog=%s",
                           filterName, refCatName)
             colorterm = self.config.colorterms.getColorterm(

tests/test_jointcal.py

@@ -19,6 +19,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
+import itertools
 import unittest
 from unittest import mock
 
@@ -46,15 +47,35 @@ def setup_module(module):
 
 def make_fake_refcat(center, flux, filterName):
     """Make a fake reference catalog."""
-    schema = LoadIndexedReferenceObjectsTask.makeMinimalSchema([filterName])
+    schema = LoadIndexedReferenceObjectsTask.makeMinimalSchema([filterName],
+                                                               addProperMotion=True)
     catalog = lsst.afw.table.SimpleCatalog(schema)
     record = catalog.addNew()
     record.setCoord(center)
     record[filterName + '_flux'] = flux
     record[filterName + '_fluxErr'] = flux*0.1
+    record['pm_ra'] = lsst.geom.Angle(1)
+    record['pm_dec'] = lsst.geom.Angle(2)
+    record['epoch'] = 65432.1
     return catalog
 
 
+def make_fake_wcs():
+    """Return two simple SkyWcs objects, with slightly different sky positions.
+
+    Use the same pixel origins as the cfht_minimal data, but put the sky origin
+    at RA=0
+    """
+    crpix = lsst.geom.Point2D(931.517869, 2438.572109)
+    cd = np.array([[5.19513851e-05, -2.81124812e-07],
+                   [-3.25186974e-07, -5.19112119e-05]])
+    crval1 = lsst.geom.SpherePoint(0.01, -0.01, lsst.geom.degrees)
+    crval2 = lsst.geom.SpherePoint(-0.01, 0.01, lsst.geom.degrees)
+    wcs1 = lsst.afw.geom.makeSkyWcs(crpix, crval1, cd)
+    wcs2 = lsst.afw.geom.makeSkyWcs(crpix, crval2, cd)
+    return wcs1, wcs2
+
+
 class JointcalTestBase:
     def setUp(self):
         struct = lsst.jointcal.testUtils.createTwoFakeCcdImages(100, 100)
@@ -174,7 +195,7 @@ def test_writeChi2(self):
 class TestJointcalLoadRefCat(JointcalTestBase, lsst.utils.tests.TestCase):
 
     def _make_fake_refcat(self):
-        """Make a fake reference catalog and the bits necessary to use it."""
+        """Mock a fake reference catalog and the bits necessary to use it."""
         center = lsst.geom.SpherePoint(30, -30, lsst.geom.degrees)
         flux = 10
         radius = 1 * lsst.geom.degrees
@@ -195,6 +216,8 @@ def test_load_reference_catalog(self):
         config.astrometryReferenceErr = 0.1  # our test refcats don't have coord errors
         jointcal = lsst.jointcal.JointcalTask(config=config, butler=self.butler)
 
+        # NOTE: we cannot test application of proper motion here, because we
+        # mock the refObjLoader, so the real loader is never called.
         refCat, fluxField = jointcal._load_reference_catalog(refObjLoader,
                                                              jointcal.astrometryReferenceSelector,
                                                              center,
@@ -325,14 +348,7 @@ def testPointsRA0(self):
         """Test for CcdImages crossing RA=0; this demonstrates a fix for
         the bug described in DM-19802.
         """
-        # Use the same pixel origins as the cfht_minimal data, but put the sky origin at RA=0
-        crpix = lsst.geom.Point2D(931.517869, 2438.572109)
-        cd = np.array([[5.19513851e-05, -2.81124812e-07],
-                      [-3.25186974e-07, -5.19112119e-05]])
-        crval1 = lsst.geom.SpherePoint(0.01, -0.01, lsst.geom.degrees)
-        crval2 = lsst.geom.SpherePoint(-0.01, 0.01, lsst.geom.degrees)
-        wcs1 = lsst.afw.geom.makeSkyWcs(crpix, crval1, cd)
-        wcs2 = lsst.afw.geom.makeSkyWcs(crpix, crval2, cd)
+        wcs1, wcs2 = make_fake_wcs()
 
         # Put the visit boresights at the WCS origin, for consistency
         visitInfo1 = lsst.afw.image.VisitInfo(exposureId=30577512,
@@ -348,6 +364,40 @@ def testPointsRA0(self):
                          struct.skyWcs[0], struct.skyWcs[1], struct.bbox)
 
 
+class TestJointcalComputePMDate(JointcalTestBase, lsst.utils.tests.TestCase):
+    """Tests of jointcal._compute_proper_motion_epoch()"""
+    def test_compute_proper_motion_epoch(self):
+        mjds = np.array((65432.1, 66666, 65555, 64322.2))
+
+        wcs1, wcs2 = make_fake_wcs()
+        visitInfo1 = lsst.afw.image.VisitInfo(exposureId=30577512,
+                                              date=DateTime(mjds[0]),
+                                              boresightRaDec=wcs1.getSkyOrigin())
+        visitInfo2 = lsst.afw.image.VisitInfo(exposureId=30621144,
+                                              date=DateTime(mjds[1]),
+                                              boresightRaDec=wcs2.getSkyOrigin())
+        visitInfo3 = lsst.afw.image.VisitInfo(exposureId=30577513,
+                                              date=DateTime(mjds[2]),
+                                              boresightRaDec=wcs1.getSkyOrigin())
+        visitInfo4 = lsst.afw.image.VisitInfo(exposureId=30621145,
+                                              date=DateTime(mjds[3]),
+                                              boresightRaDec=wcs2.getSkyOrigin())
+
+        struct1 = lsst.jointcal.testUtils.createTwoFakeCcdImages(fakeWcses=[wcs1, wcs2],
+                                                                 fakeVisitInfos=[visitInfo1, visitInfo2])
+        struct2 = lsst.jointcal.testUtils.createTwoFakeCcdImages(fakeWcses=[wcs1, wcs2],
+                                                                 fakeVisitInfos=[visitInfo3, visitInfo4])
+        ccdImageList = list(itertools.chain(struct1.ccdImageList, struct2.ccdImageList))
+        associations = lsst.jointcal.Associations()
+        for ccdImage in ccdImageList:
+            associations.addCcdImage(ccdImage)
+        associations.computeCommonTangentPoint()
+
+        jointcal = lsst.jointcal.JointcalTask(config=self.config, butler=self.butler)
+        result = jointcal._compute_proper_motion_epoch(ccdImageList)
+        self.assertEqual(result.mjd, mjds.mean())
+
+
 class MemoryTester(lsst.utils.tests.MemoryTestCase):
     pass
 

tests/test_jointcal_cfht.py

@@ -52,10 +52,10 @@ def setUpClass(cls):
             raise unittest.SkipTest("obs_cfht not setup")
 
     def setUp(self):
-        # We don't want the absolute astrometry to become significantly worse
-        # than the single-epoch astrometry (about 0.040").
+        # NOTE: refcat-comparison RMS error is worse now, because the
+        # comparison code is not applying the proper motion data.
         # See Readme for an explanation of this empirical value.
-        self.dist_rms_absolute = 56e-3*u.arcsecond
+        self.dist_rms_absolute = 70e-3*u.arcsecond
 
         do_plot = False
 
@@ -99,8 +99,8 @@ def test_jointcalTask_2_visits(self):
                        'selected_photometry_fittedStars': 2232,
                        'selected_astrometry_ccdImages': 12,
                        'selected_photometry_ccdImages': 12,
-                       'astrometry_final_chi2': 1243.59,
-                       'astrometry_final_ndof': 2446,
+                       'astrometry_final_chi2': 1146.81,
+                       'astrometry_final_ndof': 2486,
                        'photometry_final_chi2': 11624.3,
                        'photometry_final_ndof': 2778
                        }
@@ -137,8 +137,8 @@ def setup_jointcalTask_2_visits_constrainedAstrometry(self):
                    'associated_astrometry_fittedStars': 2272,
                    'selected_astrometry_fittedStars': 1229,
                    'selected_astrometry_ccdImages': 12,
-                   'astrometry_final_chi2': 1320.64,
-                   'astrometry_final_ndof': 2532,
+                   'astrometry_final_chi2': 1204.72,
+                   'astrometry_final_ndof': 2576,
                    }
 
         return dist_rms_relative, metrics
@@ -147,7 +147,6 @@ def test_jointcalTask_2_visits_constrainedAstrometry_no_photometry(self):
         dist_rms_relative, metrics = self.setup_jointcalTask_2_visits_constrainedAstrometry()
         self.config.writeInitialModel = True  # write the initial models
         # use a temporary directory for debug output, to prevent test collisions
-        # use a temporary directory for debug output, to prevent test collisions
         with tempfile.TemporaryDirectory() as tempdir:
             self.config.debugOutputPath = tempdir
 
@@ -169,8 +168,8 @@ def test_jointcalTask_2_visits_constrainedAstrometry_4sigma_outliers(self):
         """
         dist_rms_relative, metrics = self.setup_jointcalTask_2_visits_constrainedAstrometry()
         self.config.outlierRejectSigma = 4
-        metrics['astrometry_final_chi2'] = 1006.02
-        metrics['astrometry_final_ndof'] = 2386
+        metrics['astrometry_final_chi2'] = 849.75
+        metrics['astrometry_final_ndof'] = 2390
 
         self._testJointcalTask(2, dist_rms_relative, self.dist_rms_absolute, None, metrics=metrics)
 
@@ -180,8 +179,8 @@ def test_jointcalTask_2_visits_constrainedAstrometry_astrometryReferenceUncertai
         test_config = os.path.join(lsst.utils.getPackageDir('jointcal'),
                                    'tests/config/astrometryReferenceErr-config.py')
         self.configfiles.append(test_config)
-        metrics['astrometry_final_chi2'] = 1456.09
-        metrics['astrometry_final_ndof'] = 2170
+        metrics['astrometry_final_chi2'] = 1275.70
+        metrics['astrometry_final_ndof'] = 2062
 
         self._testJointcalTask(2, dist_rms_relative, self.dist_rms_absolute, None, metrics=metrics)
 

tests/test_jointcal_decam.py

@@ -96,8 +96,8 @@ def test_jointcalTask_2_visits(self):
                    'selected_photometry_fittedStars': 4637,
                    'selected_astrometry_ccdImages': 15,
                    'selected_photometry_ccdImages': 15,
-                   'astrometry_final_chi2': 1.15399,
-                   'astrometry_final_ndof': 518,
+                   'astrometry_final_chi2': 1.21494,
+                   'astrometry_final_ndof': 560,
                    'photometry_final_chi2': 16529,
                    'photometry_final_ndof': 3634,
                    }

tests/test_jointcal_hsc.py

@@ -107,7 +107,7 @@ def test_jointcalTask_10_visits_simple_astrometry_no_photometry(self):
         self.jointcalStatistics.do_photometry = False
 
         # See Readme for an explanation of these empirical values.
-        dist_rms_absolute = 22e-3*u.arcsecond
+        dist_rms_absolute = 23e-3*u.arcsecond
         dist_rms_relative = 13e-3*u.arcsecond
         pa1 = None
         metrics = {'collected_astrometry_refStars': 1316,
@@ -242,7 +242,7 @@ def test_jointcalTask_3_visits_simple_astrometry_no_photometry(self):
                    'associated_astrometry_fittedStars': 3199,
                    'selected_astrometry_fittedStars': 1282,
                    'selected_astrometry_ccdImages': 9,
-                   'astrometry_final_chi2': 1013.22,
+                   'astrometry_final_chi2': 1013.20,
                    'astrometry_final_ndof': 2900,
                    }
         pa1 = None
@@ -266,7 +266,7 @@ def test_jointcalTask_3_visits_simple_astrometry_no_photometry_min_measurements_
                    'associated_astrometry_fittedStars': 3199,
                    'selected_astrometry_fittedStars': 432,
                    'selected_astrometry_ccdImages': 9,
-                   'astrometry_final_chi2': 360.891,
+                   'astrometry_final_chi2': 360.884,
                    'astrometry_final_ndof': 1294,
                    }
         pa1 = None