Merge pull request #266 from lsst/tickets/DM-13054

DM-13054: Refactor colorterm calculations to use in jointcal

python/lsst/pipe/tasks/colorterms.py

@@ -22,7 +22,9 @@
 import fnmatch
 
 import numpy as np
+import astropy.units as u
 
+from lsst.afw.image import abMagErrFromFluxErr
 import lsst.pex.exceptions as pexExcept
 from lsst.pex.config import Config, Field, ConfigDictField
 from lsst.afw.image import Filter
@@ -58,6 +60,60 @@ class Colorterm(Config):
     c1 = Field(dtype=float, default=0.0, doc="First-order parameter")
     c2 = Field(dtype=float, default=0.0, doc="Second-order parameter")
 
+    def getCorrectedMagnitudes(self, refCat, filterName):
+        """Return the colorterm corrected magnitudes for a given filter.
+
+        Parameters
+        ----------
+        refCat : `lsst.afw.table.SimpleCatalog`
+            The reference catalog to apply color corrections to.
+        filterName : `str`
+            The camera filter to correct the reference catalog into.
+
+        Returns
+        -------
+        RefMag : `np.ndarray`
+            The corrected AB magnitudes.
+        RefMagErr : `np.ndarray`
+            The corrected AB magnitude errors.
+
+        Raises
+        ------
+        KeyError
+            Raised if the reference catalog does not have a flux uncertainty
+            for that filter.
+
+        Notes
+        -----
+        WARNING: I do not know that we can trust the propagation of magnitude
+        errors returned by this method. They need more thorough tests.
+        """
+
+        def getFluxes(fluxField):
+            """Get the flux and fluxErr of this field from refCat."""
+            fluxKey = refCat.schema.find(fluxField).key
+            refFlux = refCat[fluxKey]
+            try:
+                fluxErrKey = refCat.schema.find(fluxField + "Err").key
+                refFluxErr = refCat[fluxErrKey]
+            except KeyError as e:
+                raise KeyError("Reference catalog does not have flux uncertainties for %s" % fluxField) from e
+
+            return refFlux, refFluxErr
+
+        primaryFlux, primaryErr = getFluxes(self.primary + "_flux")
+        secondaryFlux, secondaryErr = getFluxes(self.secondary + "_flux")
+
+        primaryMag = u.Quantity(primaryFlux, u.Jy).to_value(u.ABmag)
+        secondaryMag = u.Quantity(secondaryFlux, u.Jy).to_value(u.ABmag)
+
+        refMag = self.transformMags(primaryMag, secondaryMag)
+        refFluxErrArr = self.propagateFluxErrors(primaryErr, secondaryErr)
+
+        refMagErr = abMagErrFromFluxErr(refFluxErrArr, primaryFlux)
+
+        return refMag, refMagErr
+
     def transformSource(self, source):
         """!Transform the brightness of a source
 

python/lsst/pipe/tasks/photoCal.py

@@ -24,10 +24,12 @@
 import sys
 
 import numpy as np
+import astropy.units as u
 
 import lsst.pex.config as pexConf
 import lsst.pipe.base as pipeBase
 from lsst.afw.image import abMagFromFlux, abMagErrFromFluxErr, Calib
+import lsst.afw.table as afwTable
 from lsst.meas.astrom import DirectMatchTask, DirectMatchConfigWithoutLoader
 import lsst.afw.display.ds9 as ds9
 from lsst.meas.algorithms import getRefFluxField, ReserveSourcesTask
@@ -326,68 +328,50 @@ def extractMagArrays(self, matches, filterName, sourceKeys):
         if applyColorTerms:
             self.log.info("Applying color terms for filterName=%r, config.photoCatName=%s because %s",
                           filterName, self.config.photoCatName, applyCTReason)
-            ct = self.config.colorterms.getColorterm(
+            colorterm = self.config.colorterms.getColorterm(
                 filterName=filterName, photoCatName=self.config.photoCatName, doRaise=True)
-        else:
-            self.log.info("Not applying color terms because %s", applyCTReason)
-            ct = None
+            refCat = afwTable.SimpleCatalog(matches[0].first.schema)
 
-        if ct:                          # we have a color term to worry about
-            fluxFieldList = [getRefFluxField(refSchema, filt) for filt in (ct.primary, ct.secondary)]
-            missingFluxFieldList = []
-            for fluxField in fluxFieldList:
-                try:
-                    refSchema.find(fluxField).key
-                except KeyError:
-                    missingFluxFieldList.append(fluxField)
+            # extract the matched refCat as a Catalog for the colorterm code
+            refCat.reserve(len(matches))
+            for x in matches:
+                record = refCat.addNew()
+                record.assign(x.first)
 
-            if missingFluxFieldList:
-                self.log.warn("Source catalog does not have fluxes for %s; ignoring color terms",
-                              " ".join(missingFluxFieldList))
-                ct = None
+            refMagArr, refMagErrArr = colorterm.getCorrectedMagnitudes(refCat, filterName)
+            fluxFieldList = [getRefFluxField(refSchema, filt) for filt in (colorterm.primary,
+                                                                           colorterm.secondary)]
+        else:
+            # no colorterms to apply
+            self.log.info("Not applying color terms because %s", applyCTReason)
+            colorterm = None
 
-        if not ct:
             fluxFieldList = [getRefFluxField(refSchema, filterName)]
-
-        refFluxArrList = []  # list of ref arrays, one per flux field
-        refFluxErrArrList = []  # list of ref flux arrays, one per flux field
-        for fluxField in fluxFieldList:
+            fluxField = getRefFluxField(refSchema, filterName)
             fluxKey = refSchema.find(fluxField).key
             refFluxArr = np.array([m.first.get(fluxKey) for m in matches])
+
             try:
                 fluxErrKey = refSchema.find(fluxField + "Err").key
                 refFluxErrArr = np.array([m.first.get(fluxErrKey) for m in matches])
             except KeyError:
-                # Reference catalogue may not have flux uncertainties; HACK
+                # Reference catalogue may not have flux uncertainties; HACK DM-2308
                 self.log.warn("Reference catalog does not have flux uncertainties for %s; using sqrt(flux).",
                               fluxField)
                 refFluxErrArr = np.sqrt(refFluxArr)
 
-            refFluxArrList.append(refFluxArr)
-            refFluxErrArrList.append(refFluxErrArr)
-
-        if ct:                          # we have a color term to worry about
-            refMagArr1 = np.array([abMagFromFlux(rf1) for rf1 in refFluxArrList[0]])  # primary
-            refMagArr2 = np.array([abMagFromFlux(rf2) for rf2 in refFluxArrList[1]])  # secondary
-
-            refMagArr = ct.transformMags(refMagArr1, refMagArr2)
-            refFluxErrArr = ct.propagateFluxErrors(refFluxErrArrList[0], refFluxErrArrList[1])
-        else:
-            refMagArr = np.array([abMagFromFlux(rf) for rf in refFluxArrList[0]])
+            refMagArr = u.Quantity(refFluxArr, u.Jy).to_value(u.ABmag)
+            refMagErrArr = abMagErrFromFluxErr(refFluxErrArr, refFluxArr)
 
+        # compute the source catalog magnitudes and errors
         srcMagArr = np.array([abMagFromFlux(sf) for sf in srcInstFluxArr])
-
         # Fitting with error bars in both axes is hard
         # for now ignore reference flux error, but ticket DM-2308 is a request for a better solution
-        magErrArr = np.array([abMagErrFromFluxErr(fe, sf)
-                             for fe, sf in zip(srcInstFluxErrArr, srcInstFluxArr)])
+        magErrArr = abMagErrFromFluxErr(srcInstFluxErrArr, srcInstFluxArr)
         if self.config.magErrFloor != 0.0:
             magErrArr = (magErrArr**2 + self.config.magErrFloor**2)**0.5
 
-        srcMagErrArr = np.array([abMagErrFromFluxErr(sfe, sf)
-                                for sfe, sf in zip(srcInstFluxErrArr, srcInstFluxArr)])
-        refMagErrArr = np.array([abMagErrFromFluxErr(rfe, rf)
-                                for rfe, rf in zip(refFluxErrArr, refFluxArr)])
+        srcMagErrArr = abMagErrFromFluxErr(srcInstFluxErrArr, srcInstFluxArr)
 
         good = np.isfinite(srcMagArr) & np.isfinite(refMagArr)
 

tests/test_colorterm.py

@@ -20,20 +20,13 @@
 # see <http://www.lsstcorp.org/LegalNotices/>.
 #
 
-"""
-Tests for Colorterm
-
-Run with:
-   python colorterm.py
-or
-   python
-   >>> import colorterm; colorterm.run()
-"""
-
 import unittest
 import pickle
 
+import astropy.units as u
+
 import lsst.utils.tests
+from lsst.meas.algorithms import LoadReferenceObjectsTask
 from lsst.pipe.tasks.colorterms import Colorterm, ColortermDict, ColortermLibrary, ColortermNotFoundError
 
 # From the last page of http://www.naoj.org/staff/nakata/suprime/illustration/colorterm_report_ver3.pdf
@@ -117,6 +110,37 @@ def testPickle(self):
         self.assertEqual(colorterms, self.colorterms)
 
 
+def make_fake_refcat(center, flux):
+    """Make a fake reference catalog."""
+    filters = ['f1', 'f2']
+    schema = LoadReferenceObjectsTask.makeMinimalSchema(filters)
+    catalog = lsst.afw.table.SimpleCatalog(schema)
+    record = catalog.addNew()
+    record.setCoord(center)
+    record[filters[0] + '_flux'] = flux
+    record[filters[0] + '_fluxErr'] = flux*0.1
+    record[filters[1] + '_flux'] = flux*10
+    record[filters[1] + '_fluxErr'] = flux*10*0.1
+    return catalog
+
+
+class ApplyColortermsTestCase(lsst.utils.tests.TestCase):
+    def setUp(self):
+        self.colorterm = Colorterm(primary="f1", secondary="f2", c0=2.0, c1=3.0)
+
+    def testGetCorrectedMagnitudes(self):
+        center = lsst.geom.SpherePoint(30, -30, lsst.geom.degrees)
+        flux = 100
+        refCat = make_fake_refcat(center, flux)
+
+        expectMag = self.colorterm.transformMags((u.Jy*refCat['f1_flux']).to_value(u.ABmag),
+                                                 (u.Jy*refCat['f2_flux']).to_value(u.ABmag))
+        refMag, refMagErr = self.colorterm.getCorrectedMagnitudes(refCat, 'r')
+        self.assertEqual(refMag, expectMag)
+        # TODO DM-17692: Not testing the returned errors, as I do not trust `propagateFluxErrors()`
+        # and there is some interesting logic involved in how the errors are propagated.
+
+
 class MemoryTester(lsst.utils.tests.MemoryTestCase):
     pass
 

tests/test_photoCal.py

@@ -113,7 +113,6 @@ def _runTask(self):
         task = PhotoCalTask(self.refObjLoader, config=self.config, schema=self.srcCat.schema)
         pCal = task.run(exposure=self.exposure, sourceCat=self.srcCat)
         matches = pCal.matches
-        print("Ref flux fields list =", pCal.arrays.refFluxFieldList)
         refFluxField = pCal.arrays.refFluxFieldList[0]
 
         # These are *all* the matches; we don't really expect to do that well.