Merge pull request #196 from lsst/tickets/DM-24472

DM-24472: Regenerate Gaia DR2 catalogs to correct coordinate error fields

doc/lsst.meas.algorithms/creating-a-reference-catalog.rst

@@ -69,6 +69,7 @@ This is an example configuration that was used to ingest the Gaia DR2 catalog:
 
     config.parallax_name = "parallax"
     config.parallax_err_name = "parallax_error"
+    config.coord_err_unit = "milliarcsecond"
 
     config.pm_ra_name = "pmra"
     config.pm_ra_err_name = "pmra_error"

python/lsst/meas/algorithms/ingestIndexManager.py

@@ -80,6 +80,8 @@ def __init__(self, filenames, config, file_reader, indexer,
         self.htmRange = htmRange
         self.addRefCatMetadata = addRefCatMetadata
         self.log = log
+        # cache this to speed up coordinate conversions
+        self.coord_err_unit = u.Unit(self.config.coord_err_unit)
 
     def run(self, inputFiles):
         """Index a set of input files from a reference catalog, and write the
@@ -120,12 +122,13 @@ def _ingestOneFile(self, filename, fileLocks):
         global FILE_PROGRESS
         inputData = self.file_reader.run(filename)
         fluxes = self._getFluxes(inputData)
+        coordErr = self._getCoordErr(inputData)
         matchedPixels = self.indexer.indexPoints(inputData[self.config.ra_name],
                                                  inputData[self.config.dec_name])
         pixel_ids = set(matchedPixels)
         for pixelId in pixel_ids:
             with fileLocks[pixelId]:
-                self._doOnePixel(inputData, matchedPixels, pixelId, fluxes)
+                self._doOnePixel(inputData, matchedPixels, pixelId, fluxes, coordErr)
         with FILE_PROGRESS.get_lock():
             oldPercent = 100 * FILE_PROGRESS.value / self.nInputFiles
             FILE_PROGRESS.value += 1
@@ -137,7 +140,7 @@ def _ingestOneFile(self, filename, fileLocks):
                               self.nInputFiles,
                               percent)
 
-    def _doOnePixel(self, inputData, matchedPixels, pixelId, fluxes):
+    def _doOnePixel(self, inputData, matchedPixels, pixelId, fluxes, coordErr):
         """Process one HTM pixel, appending to an existing catalog or creating
         a new catalog, as needed.
 
@@ -152,6 +155,9 @@ def _doOnePixel(self, inputData, matchedPixels, pixelId, fluxes):
         fluxes : `dict` [`str`, `numpy.ndarray`]
             The values that will go into the flux and fluxErr fields in the
             output catalog.
+        coordErr : `dict` [`str`, `numpy.ndarray`]
+            The values that will go into the coord_raErr, coord_decErr, and
+            coord_ra_dec_Cov fields in the output catalog (in radians).
         """
         idx = np.where(matchedPixels == pixelId)[0]
         catalog = self.getCatalog(pixelId, self.schema, len(idx))
@@ -162,9 +168,14 @@ def _doOnePixel(self, inputData, matchedPixels, pixelId, fluxes):
         with COUNTER.get_lock():
             self._setIds(inputData[idx], catalog)
 
+        # set fluxes from the pre-computed array
         for name, array in fluxes.items():
             catalog[self.key_map[name]][-len(idx):] = array[idx]
 
+        # set coordinate errors from the pre-computed array
+        for name, array in coordErr.items():
+            catalog[name][-len(idx):] = array[idx]
+
         catalog.writeFits(self.filenames[pixelId])
 
     def _setIds(self, inputData, catalog):
@@ -238,22 +249,32 @@ def computeCoord(row, ra_name, dec_name):
         """
         return lsst.geom.SpherePoint(row[ra_name], row[dec_name], lsst.geom.degrees)
 
-    def _setCoordErr(self, record, row):
-        """Set coordinate error in a record of an indexed catalog.
-
-        The errors are read from the specified columns, and installed
-        in the appropriate columns of the output.
+    def _getCoordErr(self, inputData, ):
+        """Compute the ra/dec error fields that will go into the output catalog.
 
         Parameters
         ----------
-        record : `lsst.afw.table.SimpleRecord`
-            Row from indexed catalog to modify.
-        row : `numpy.ndarray`
-            Row from catalog being ingested.
+        inputData : `numpy.ndarray`
+            The input data to compute fluxes for.
+
+        Returns
+        -------
+        coordErr : `dict` [`str`, `numpy.ndarray`]
+            The values that will go into the coord_raErr, coord_decErr, fields
+            in the output catalog (in radians).
+
+        Notes
+        -----
+        This does not currently handle the ra/dec covariance field,
+        ``coord_ra_dec_Cov``. That field may require extra work, as its units
+        may be more complicated in external catalogs.
         """
-        if self.config.ra_err_name:  # IngestIndexedReferenceConfig.validate ensures all or none
-            record.set(self.key_map["coord_raErr"], np.radians(row[self.config.ra_err_name]))
-            record.set(self.key_map["coord_decErr"], np.radians(row[self.config.dec_err_name]))
+        result = {}
+        result['coord_raErr'] = u.Quantity(inputData[self.config.ra_err_name],
+                                           self.coord_err_unit).to_value(u.radian)
+        result['coord_decErr'] = u.Quantity(inputData[self.config.dec_err_name],
+                                            self.coord_err_unit).to_value(u.radian)
+        return result
 
     def _setFlags(self, record, row):
         """Set flags in an output record.
@@ -372,7 +393,6 @@ def _fillRecord(self, record, row):
         """
         record.setCoord(self.computeCoord(row, self.config.ra_name, self.config.dec_name))
 
-        self._setCoordErr(record, row)
         self._setFlags(record, row)
         self._setProperMotion(record, row)
         self._setParallax(record, row)

python/lsst/meas/algorithms/ingestIndexReferenceTask.py

@@ -26,6 +26,8 @@
 
 import os.path
 
+import astropy.units
+
 import lsst.pex.config as pexConfig
 import lsst.pipe.base as pipeBase
 import lsst.geom
@@ -126,11 +128,11 @@ class IngestIndexedReferenceConfig(pexConfig.Config):
     )
     ra_name = pexConfig.Field(
         dtype=str,
-        doc="Name of RA column",
+        doc="Name of RA column (values in decimal degrees)",
     )
     dec_name = pexConfig.Field(
         dtype=str,
-        doc="Name of Dec column",
+        doc="Name of Dec column (values in decimal degrees)",
     )
     ra_err_name = pexConfig.Field(
         dtype=str,
@@ -142,6 +144,11 @@ class IngestIndexedReferenceConfig(pexConfig.Config):
         doc="Name of Dec error column",
         optional=True,
     )
+    coord_err_unit = pexConfig.Field(
+        dtype=str,
+        doc="Unit of RA/Dec error fields (astropy.unit.Unit compatible)",
+        optional=True
+    )
     mag_column_list = pexConfig.ListField(
         dtype=str,
         doc="The values in the reference catalog are assumed to be in AB magnitudes. "
@@ -259,7 +266,13 @@ def assertAllOrNone(*names):
             raise ValueError(
                 "mag_err_column_map specified, but keys do not match mag_column_list: {} != {}".format(
                     sorted(self.mag_err_column_map.keys()), sorted(self.mag_column_list)))
-        assertAllOrNone("ra_err_name", "dec_err_name")
+        assertAllOrNone("ra_err_name", "dec_err_name", "coord_err_unit")
+        if self.coord_err_unit is not None:
+            result = astropy.units.Unit(self.coord_err_unit, parse_strict='silent')
+            if isinstance(result, astropy.units.UnrecognizedUnit):
+                msg = f"{self.coord_err_unit} is not a valid astropy unit string."
+                raise pexConfig.FieldValidationError(IngestIndexedReferenceConfig.coord_err_unit, self, msg)
+
         assertAllOrNone("epoch_name", "epoch_format", "epoch_scale")
         assertAllOrNone("pm_ra_name", "pm_dec_name")
         assertAllOrNone("pm_ra_err_name", "pm_dec_err_name")
@@ -281,12 +294,12 @@ class IngestIndexedReferenceTask(pipeBase.CmdLineTask):
 
     For producing catalogs this task makes the following assumptions
     about the input catalogs:
-    - RA, Dec, RA error and Dec error are all in decimal degrees.
+    - RA, Dec are in decimal degrees.
     - Epoch is available in a column, in a format supported by astropy.time.Time.
     - There are no off-diagonal covariance terms, such as covariance
-        between RA and Dec, or between PM RA and PM Dec. Gaia is a well
-        known example of a catalog that has such terms, and thus should not
-        be ingested with this task.
+      between RA and Dec, or between PM RA and PM Dec. Support for such
+     covariance would have to be added to to the config, including consideration
+     of the units in the input catalog.
 
     Parameters
     ----------