Skip to content

Commit be3b003

Browse files
authored
Merge pull request #232 from lsst/tickets/DM-15751
DM-15751: Allow using jointcal output to make warps
2 parents 83da4df + e720943 commit be3b003

File tree

4 files changed

+210
-35
lines changed

4 files changed

+210
-35
lines changed

python/lsst/pipe/tasks/coaddBase.py

+7-34
Original file line numberDiff line numberDiff line change
@@ -31,11 +31,6 @@
3131
from .coaddInputRecorder import CoaddInputRecorderTask
3232
from .scaleVariance import ScaleVarianceTask
3333

34-
try:
35-
from lsst.meas.mosaic import applyMosaicResults
36-
except ImportError:
37-
applyMosaicResults = None
38-
3934
__all__ = ["CoaddBaseTask", "getSkyInfo"]
4035

4136

@@ -72,7 +67,13 @@ class CoaddBaseConfig(pexConfig.Config):
7267
modelPsf = measAlg.GaussianPsfFactory.makeField(doc="Model Psf factory")
7368
doApplyUberCal = pexConfig.Field(
7469
dtype=bool,
75-
doc="Apply meas_mosaic ubercal results to input calexps?",
70+
doc="Apply jointcal WCS and PhotoCalib results to input calexps?",
71+
default=False
72+
)
73+
useMeasMosaic = pexConfig.Field(
74+
dtype=bool,
75+
doc="Use meas_mosaic's applyMosaicResultsExposure() to do the photometric "
76+
"calibration/wcs update (deprecated).",
7677
default=False
7778
)
7879
matchingKernelSize = pexConfig.Field(
@@ -140,34 +141,6 @@ def getSkyInfo(self, patchRef):
140141
"""
141142
return getSkyInfo(coaddName=self.config.coaddName, patchRef=patchRef)
142143

143-
def getCalExp(self, dataRef, bgSubtracted):
144-
"""!Return one "calexp" calibrated exposure
145-
146-
@param[in] dataRef a sensor-level data reference
147-
@param[in] bgSubtracted return calexp with background subtracted? If False get the
148-
calexp's background background model and add it to the calexp.
149-
@return calibrated exposure
150-
151-
If config.doApplyUberCal, meas_mosaic calibrations will be applied to
152-
the returned exposure using applyMosaicResults.
153-
"""
154-
exposure = dataRef.get("calexp", immediate=True)
155-
if not bgSubtracted:
156-
background = dataRef.get("calexpBackground", immediate=True)
157-
mi = exposure.getMaskedImage()
158-
mi += background.getImage()
159-
del mi
160-
if not self.config.doApplyUberCal:
161-
return exposure
162-
if applyMosaicResults is None:
163-
raise RuntimeError(
164-
"Cannot use improved calibrations for %s because meas_mosaic could not be imported."
165-
% dataRef.dataId
166-
)
167-
else:
168-
applyMosaicResults(dataRef, calexp=exposure)
169-
return exposure
170-
171144
def getCoaddDatasetName(self, warpType="direct"):
172145
"""Return coadd name for given warpType and task config
173146

python/lsst/pipe/tasks/makeCoaddTempExp.py

+67-1
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@
2222
import numpy
2323

2424
import lsst.pex.config as pexConfig
25+
import lsst.daf.persistence as dafPersist
2526
import lsst.afw.image as afwImage
2627
import lsst.coadd.utils as coaddUtils
2728
import lsst.pipe.base as pipeBase
@@ -34,6 +35,14 @@
3435
__all__ = ["MakeCoaddTempExpTask"]
3536

3637

38+
class MissingExposureError(Exception):
39+
"""Raised when data cannot be retrieved for an exposure.
40+
When processing patches, sometimes one exposure is missing; this lets us
41+
distinguish bewteen that case, and other errors.
42+
"""
43+
pass
44+
45+
3746
class MakeCoaddTempExpConfig(CoaddBaseTask.ConfigClass):
3847
"""Config for MakeCoaddTempExpTask
3948
"""
@@ -364,7 +373,7 @@ def run(self, calexpRefList, skyInfo, visitId=0):
364373
# which do.
365374
calExpRef = calExpRef.butlerSubset.butler.dataRef("calexp", dataId=calExpRef.dataId,
366375
tract=skyInfo.tractInfo.getId())
367-
calExp = self.getCalExp(calExpRef, bgSubtracted=self.config.bgSubtracted)
376+
calExp = self.getCalibratedExposure(calExpRef, bgSubtracted=self.config.bgSubtracted)
368377
except Exception as e:
369378
self.log.warn("Calexp %s not found; skipping it: %s", calExpRef.dataId, e)
370379
continue
@@ -430,6 +439,63 @@ def run(self, calexpRefList, skyInfo, visitId=0):
430439
result = pipeBase.Struct(exposures=coaddTempExps)
431440
return result
432441

442+
def getCalibratedExposure(self, dataRef, bgSubtracted):
443+
"""Return one calibrated Exposure, possibly with an updated SkyWcs.
444+
445+
@param[in] dataRef a sensor-level data reference
446+
@param[in] bgSubtracted return calexp with background subtracted? If False get the
447+
calexp's background background model and add it to the calexp.
448+
@return calibrated exposure
449+
450+
@raises MissingExposureError If data for the exposure is not available.
451+
452+
If config.doApplyUberCal, the exposure will be photometrically
453+
calibrated via the `jointcal_photoCalib` dataset and have its SkyWcs
454+
updated to the `jointcal_wcs`, otherwise it will be calibrated via the
455+
Exposure's own Calib and have the original SkyWcs.
456+
"""
457+
try:
458+
exposure = dataRef.get("calexp", immediate=True)
459+
except dafPersist.NoResults as e:
460+
raise MissingExposureError('Exposure not found: %s ' % str(e)) from e
461+
462+
if not bgSubtracted:
463+
background = dataRef.get("calexpBackground", immediate=True)
464+
mi = exposure.getMaskedImage()
465+
mi += background.getImage()
466+
del mi
467+
468+
if self.config.doApplyUberCal:
469+
if self.config.useMeasMosaic:
470+
from lsst.meas.mosaic import applyMosaicResultsExposure
471+
# NOTE: this changes exposure in-place, updating its Calib and Wcs.
472+
# Save the calibration error, as it gets overwritten with zero.
473+
fluxMag0Err = exposure.getCalib().getFluxMag0()[1]
474+
try:
475+
applyMosaicResultsExposure(dataRef, calexp=exposure)
476+
except dafPersist.NoResults as e:
477+
raise MissingExposureError('Mosaic calibration not found: %s ' % str(e)) from e
478+
fluxMag0 = exposure.getCalib().getFluxMag0()[0]
479+
photoCalib = afwImage.PhotoCalib(1.0/fluxMag0,
480+
fluxMag0Err/fluxMag0**2,
481+
exposure.getBBox())
482+
else:
483+
photoCalib = dataRef.get("jointcal_photoCalib")
484+
skyWcs = dataRef.get("jointcal_wcs")
485+
exposure.setWcs(skyWcs)
486+
else:
487+
fluxMag0 = exposure.getCalib().getFluxMag0()
488+
photoCalib = afwImage.PhotoCalib(1.0/fluxMag0[0],
489+
fluxMag0[1]/fluxMag0[0]**2,
490+
exposure.getBBox())
491+
492+
exposure.maskedImage = photoCalib.calibrateImage(exposure.maskedImage)
493+
exposure.maskedImage /= photoCalib.getCalibrationMean()
494+
exposure.setCalib(afwImage.Calib(1/photoCalib.getCalibrationMean()))
495+
# TODO: The images will have a calibration of 1.0 everywhere once RFC-545 is implemented.
496+
# exposure.setCalib(afwImage.Calib(1.0))
497+
return exposure
498+
433499
@staticmethod
434500
def _prepareEmptyExposure(skyInfo):
435501
"""Produce an empty exposure for a given patch"""

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

+3
Original file line numberDiff line numberDiff line change
@@ -189,6 +189,9 @@ def buildInputImages(self, butler, obsCatalog=None, truthCatalog=None, tract=0):
189189
visit = obsRecord.getI(visitKey)
190190
self.log.info("Generating image for visit={visit}, ccd={ccd}".format(ccd=ccd, visit=visit))
191191
exposure = lsst.afw.image.ExposureF(obsRecord.getBBox())
192+
# Apply a tiny offset to the images, so that they have non-zero background.
193+
# If the image background is identically zero, the calculated variance will be NaN.
194+
exposure.maskedImage.image.array += 1e-8
192195
exposure.setCalib(obsRecord.getCalib())
193196
exposure.setWcs(obsRecord.getWcs())
194197
exposure.setPsf(obsRecord.getPsf())

tests/test_makeCoaddTempExp.py

+133
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,133 @@
1+
# This file is part of pipe_tasks.
2+
#
3+
# Developed for the LSST Data Management System.
4+
# This product includes software developed by the LSST Project
5+
# (https://www.lsst.org).
6+
# See the COPYRIGHT file at the top-level directory of this distribution
7+
# for details of code ownership.
8+
#
9+
# This program is free software: you can redistribute it and/or modify
10+
# it under the terms of the GNU General Public License as published by
11+
# the Free Software Foundation, either version 3 of the License, or
12+
# (at your option) any later version.
13+
#
14+
# This program is distributed in the hope that it will be useful,
15+
# but WITHOUT ANY WARRANTY; without even the implied warranty of
16+
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17+
# GNU General Public License for more details.
18+
#
19+
# You should have received a copy of the GNU General Public License
20+
# along with this program. If not, see <https://www.gnu.org/licenses/>.
21+
22+
import unittest
23+
import unittest.mock
24+
25+
import numpy as np
26+
27+
import lsst.utils.tests
28+
29+
import lsst.afw.image
30+
from lsst.daf.persistence import NoResults, ButlerDataRef
31+
from lsst.pipe.tasks.makeCoaddTempExp import (MakeCoaddTempExpTask,
32+
MakeCoaddTempExpConfig,
33+
MissingExposureError)
34+
35+
36+
class GetCalibratedExposureTestCase(lsst.utils.tests.TestCase):
37+
def setUp(self):
38+
np.random.seed(10)
39+
40+
self.config = MakeCoaddTempExpConfig()
41+
42+
# TODO DM-10156: once Calib is replaced, this will be much cleaner
43+
meanCalibration = 1e-4
44+
calibrationErr = 1e-5
45+
self.exposurePhotoCalib = lsst.afw.image.PhotoCalib(meanCalibration, calibrationErr)
46+
# a "jointcal_photoCalib" calibration to return
47+
self.jointcalPhotoCalib = lsst.afw.image.PhotoCalib(1e-6, 1e-8)
48+
49+
crpix = lsst.geom.Point2D(0, 0)
50+
crval = lsst.geom.SpherePoint(0, 45, lsst.geom.degrees)
51+
cdMatrix = lsst.afw.geom.makeCdMatrix(scale=1.0*lsst.geom.arcseconds)
52+
self.skyWcs = lsst.afw.geom.makeSkyWcs(crpix, crval, cdMatrix)
53+
jointcalCdMatrix = lsst.afw.geom.makeCdMatrix(scale=0.9*lsst.geom.arcseconds)
54+
# a "jointcal_wcs" skyWcs to return
55+
self.jointcalSkyWcs = lsst.afw.geom.makeSkyWcs(crpix, crval, jointcalCdMatrix)
56+
57+
self.exposure = lsst.afw.image.ExposureF(10, 10)
58+
self.exposure.maskedImage.image.array = np.random.random((10, 10)).astype(np.float32) * 1000
59+
self.exposure.maskedImage.variance.array = np.random.random((10, 10)).astype(np.float32)
60+
# mask at least one pixel
61+
self.exposure.maskedImage.mask[5, 5] = 3
62+
# set the Calib and Wcs objects of this exposure.
63+
self.exposure.getCalib().setFluxMag0(1/meanCalibration, calibrationErr/meanCalibration**2)
64+
self.exposure.setWcs(self.skyWcs)
65+
66+
# set to True in a test to raise NoResults for get('calexp')
67+
self.raiseOnGetCalexp = False
68+
69+
def mockGet(datasetType, dataId=None, immediate=None):
70+
"""Minimally fake a butler.get()."""
71+
if "calexp" in datasetType:
72+
if self.raiseOnGetCalexp:
73+
raise NoResults("failed!", datasetType, dataId)
74+
else:
75+
return self.exposure.clone()
76+
if "jointcal_photoCalib" in datasetType:
77+
return self.jointcalPhotoCalib
78+
if "jointcal_wcs" in datasetType:
79+
return self.jointcalSkyWcs
80+
81+
self.dataRef = unittest.mock.Mock(spec=ButlerDataRef)
82+
self.dataRef.get.side_effect = mockGet
83+
self.dataRef.dataId = {"ccd": 10000, "visit": 1}
84+
85+
def test_getCalibratedExposureGetCalexpRaises(self):
86+
"""If get('calexp') raises NoResults, we should get a usefully
87+
chained exception.
88+
"""
89+
task = MakeCoaddTempExpTask(config=self.config)
90+
91+
self.raiseOnGetCalexp = True
92+
93+
with self.assertRaises(MissingExposureError) as cm:
94+
task.getCalibratedExposure(self.dataRef, True)
95+
self.assertIsInstance(cm.exception.__cause__, NoResults)
96+
self.assertIn('Exposure not found', str(cm.exception))
97+
98+
def test_getCalibratedExposure(self):
99+
task = MakeCoaddTempExpTask(config=self.config)
100+
101+
expect = self.exposurePhotoCalib.calibrateImage(self.exposure.maskedImage)
102+
expect /= self.exposurePhotoCalib.getCalibrationMean()
103+
result = task.getCalibratedExposure(self.dataRef, True)
104+
105+
self.assertMaskedImagesEqual(result.maskedImage, expect)
106+
# TODO: once RFC-545 is implemented, this should be 1.0
107+
self.assertEqual(result.getCalib().getFluxMag0()[0], 1/self.exposurePhotoCalib.getCalibrationMean())
108+
self.assertEqual(result.getWcs(), self.skyWcs)
109+
110+
def test_getCalibratedExposureJointcal(self):
111+
self.config.doApplyUberCal = True
112+
task = MakeCoaddTempExpTask(config=self.config)
113+
114+
expect = self.jointcalPhotoCalib.calibrateImage(self.exposure.maskedImage)
115+
expect /= self.jointcalPhotoCalib.getCalibrationMean()
116+
result = task.getCalibratedExposure(self.dataRef, True)
117+
self.assertMaskedImagesEqual(result.maskedImage, expect)
118+
# TODO: once RFC-545 is implemented, this should be 1.0
119+
self.assertEqual(result.getCalib().getFluxMag0()[0], 1/self.jointcalPhotoCalib.getCalibrationMean())
120+
self.assertEqual(result.getWcs(), self.jointcalSkyWcs)
121+
122+
123+
def setup_module(module):
124+
lsst.utils.tests.init()
125+
126+
127+
class MatchMemoryTestCase(lsst.utils.tests.MemoryTestCase):
128+
pass
129+
130+
131+
if __name__ == "__main__":
132+
lsst.utils.tests.init()
133+
unittest.main()

0 commit comments

Comments
 (0)