Merge pull request #183 from lsst/tickets/DM-12374

DM-12374: propagate transmission curves in coaddition

python/lsst/pipe/tasks/assembleCoadd.py

@@ -139,6 +139,11 @@ class AssembleCoaddConfig(CoaddBaseTask.ConfigClass):
         doc="Configuration for CoaddPsf",
         dtype=measAlg.CoaddPsfConfig,
     )
+    doAttachTransmissionCurve = pexConfig.Field(
+        dtype=bool, default=False, optional=False,
+        doc=("Attach a piecewise TransmissionCurve for the coadd? "
+             "(requires all input Exposures to have TransmissionCurves).")
+    )
 
     def setDefaults(self):
         CoaddBaseTask.ConfigClass.setDefaults(self)
@@ -574,6 +579,9 @@ def assembleMetadata(self, coaddExposure, tempExpRefList, weightList):
         apCorrMap = measAlg.makeCoaddApCorrMap(coaddInputs.ccds, coaddExposure.getBBox(afwImage.PARENT),
                                                coaddExposure.getWcs())
         coaddExposure.getInfo().setApCorrMap(apCorrMap)
+        if self.config.doAttachTransmissionCurve:
+            transmissionCurve = measAlg.makeCoaddTransmissionCurve(coaddExposure.getWcs(), coaddInputs.ccds)
+            coaddExposure.getInfo().setTransmissionCurve(transmissionCurve)
 
     def assembleSubregion(self, coaddExposure, bbox, tempExpRefList, imageScalerList, weightList,
                           altMaskList, statsFlags, statsCtrl, nImage=None):

python/lsst/pipe/tasks/coaddInputRecorder.py

@@ -39,6 +39,7 @@ class CoaddInputRecorderConfig(pexConfig.Config):
     or the schemas created by earlier tasks (like MakeCoaddTempExpTask) will not contain
     the fields filled by later tasks (like AssembleCoaddTask).
     """
+
     saveEmptyCcds = pexConfig.Field(
         dtype=bool, default=False, optional=False,
         doc=("Add records for CCDs we iterated over but did not add a coaddTempExp"
@@ -62,8 +63,8 @@ class CoaddInputRecorderConfig(pexConfig.Config):
 
 
 class CoaddTempExpInputRecorder(object):
-    """A helper class for CoaddInputRecorderTask, managing the CoaddInputs object for that a single
-    CoaddTempExp.  This will contain single 'visit' record for the CoaddTempExp and a number of 'ccd'
+    """A helper class for CoaddInputRecorderTask, managing the CoaddInputs object for that single
+    CoaddTempExp.  This will contain a single 'visit' record for the CoaddTempExp and a number of 'ccd'
     records.
 
     Should generally be created by calling CoaddInputRecorderTask.makeCoaddTempExp().
@@ -149,6 +150,8 @@ def _setExposureInfoInRecord(self, exposure, record):
         record.setValidPolygon(info.getValidPolygon())
         record.setVisitInfo(info.getVisitInfo())
         record.setBBox(exposure.getBBox())
+        record.setTransmissionCurve(info.getTransmissionCurve())
+
 
 class CoaddInputRecorderTask(pipeBase.Task):
     """Subtask that handles filling a CoaddInputs object for a coadd exposure, tracking the CCDs and

python/lsst/pipe/tasks/mocks/mockCoadd.py

@@ -195,6 +195,7 @@ def buildInputImages(self, butler, obsCatalog=None, truthCatalog=None, tract=0):
             exposure.setWcs(obsRecord.getWcs())
             exposure.setPsf(obsRecord.getPsf())
             exposure.getInfo().setApCorrMap(obsRecord.getApCorrMap())
+            exposure.getInfo().setTransmissionCurve(obsRecord.getTransmissionCurve())
             for truthRecord in truthCatalog:
                 status = self.mockObject.drawSource(truthRecord, exposure, buffer=self.config.edgeBuffer)
                 if status:
@@ -251,6 +252,7 @@ def makeCoaddTask(self, cls, assemblePsfMatched=False):
                 config.doWrite = False
             if cls == CompareWarpAssembleCoaddTask:
                 config.assembleStaticSkyModel.select.retarget(MockSelectImagesTask)
+            config.doAttachTransmissionCurve = True
         return cls(config=config)
 
     def iterPatchRefs(self, butler, tractInfo):

python/lsst/pipe/tasks/mocks/mockObservation.py

@@ -33,6 +33,7 @@
 import lsst.afw.detection
 import lsst.pipe.base
 from lsst.meas.base.apCorrRegistry import getApCorrNameSet
+from lsst.meas.algorithms.testUtils import makeRandomTransmissionCurve
 
 
 class MockObservationConfig(lsst.pex.config.Config):
@@ -131,6 +132,7 @@ def run(self, butler, n, tractInfo, camera, catalog=None):
                 record.setVisitInfo(visitInfo)
                 record.setPsf(self.buildPsf(detector))
                 record.setApCorrMap(self.buildApCorrMap(detector))
+                record.setTransmissionCurve(self.buildTransmissionCurve(detector))
                 record.setBBox(detector.getBBox())
                 detectorId = detector.getId()
                 obj = butler.get("ccdExposureId", visit=visit, ccd=detectorId, immediate=True)
@@ -247,3 +249,8 @@ def makeRandomBoundedField():
             apCorrMap.set(name + "_flux", makeRandomBoundedField())
             apCorrMap.set(name + "_fluxSigma", makeRandomBoundedField())
         return apCorrMap
+
+    def buildTransmissionCurve(self, detector):
+        """Build a random spacially-varying TransmissionCurve."""
+        bbox = detector.getBBox()
+        return makeRandomTransmissionCurve(rng=self.rng, maxRadius=max(bbox.getWidth(), bbox.getHeight()))

tests/nopytest_test_coadds.py

@@ -56,6 +56,7 @@
 import os
 import numbers
 from collections import Iterable
+import numpy as np
 
 import lsst.utils.tests
 import lsst.afw.math
@@ -286,6 +287,7 @@ def testTempExpInputs(self, tract=0):
                         self.assertEqual(obsRecord.getId(), ccdRecord.getId())
                         self.assertEqual(obsRecord.getWcs(), ccdRecord.getWcs())
                         self.assertEqual(obsRecord.getBBox(), ccdRecord.getBBox())
+                        self.assertIsNotNone(ccdRecord.getTransmissionCurve())
                         self.comparePsfs(obsRecord.getPsf(), ccdRecord.getPsf())
                 self.assertTrue(foundOneTempExp)
 
@@ -311,6 +313,7 @@ def testCoaddInputs(self, tract=0):
                     self.assertEqual(obsRecord.getBBox(), ccdRecord.getBBox())
                     self.assertEqual(obsRecord.get("filter"), ccdRecord.get("filter"))
                     self.comparePsfs(obsRecord.getPsf(), ccdRecord.getPsf())
+                    self.assertIsNotNone(ccdRecord.getTransmissionCurve())
                     self.assertIsNotNone(coaddInputs.visits.find(ccdRecord.getL(ccdVisitKey)))
                 for visitRecord in coaddInputs.visits:
                     nCcds = len([ccdRecord for ccdRecord in coaddInputs.ccds
@@ -399,6 +402,41 @@ def testCoaddPsf(self, tract=0):
             print("WARNING: CoaddPsf test inconclusive (this can occur randomly, but very rarely; "
                   "first try running the test again)")
 
+    def testCoaddTransmissionCurves(self, tract=0):
+        """Test that coadded TransmissionCurves agree with those of the inputs."""
+        skyMap = self.butler.get(self.mocksTask.config.coaddName + "Coadd_skyMap", immediate=True)
+        tractInfo = skyMap[tract]
+        truthCatalog = self.butler.get("truth", tract=tract, immediate=True)
+        wavelengths = np.linspace(4000, 7000, 10)
+        for dataProduct in self.coaddNameList:
+            nTested = 0
+            for patchRef in self.mocksTask.iterPatchRefs(self.butler, tractInfo):
+                coaddExp = patchRef.get(self.mocksTask.config.coaddName + dataProduct, immediate=True)
+                coaddWcs = coaddExp.getWcs()
+                coaddTransmissionCurve = coaddExp.getInfo().getTransmissionCurve()
+                coaddBBox = lsst.afw.geom.Box2D(coaddExp.getBBox())
+                inputs = coaddExp.getInfo().getCoaddInputs().ccds
+                for truthRecord in truthCatalog:
+                    coaddPosition = coaddWcs.skyToPixel(truthRecord.getCoord())
+                    if not coaddBBox.contains(coaddPosition):
+                        continue
+                    summedThroughput = np.zeros(wavelengths.shape, dtype=float)
+                    weightSum = 0.0
+                    for sensorRecord in inputs.subsetContaining(truthRecord.getCoord(),
+                                                                includeValidPolygon=True):
+                        sensorPosition = sensorRecord.getWcs().skyToPixel(truthRecord.getCoord())
+                        sensorTransmission = sensorRecord.getTransmissionCurve()
+                        weight = sensorRecord.get("weight")
+                        summedThroughput += sensorTransmission.sampleAt(sensorPosition, wavelengths)*weight
+                        weightSum += weight
+                    if weightSum == 0.0:
+                        continue
+                    summedThroughput /= weightSum
+                    coaddThroughput = coaddTransmissionCurve.sampleAt(coaddPosition, wavelengths)
+                    self.assertFloatsAlmostEqual(coaddThroughput, summedThroughput, rtol=1E-10)
+                    nTested += 1
+            self.assertGreater(nTested, 5)
+
     def testSchemaConsistency(self):
         """Test that _schema catalogs are consistent with the data catalogs.
         """