Add astrometryReferenceErr config and refcat coord errs

This allows us to change the reference catalog coordinate errors if they are
not specified in the refcat itself.

NOTE: Unfortunately, none of the current testdata_jointcal refcats have
coordinate errors, so I can't properly test the code that uses the refcat
errors.

include/lsst/jointcal/Associations.h

@@ -135,11 +135,15 @@ class Associations {
      * @param[in]  matchCut       Separation radius to match fitted and
      *                            reference stars.
      * @param      fluxField      The field name in refCat to get the flux from.
+     * @param      refCoordinateErr Error on reference catalog coordinates (mas). If not NaN, this
+     *                              overrides the `coord_*_err` values in the reference catalog itself.
+     *                              This value is divided by cos(dec) before being used for ra_err.
      * @param      rejectBadFluxes  Reject reference sources with flux=NaN or 0 and/or fluxErr=NaN or 0.
      *                              Typically false for astrometry and true for photometry.
      */
     void collectRefStars(afw::table::SimpleCatalog &refCat, afw::geom::Angle matchCut,
-                         std::string const &fluxField, bool rejectBadFluxes = false);
+                         std::string const &fluxField, float const refCoordinateErr,
+                         bool rejectBadFluxes = false);
 
     //! Sends back the fitted stars coordinates on the sky FittedStarsList::inTangentPlaneCoordinates keeps
     //! track of that.

python/lsst/jointcal/associations.cc

@@ -46,7 +46,7 @@ void declareAssociations(py::module &mod) {
     cls.def("associateCatalogs", &Associations::associateCatalogs, "matchCutInArcsec"_a = 0,
             "useFittedList"_a = false, "enlargeFittedList"_a = true);
     cls.def("collectRefStars", &Associations::collectRefStars, "refCat"_a, "matchCut"_a, "fluxField"_a,
-            "rejectBadFluxes"_a = false);
+            "refCoordinateErr"_a, "rejectBadFluxes"_a = false);
     cls.def("deprojectFittedStars", &Associations::deprojectFittedStars);
     cls.def("nCcdImagesValidForFit", &Associations::nCcdImagesValidForFit);
     cls.def("nFittedStarsWithAssociatedRefStar", &Associations::nFittedStarsWithAssociatedRefStar);

python/lsst/jointcal/jointcal.py

@@ -290,6 +290,13 @@ class JointcalConfig(pexConfig.Config):
         target=ReferenceSourceSelectorTask,
         doc="How to down-select the loaded photometry reference catalog.",
     )
+    astrometryReferenceErr = pexConfig.Field(
+        doc="Uncertainty on reference catalog coordinates [mas] if they do not have `coord_*_err` fields."
+            " If None, then raise an exception if the reference catalog is missing coordinate errors.",
+        dtype=float,
+        default=None,
+        optional=True
+    )
     writeInitMatrix = pexConfig.Field(
         dtype=bool,
         doc="Write the pre/post-initialization Hessian and gradient to text files, for debugging."
@@ -606,8 +613,16 @@ def _do_load_refcat_and_fit(self, associations, defaultFilter, center, radius,
                                                          center, radius, defaultFilter,
                                                          applyColorterms=applyColorterms)
 
-        associations.collectRefStars(refCat, self.config.matchCut*afwGeom.arcseconds,
-                                     fluxField, reject_bad_fluxes)
+        if self.config.astrometryReferenceErr is None:
+            refCoordErr = float('nan')
+        else:
+            refCoordErr = self.config.astrometryReferenceErr
+
+        associations.collectRefStars(refCat,
+                                     self.config.matchCut*afwGeom.arcseconds,
+                                     fluxField,
+                                     refCoordinateErr=refCoordErr,
+                                     rejectBadFluxes=reject_bad_fluxes)
         add_measurement(self.job, 'jointcal.collected_%s_refStars' % name,
                         associations.refStarListSize())
 
@@ -671,6 +686,17 @@ def _load_reference_catalog(self, refObjLoader, referenceSelector, center, radiu
         else:
             refCat = selected.sourceCat
 
+        if self.config.astrometryReferenceErr is None and 'coord_ra_err' not in refCat.schema:
+            #####################
+            # TODO: is KeyError really the best thing to raise here?
+            # RuntimeError also feels too generic, but may be better?
+            #####################
+            raise KeyError("Reference catalog does not contain coordinate errors, and"
+                           "config.astrometryReferenceErr not supplied.")
+        if self.config.astrometryReferenceErr is not None and 'coord_ra_err' in refCat.schema:
+            self.log.warn("Overriding reference catalog coordinate errors with %f/coordinate [mas]",
+                          self.config.astrometryReferenceErr)
+
         if applyColorterms:
             try:
                 refCatName = refObjLoader.ref_dataset_name

src/Associations.cc

@@ -23,6 +23,7 @@
  */
 
 #include <cmath>
+#include <math.h>  // for isnan
 #include <iostream>
 #include <limits>
 #include <sstream>
@@ -170,15 +171,24 @@ void Associations::associateCatalogs(const double matchCutInArcSec, const bool u
 }
 
 void Associations::collectRefStars(afw::table::SimpleCatalog &refCat, afw::geom::Angle matchCut,
-                                   std::string const &fluxField, bool rejectBadFluxes) {
+                                   std::string const &fluxField, float const refCoordinateErr,
+                                   bool rejectBadFluxes) {
     if (refCat.size() == 0) {
         throw(LSST_EXCEPT(pex::exceptions::InvalidParameterError,
                           " reference catalog is empty : stop here "));
     }
 
     afw::table::CoordKey coordKey = refCat.getSchema()["coord"];
+    // Handle reference catalogs that don't have position errors.
+    afw::table::Key<double> raErrKey;
+    afw::table::Key<double> decErrKey;
+    if (isnan(refCoordinateErr)) {
+        raErrKey = refCat.getSchema().find<double>("coord_ra_err").key;
+        decErrKey = refCat.getSchema().find<double>("coord_dec_err").key;
+    }
+
     auto fluxKey = refCat.getSchema().find<double>(fluxField).key;
-    // Don't blow up if the reference catalog doesn't contain errors.
+    // Handle reference catalogs that don't have flux errors.
     afw::table::Key<double> fluxErrKey;
     try {
         fluxErrKey = refCat.getSchema().find<double>(fluxField + "Err").key;
@@ -204,12 +214,15 @@ void Associations::collectRefStars(afw::table::SimpleCatalog &refCat, afw::geom:
         double dec = lsst::afw::geom::radToDeg(coord.getLatitude());
         auto star = std::make_shared<RefStar>(ra, dec, flux, fluxErr);
 
-        // TODO DM-10826: RefCats aren't guaranteed to have position errors.
-        // TODO: Need to devise a way to check whether the refCat has position errors
-        // TODO: and use them instead, if available.
-        // cook up errors: 100 mas per cooordinate
-        star->vx = std::pow(0.1 / 3600 / cos(coord.getLatitude()), 2);
-        star->vy = std::pow(0.1 / 3600, 2);
+        if (isnan(refCoordinateErr)) {
+            star->vx = record->get(raErrKey);
+            star->vy = record->get(decErrKey);
+        } else {
+            // Add fake errors when the catalog has none
+            star->vx = std::pow(refCoordinateErr / 3600. / cos(coord.getLatitude()), 2);
+            star->vy = std::pow(refCoordinateErr / 3600, 2);
+        }
+        // TODO: cook up a covariance as none of our current refcats have it
         star->vxy = 0.;
 
         // Reject sources with non-finite fluxes and flux errors, and fluxErr=0 (which gives chi2=inf).

tests/config/astrometryReferenceErr-None-config.py

@@ -0,0 +1 @@
+config.astrometryReferenceErr = None

tests/config/astrometryReferenceErr-config.py

@@ -0,0 +1 @@
+config.astrometryReferenceErr = 0.01

tests/config/config.py

@@ -1,2 +1,3 @@
 config.astrometryRefObjLoader.ref_dataset_name = "sdss-dr9-fink-v5b"
 config.photometryRefObjLoader.ref_dataset_name = "sdss-dr9-fink-v5b"
+config.astrometryReferenceErr = 0.1  # milli-arcsecond

tests/test_ccdImage.py

@@ -86,7 +86,7 @@ def checkCountStars(self, ccdImage, nStars):
         skyWcs = ccdImage.getReadWcs().getSkyWcs()
         self.refCat = testUtils.createFakeCatalog(nStars, self.bbox, "refFlux", skyWcs=skyWcs, refCat=True)
         # associate the reference stars
-        self.associations.collectRefStars(self.refCat, matchCut, 'refFlux_instFlux')
+        self.associations.collectRefStars(self.refCat, matchCut, 'refFlux_instFlux', 0.1)
         measStars, refStars = ccdImage.countStars()
         self.assertEqual(measStars, nStars)
         self.assertEqual(refStars, nStars)

tests/test_jointcal.py

@@ -190,6 +190,7 @@ def test_load_reference_catalog(self):
         refObjLoader, center, radius, filterName, fakeRefCat = self._make_fake_refcat()
 
         config = lsst.jointcal.jointcal.JointcalConfig()
+        config.astrometryReferenceErr = 0.1  # our test refcats don't have coord errors
         jointcal = lsst.jointcal.JointcalTask(config=config, butler=self.butler)
 
         refCat, fluxField = jointcal._load_reference_catalog(refObjLoader,
@@ -210,6 +211,7 @@ def test_load_reference_catalog_subselect(self):
         refObjLoader, center, radius, filterName, fakeRefCat = self._make_fake_refcat()
 
         config = lsst.jointcal.jointcal.JointcalConfig()
+        config.astrometryReferenceErr = 0.1  # our test refcats don't have coord errors
         config.astrometryReferenceSelector.doSignalToNoise = True
         config.astrometryReferenceSelector.signalToNoise.minimum = 1e10
         config.astrometryReferenceSelector.signalToNoise.fluxField = "fake_flux"

tests/test_jointcal_cfht.py

@@ -154,6 +154,29 @@ def test_jointcalTask_2_visits_constrainedAstrometry_4sigma_outliers(self):
 
         self._testJointcalTask(2, dist_rms_relative, self.dist_rms_absolute, None, metrics=metrics)
 
+    def test_jointcalTask_2_visits_constrainedAstrometry_astrometryReferenceUncertainty_smaller(self):
+        """Test with a smaller fake reference uncertainty: chi2 will be higher."""
+        dist_rms_relative, metrics = self.setup_jointcalTask_2_visits_constrainedAstrometry()
+        test_config = os.path.join(lsst.utils.getPackageDir('jointcal'),
+                                   'tests/config/astrometryReferenceErr-config.py')
+        self.configfiles.append(test_config)
+        metrics['astrometry_final_chi2'] = 11522.9
+        metrics['astrometry_final_ndof'] = 3406
+
+        self._testJointcalTask(2, dist_rms_relative, self.dist_rms_absolute, None, metrics=metrics)
+
+    def test_jointcalTask_2_visits_constrainedAstrometry_astrometryReferenceUncertainty_None_fails(self):
+        """The default `None` should fail for the existing refcats that have no position errors."""
+        dist_rms_relative, metrics = self.setup_jointcalTask_2_visits_constrainedAstrometry()
+        # This is the default, but we override it in tests/config/config.py,
+        # because none of the existing test refcats have errors. So we have to
+        # re-override it here.
+        test_config = os.path.join(lsst.utils.getPackageDir('jointcal'),
+                                   'tests/config/astrometryReferenceErr-None-config.py')
+        self.configfiles.append(test_config)
+        with self.assertRaises(KeyError):
+            self._testJointcalTask(2, None, None, None)
+
     def setup_jointcalTask_2_visits_constrainedPhotometry(self):
         """Set default values for the constrainedPhotometry tests, and make
         the differences between each test and the defaults more obvious.