From 0994555f1e00b7095e63ae2bf93e29011dd3d395 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Fri, 16 Oct 2026 12:30:52 -0700 Subject: [PATCH 1/3] Record the template instrument in GetTemplateTask A template is built from coadds of a single instrument: the science image's, if any coadd was made from it, and otherwise the one most coadds were made from. Coadds from any other instrument are dropped. --- python/lsst/ip/diffim/getTemplate.py | 150 ++++++++++++++---- tests/test_getTemplate.py | 219 +++++++++++++++++++++++++-- 2 files changed, 330 insertions(+), 39 deletions(-) diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index 92d18dd7..76290066 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -59,6 +59,11 @@ _PSF_STAMP_SIGMAS = 5 """Width, in N sigma, of the Gaussian approximation to the CoaddPsf.""" +_UNKNOWN_INSTRUMENT = "UNKNOWN" +"""Instrument recorded for a coadd or template made from observations of an +unknown instrument. +""" + class GetTemplateConnections( pipeBase.PipelineTaskConnections, @@ -267,6 +272,12 @@ def __init__(self, *args, **kwargs): type=float, doc="Weight for each exposure, used to make the CoaddPsf; should always be 1.", ) + self.schema.addField( + "instrument", + type=str, + size=64, + doc="Instrument whose observations were coadded to make the template, or UNKNOWN.", + ) self.makeSubtask("varianceBackground") def runQuantum(self, butlerQC, inputRefs, outputRefs): @@ -290,6 +301,7 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): dataIds=results.dataIds, physical_filter=physical_filter, visit=outputRefs.template.dataId["visit"], + instrument=butlerQC.quantum.dataId["instrument"], ) if self.config.output_image_type == "future": dataId = butlerQC.quantum.dataId @@ -305,9 +317,9 @@ def convert_outputs_to_future(self, template, coadd_refs, detector, exposureReco """Convert a template image to DifferenceImage format. This replaces ``template`` with an `lsst.images.DifferenceImage` - instance and records the coadds contributed to it. The Gaussian PSF - set by `run` is converted to an - `lsst.images.psfs.GaussianPointSpreadFunction`. + instance and records the coadds contributed to it, along with the + instrument they were made from. The Gaussian PSF set by `run` is + converted to an `lsst.images.psfs.GaussianPointSpreadFunction`. Parameters ---------- @@ -340,7 +352,10 @@ def convert_outputs_to_future(self, template, coadd_refs, detector, exposureReco """ # The template PSF is a Gaussian, so rebuild the CoaddPsf from the # template inputs to record which coadds contributed. - coaddPsf = self._makePsf(template, template.getInfo().getCoaddInputs().ccds, template.getWcs()) + catalog = template.getInfo().getCoaddInputs().ccds + coaddPsf = self._makePsf(template, catalog, template.getWcs()) + # `run` records the template's instrument on every input. + instrument = catalog[0]["instrument"] template.setDetector(detector) skymaps = {ref.dataId["skymap"] for ref in coadd_refs} if len(skymaps) != 1: @@ -364,6 +379,7 @@ def convert_outputs_to_future(self, template, coadd_refs, detector, exposureReco skymaps.pop(), butler_info, log=self.log, + instrument=instrument, ) return convertedTemplate @@ -477,7 +493,7 @@ def getExposures(self, coaddExposureHandles, bbox, skymap, wcs): return pipeBase.Struct(coaddExposures=coaddExposures, dataIds=dataIds) @timeMethod - def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visit=None): + def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visit=None, instrument=None): """Warp coadds from multiple tracts and patches to form a template to subtract from a science image. @@ -488,6 +504,12 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi The PSF of the template is a Gaussian with the same width as the CoaddPsf evaluated at the center of the detector. + The template is built from coadds of one instrument: the science + image's ``instrument`` if any coadd was made from it, and otherwise the + one that most of the coadds were made from. Coadds made from another + instrument are dropped, and so are coadds that record no instrument + unless none of the coadds records one. + Parameters ---------- coaddExposureHandles : `dict` [`int`, `list` of \ @@ -509,6 +531,9 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi If supplied, over-write the visit ID in the template's visitInfo so that downstream source injection tasks can link the template and science image for the visit. + instrument : `str`, optional + Instrument of the science image. Coadds made from it are used in + preference to coadds from any other instrument. Returns ------- @@ -525,29 +550,29 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi NoWorkFound If no coadds are found with sufficient un-masked pixels. """ - band, photoCalib = self._checkInputs(dataIds, coaddExposureHandles) + band, photoCalib, templateInstrument = self._checkInputs(dataIds, coaddExposureHandles, instrument) bbox = geom.Box2I(bbox) bbox.grow(self.config.templateBorderSize) warped = {} catalogs = [] - # Determine the ratio of the original pixel area to the pixel area of - # the science image, if configured. This will only be different from 1 - # if the coadd was from a different instrument. The ratio can be - # determined once from the components of one coadd exposure, so that no - # pixels have to be read yet. Coadds comprising images from multiple - # different instruments are not supported. plateScaleFactor = None - if self.config.doCorrectVariancePlateScale and coaddExposureHandles: - plateScaleFactor = self._plateScaleFactor( - next(iter(coaddExposureHandles.values()))[0], wcs, bbox - ) for tract in coaddExposureHandles: - maskedImages, catalog, totalBox = self._makeExposureCatalog( - coaddExposureHandles[tract], dataIds[tract] + maskedImages, catalog, totalBox, coaddInputs = self._makeExposureCatalog( + coaddExposureHandles[tract], dataIds[tract], templateInstrument ) + # Determine the ratio of the original pixel area to the pixel area + # of the science image, if configured. This will only be different + # from 1 if the coadd was from a different instrument. One coadd + # stands for all of them; if the first coadd kept in this tract + # cannot give the ratio, the next tract is tried. + if self.config.doCorrectVariancePlateScale and plateScaleFactor is None: + plateScaleFactor = self._plateScaleFactor(coaddInputs, wcs, bbox) + if len(catalog) == 0: + # Every coadd in this tract was from another instrument. + continue warpedBox = computeWarpedBBox(catalog[0].wcs, bbox, wcs) warpedBox.grow(5) # to ensure we catch all relevant input pixels # Combine images from individual patches together. @@ -677,16 +702,14 @@ def _correctVariance(self, template, plateScaleFactor): scale, ) - def _plateScaleFactor(self, coaddHandle, wcs, bbox): + def _plateScaleFactor(self, coaddInputs, wcs, bbox): """Compute the variance correction for the total change in pixel area between the coadd's constituent images and the science image. Parameters ---------- - coaddHandle : `lsst.daf.butler.DeferredDatasetHandle` of \ - `lsst.afw.image.Exposure` - Handle to one of the input coadd patches. Only its ``coaddInputs`` - component is read; the pixels are left alone. + coaddInputs : `lsst.afw.image.CoaddInputs` or `None` + Input records of one of the coadd patches. wcs : `lsst.afw.geom.SkyWcs` WCS of the science image the template is being built for. bbox : `lsst.geom.Box2I` @@ -700,7 +723,6 @@ def _plateScaleFactor(self, coaddHandle, wcs, bbox): from the coadd's instrument to the science instrument, or `None` if this patch carries no input record from which it could be computed. """ - coaddInputs = coaddHandle.get(component="coaddInputs") if coaddInputs is None: return None scienceScale = wcs.getPixelScale(geom.Box2D(bbox).getCenter()).asArcseconds() @@ -744,9 +766,10 @@ def checkHighVariance(self, template): template.mask.array[highVariancePix] |= 2**highVarianceMaskPlaneBit @staticmethod - def _checkInputs(dataIds, coaddExposures): + def _checkInputs(dataIds, coaddExposures, scienceInstrument=None): """Check that the all the dataIds are from the same band and that - the exposures all have the same photometric calibration. + the exposures all have the same photometric calibration, and choose + the instrument of the template. Parameters ---------- @@ -757,6 +780,8 @@ def _checkInputs(dataIds, coaddExposures): `lsst.afw.image.Exposure` or `lsst.afw.image.Exposure`]] Coadds to be mosaicked. + scienceInstrument : `str`, optional + Instrument of the science image. Returns ------- @@ -764,6 +789,11 @@ def _checkInputs(dataIds, coaddExposures): Filter band of all the input exposures. photoCalib : `lsst.afw.image.PhotoCalib` Photometric calibration of all of the input exposures. + instrument : `str` + ``scienceInstrument`` if any coadd was made from it, otherwise the + instrument that the most coadds were made from, or ``"UNKNOWN"`` if + no coadd records one. A tie goes to the instrument whose name sorts + first. Raises ------ @@ -784,9 +814,48 @@ def _checkInputs(dataIds, coaddExposures): msg = f"GetTemplateTask called with exposures with different photoCalibs: {photoCalibs}" raise RuntimeError(msg) photoCalib = photoCalibs[0] - return band, photoCalib + instrumentCounts = collections.Counter( + GetTemplateTask._getInstrument(exposure.get(component="coaddInputs")) + for exposures in coaddExposures.values() + for exposure in exposures + ) + # Drop instrument labels that are empty or set to UNKNOWN from the sort + del instrumentCounts[""] + del instrumentCounts["UNKNOWN"] + if scienceInstrument in instrumentCounts: + return band, photoCalib, scienceInstrument + if not instrumentCounts: + return band, photoCalib, "UNKNOWN" + instrument = max(sorted(instrumentCounts), key=instrumentCounts.get) + return band, photoCalib, instrument - def _makeExposureCatalog(self, exposureRefs, dataIds): + @staticmethod + def _getInstrument(coaddInputs): + """Return the instrument whose observations were coadded to make a + coadd. + + Parameters + ---------- + coaddInputs : `lsst.afw.image.CoaddInputs` or `None` + Input records of the coadd. + + Returns + ------- + instrument : `str` + Instrument label of the coadd, or ``"UNKNOWN"`` if it records none. + + Notes + ----- + Every input of a coadd comes from the same instrument, so the first + input record stands for all of them. + """ + if coaddInputs is None or len(coaddInputs.ccds) == 0: + return _UNKNOWN_INSTRUMENT + visitInfo = coaddInputs.ccds[0].getVisitInfo() + label = visitInfo.getInstrumentLabel() if visitInfo is not None else "" + return label or _UNKNOWN_INSTRUMENT + + def _makeExposureCatalog(self, exposureRefs, dataIds, instrument=_UNKNOWN_INSTRUMENT): """Make an exposure catalog for one tract. Parameters @@ -797,6 +866,10 @@ def _makeExposureCatalog(self, exposureRefs, dataIds): dataIds : `list` [`lsst.daf.butler.DataCoordinate`] Data ids of each of the included exposures; must have "tract" and "patch" entries. + instrument : `str`, optional + Instrument of the template, recorded in the catalog. Exposures + made from another instrument are dropped, as are exposures that + record no instrument unless ``instrument`` is ``"UNKNOWN"``. Returns ------- @@ -806,14 +879,30 @@ def _makeExposureCatalog(self, exposureRefs, dataIds): Catalog of metadata for each exposure totalBox : `lsst.geom.Box2I` The union of the bounding boxes of all the input exposures. + coaddInputs : `lsst.afw.image.CoaddInputs` or `None` + Input records of the first exposure kept, or `None` if none is + kept. """ catalog = afwTable.ExposureCatalog(self.schema) catalog.reserve(len(exposureRefs)) exposures = (exposureRef.get() for exposureRef in exposureRefs) images = {} totalBox = geom.Box2I() + firstCoaddInputs = None for coadd, dataId in zip(exposures, dataIds): + coaddInputs = coadd.getInfo().getCoaddInputs() + coaddInstrument = self._getInstrument(coaddInputs) + if coaddInstrument != instrument: + self.log.info( + "Dropping template coadd %s: it was made from %s, but the template is made from %s.", + dataId, + coaddInstrument, + instrument, + ) + continue + if not images: + firstCoaddInputs = coaddInputs images[dataId] = coadd.maskedImage bbox = coadd.getBBox() totalBox = totalBox.expandedTo(bbox) @@ -825,11 +914,12 @@ def _makeExposureCatalog(self, exposureRefs, dataIds): record.setValidPolygon(afwGeom.Polygon(geom.Box2D(bbox).getCorners())) record.set("tract", dataId["tract"]) record.set("patch", dataId["patch"]) + record.set("instrument", instrument) # Weight is used by CoaddPsf, but the PSFs from overlapping patches # should be very similar, so this value mostly shouldn't matter. record.set("weight", 1) - return images, catalog, totalBox + return images, catalog, totalBox, firstCoaddInputs def _merge(self, maskedImages, bbox, wcs): """Merge the images that came from one tract into one larger image, @@ -1070,6 +1160,8 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): wcs=wcs, dataIds=results.dataIds, physical_filter=physical_filter, + visit=outputRefs.template.dataId["visit"], + instrument=butlerQC.quantum.dataId["instrument"], ) butlerQC.put(outputs, outputRefs) diff --git a/tests/test_getTemplate.py b/tests/test_getTemplate.py index 1bb0a81b..5e1496bf 100644 --- a/tests/test_getTemplate.py +++ b/tests/test_getTemplate.py @@ -438,7 +438,7 @@ def testNanInputs(self, box=None, nInput=None): self.assertTrue(np.isfinite(result.template.image.array).all()) def _runCorrection(self, doCorrectVariancePlateScale=False, doScaleVariance=False, - handles=None): + handles=None, instrument=None): """Run the task on tract 0 with the named variance corrections. Everything defaults to off, so each test turns on exactly what it is @@ -452,7 +452,8 @@ def _runCorrection(self, doCorrectVariancePlateScale=False, doScaleVariance=Fals bbox=lsst.geom.Box2I(self.varianceBox), wcs=self.exposure.wcs, dataIds={0: self.dataIds[0]}, - physical_filter="a_test") + physical_filter="a_test", + instrument=instrument) return task, result.template def _makeScaledWcs(self, factor): @@ -466,7 +467,7 @@ def _makeScaledWcs(self, factor): cdMatrix) @staticmethod - def _makeCoaddInputs(records): + def _makeCoaddInputs(records, instrument="TestCam"): """Make a CoaddInputs holding the given input records. Parameters @@ -476,6 +477,10 @@ def _makeCoaddInputs(records): The WCS and bbox to record for each input. A `None` WCS makes a record that `_plateScaleFactor` has to skip. May be empty, to simulate a coadd whose plate scale cannot be reconstructed. + instrument : `str` or `None`, optional + Instrument label to record for every input. `None` leaves the + inputs without a visit info, and an empty string gives them a + visit info with no label; either way the instrument is unknown. """ ccdSchema = lsst.afw.table.ExposureTable.makeMinimalSchema() weightKey = ccdSchema.addField("weight", type=float, doc="Coadd weight") @@ -485,19 +490,21 @@ def _makeCoaddInputs(records): record = coaddInputs.ccds.addNew() record.setWcs(wcs) record.setBBox(bbox) + if instrument is not None: + record.setVisitInfo(lsst.afw.image.VisitInfo(instrumentLabel=instrument)) # Included because real coadds have it, though a single-record # correction does not use it. record.set(weightKey, 1.0) return coaddInputs - def _patchHandles(self, tract, records): + def _patchHandles(self, tract, records, instrument="TestCam"): """Return handles for a tract's patches, with their CoaddInputs - replaced by ``records``. + replaced by ``records`` from ``instrument``. """ handles = [] for ref in self.patches[tract]: coadd = ref.get() - coadd.getInfo().setCoaddInputs(self._makeCoaddInputs(records)) + coadd.getInfo().setCoaddInputs(self._makeCoaddInputs(records, instrument=instrument)) handles.append(pipeBase.InMemoryDatasetHandle( coadd, storageClass="ExposureF", copy=True, dataId=ref.dataId)) return handles @@ -543,6 +550,33 @@ def testCorrectVariancePlateScaleUsesOneRecord(self): # last (which would give 1.0) is never reached. self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 4.0, rtol=1e-6) + def testCorrectVariancePlateScaleSkipsDroppedCoadd(self): + """The plate scale comes from the first coadd that is used, not from + one dropped for its instrument. + """ + bbox = self.exposure.getBBox() + handles = (self._patchHandles(0, [(self._makeScaledWcs(2.0), bbox)], instrument="OtherCam")[:1] + + self._patchHandles(0, [(self._makeScaledWcs(1.315), bbox)])[1:]) + task, _ = self._runCorrection(doCorrectVariancePlateScale=True, handles=handles, + instrument="TestCam") + self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 1.315**2, rtol=1e-6) + + def testCorrectVariancePlateScaleNextTract(self): + """If the first coadd kept in a tract cannot give the plate scale, the + first coadd kept in the next tract is used. + """ + config = lsst.ip.diffim.GetTemplateTask.ConfigClass() + config.doCorrectVariancePlateScale = True + task = lsst.ip.diffim.GetTemplateTask(config=config) + handles = {0: self._patchHandles(0, []), + 1: self._patchHandles(1, [(self._makeScaledWcs(1.315), self.exposure.getBBox())])} + task.run(coaddExposureHandles=handles, + bbox=lsst.geom.Box2I(self.varianceBox), + wcs=self.exposure.wcs, + dataIds={0: self.dataIds[0], 1: self.dataIds[1]}, + physical_filter="a_test") + self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 1.315**2, rtol=1e-6) + def testCorrectVariancePlateScaleNeedsCoaddInputs(self): """Without usable coadd inputs the plate scale cannot be reconstructed, and the task must say so rather than silently applying @@ -552,6 +586,150 @@ def testCorrectVariancePlateScaleNeedsCoaddInputs(self): with self.assertRaisesRegex(RuntimeError, "doCorrectVariancePlateScale"): self._runCorrection(doCorrectVariancePlateScale=True, handles=handles) + def _runTract0(self, handles): + """Build a template from the given tract 0 patch handles. + """ + box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) + task = lsst.ip.diffim.GetTemplateTask() + return task.run(coaddExposureHandles={0: handles}, + bbox=box, + wcs=self.exposure.wcs, + dataIds={0: self.dataIds[0]}, + physical_filter="a_test") + + def testRunRecordsInstrument(self): + """The instrument of each coadd is recorded with the template's coadd + inputs. + """ + result = self._runTract0(self.patches[0]) + ccds = result.template.getInfo().getCoaddInputs().ccds + self.assertEqual(len(ccds), 4) + self.assertEqual({record["instrument"] for record in ccds}, {"TestCam"}) + + def testRunUnknownInstrument(self): + """A coadd whose inputs record no instrument is dropped when the + other coadds record one, and kept, with the template's instrument + ``"UNKNOWN"``, when none of them records one. + """ + records = [(self.exposure.wcs, self.exposure.getBBox())] + # Inputs with no visit info, and with a visit info but no label. + for unknown in (None, ""): + with self.subTest(label=unknown): + handles = (self._patchHandles(0, records, instrument=unknown)[:1] + + self._patchHandles(0, records)[1:]) + ccds, dropped = self._runLogged({0: handles}) + self.assertEqual(len(ccds), 3) + self.assertEqual({record["instrument"] for record in ccds}, {"TestCam"}) + self.assertEqual(len(dropped), 1) + self.assertIn("made from UNKNOWN", dropped[0]) + + # No coadd records an instrument. + result = self._runTract0(self._patchHandles(0, records, instrument=unknown)) + ccds = result.template.getInfo().getCoaddInputs().ccds + self.assertEqual(len(ccds), 4) + self.assertEqual({record["instrument"] for record in ccds}, {"UNKNOWN"}) + + def _runLogged(self, handles, dataIds=None, box=None, instrument=None): + """Build a template, capturing the messages about dropped coadds. + + Parameters + ---------- + handles : `dict` [`int`, `list`] + Patch handles, keyed by tract. + dataIds : `dict` [`int`, `list`], optional + Data ids matching ``handles``; defaults to the fixture's. + box : `lsst.geom.Box2I`, optional + Box to build the template on; defaults to one in tract 0. + instrument : `str`, optional + Instrument of the science image. + + Returns + ------- + ccds : `lsst.afw.table.ExposureCatalog` + The coadd inputs recorded on the template. + dropped : `list` [`str`] + The info messages about dropped coadds. + """ + if dataIds is None: + dataIds = {tract: self.dataIds[tract] for tract in handles} + if box is None: + box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) + task = lsst.ip.diffim.GetTemplateTask() + with self.assertLogs(task.log.name, level="INFO") as cm: + result = task.run(coaddExposureHandles=handles, + bbox=box, + wcs=self.exposure.wcs, + dataIds=dataIds, + physical_filter="a_test", + instrument=instrument) + dropped = [record for record in cm.records + if record.getMessage().startswith("Dropping template coadd")] + self.assertEqual({record.levelname for record in dropped}, {"INFO"} if dropped else set()) + return result.template.getInfo().getCoaddInputs().ccds, [record.getMessage() for record in dropped] + + def testRunDropsMinorityInstrument(self): + """Coadds made from an instrument other than the one most coadds were + made from are dropped. + """ + records = [(self.exposure.wcs, self.exposure.getBBox())] + handles = (self._patchHandles(0, records)[:3] + + self._patchHandles(0, records, instrument="OtherCam")[3:]) + + ccds, dropped = self._runLogged({0: handles}) + self.assertEqual({record["patch"] for record in ccds}, {0, 1, 2}) + self.assertEqual({record["instrument"] for record in ccds}, {"TestCam"}) + self.assertEqual(len(dropped), 1) + self.assertIn(str(self.dataIds[0][3]), dropped[0]) + self.assertIn("made from OtherCam", dropped[0]) + + def testRunPrefersScienceInstrument(self): + """Coadds made from the science image's instrument are used even when + they are in the minority; otherwise the majority instrument is used. + """ + records = [(self.exposure.wcs, self.exposure.getBBox())] + handles = (self._patchHandles(0, records, instrument="OtherCam")[:3] + + self._patchHandles(0, records)[3:]) + + ccds, dropped = self._runLogged({0: handles}, instrument="TestCam") + self.assertEqual({record["patch"] for record in ccds}, {3}) + self.assertEqual({record["instrument"] for record in ccds}, {"TestCam"}) + self.assertEqual(len(dropped), 3) + self.assertTrue(all("made from OtherCam" in message for message in dropped)) + + # No coadd was made from the science image's instrument. + ccds, dropped = self._runLogged({0: handles}, instrument="SciCam") + self.assertEqual({record["patch"] for record in ccds}, {0, 1, 2}) + self.assertEqual({record["instrument"] for record in ccds}, {"OtherCam"}) + self.assertEqual(len(dropped), 1) + + def testRunInstrumentTie(self): + """A tie goes to the instrument whose name sorts first. + """ + records = [(self.exposure.wcs, self.exposure.getBBox())] + handles = (self._patchHandles(0, records)[:2] + + self._patchHandles(0, records, instrument="OtherCam")[2:]) + + ccds, dropped = self._runLogged({0: handles}) + self.assertEqual({record["patch"] for record in ccds}, {2, 3}) + self.assertEqual({record["instrument"] for record in ccds}, {"OtherCam"}) + self.assertEqual(len(dropped), 2) + + def testRunDropsNonOverlappingCoadd(self): + """A coadd from another instrument is dropped even if it does not + overlap the template. + """ + records = [(self.exposure.wcs, self.exposure.getBBox())] + # Tract 3 does not overlap a box in the lower left of tract 0. + box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(50, 50)) + handles = {0: self.patches[0], 3: self._patchHandles(3, records, instrument="OtherCam")[:1]} + dataIds = {0: self.dataIds[0], 3: self.dataIds[3][:1]} + + ccds, dropped = self._runLogged(handles, dataIds=dataIds, box=box) + self.assertEqual({record["tract"] for record in ccds}, {0}) + self.assertEqual({record["instrument"] for record in ccds}, {"TestCam"}) + self.assertEqual(len(dropped), 1) + self.assertIn(str(self.dataIds[3][0]), dropped[0]) + def _scaleInputVariance(self, tract, factor): """Return fresh handles for one tract's patches, with their variance planes multiplied by ``factor``. @@ -639,7 +817,7 @@ def _configureAndRunTask(doScaleVariance, varianceScale=1.): factor*scaleFactor, rtol=1e-5) self.assertImagesAlmostEqual(templateLow.variance, templateOn.variance, rtol=1e-5) - def _runLegacyForFuture(self, raiseOnUndefinedMaskMap=True): + def _runLegacyForFuture(self, raiseOnUndefinedMaskMap=True, handles=None): """Build a template from tract 0 with a task configured for the future output type, without converting it. @@ -649,6 +827,8 @@ def _runLegacyForFuture(self, raiseOnUndefinedMaskMap=True): ---------- raiseOnUndefinedMaskMap : `bool`, optional Value of the task config field of the same name. + handles : `list` [`lsst.pipe.base.InMemoryDatasetHandle`], optional + Handles to the tract 0 patches to use instead of the fixture's. Returns ------- @@ -664,7 +844,7 @@ def _runLegacyForFuture(self, raiseOnUndefinedMaskMap=True): task = lsst.ip.diffim.GetTemplateTask(config=config) self.visit = 9876 box = lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)) - result = task.run(coaddExposureHandles={0: self.patches[0]}, + result = task.run(coaddExposureHandles={0: self.patches[0] if handles is None else handles}, bbox=box, wcs=self.exposure.wcs, dataIds={0: self.dataIds[0]}, @@ -675,9 +855,14 @@ def _runLegacyForFuture(self, raiseOnUndefinedMaskMap=True): self.futureTask = task return result, box - def _runFuture(self): + def _runFuture(self, handles=None): """Build a template and convert it, the way ``runQuantum`` does. + Parameters + ---------- + handles : `list` [`lsst.pipe.base.InMemoryDatasetHandle`], optional + Handles to the tract 0 patches to use instead of the fixture's. + Returns ------- result : `lsst.pipe.base.Struct` @@ -689,7 +874,7 @@ def _runFuture(self): # lsst.images needs per-amplifier raw geometry and a field angle # transform, which the trivial test detectors do not have. detector = list(lsst.afw.cameraGeom.testUtils.CameraWrapper().camera)[0] - result, box = self._runLegacyForFuture() + result, box = self._runLegacyForFuture(handles=handles) result.template = self.futureTask.convert_outputs_to_future( result.template, self._coaddRefs(0), detector=detector, exposureRecord=self._exposureRecord()) @@ -851,10 +1036,24 @@ def testConvertOutputsToFutureTemplates(self): self.assertEqual(template.skymap, "skymap") self.assertEqual(template.dataset_id, ref.id) self.assertEqual(template.dataset_run, ref.run) + self.assertEqual(template.instrument, "TestCam") self.assertFalse(template.psf_shape_flag) self.assertGreater(template.psf_shape_xx, 0) self.assertGreater(template.bounds.area, 0) + def testConvertOutputsToFutureInstrument(self): + """The template instrument is the one recorded by the coadds that + know it, or ``"UNKNOWN"`` if none of them do. + """ + records = [(self.exposure.wcs, self.exposure.getBBox())] + handles = (self._patchHandles(0, records, instrument=None)[:1] + + self._patchHandles(0, records)[1:]) + result, _, _ = self._runFuture(handles=handles) + self.assertEqual({template.instrument for template in result.template.templates}, {"TestCam"}) + + result, _, _ = self._runFuture(handles=self._patchHandles(0, records, instrument=None)) + self.assertEqual({template.instrument for template in result.template.templates}, {"UNKNOWN"}) + def testConvertOutputsToFuturePsf(self): """The Gaussian PSF set by `run` is kept, defined over the detector. """ From 6ae34e9a5d4ee75fc39cdd894212c4866a2271d4 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Fri, 16 Oct 2026 12:30:52 -0700 Subject: [PATCH 2/3] Remove the doCorrectVariancePlateScale config The template variance is now corrected for the change in pixel area whenever the template coadds were made from a different instrument than the science image, and a factor that cannot be determined logs a warning instead of raising. The field never appeared in a versioned release, only in weeklies w.2026.38 to w.2026.42, so it is removed without a deprecation period. --- python/lsst/ip/diffim/getTemplate.py | 66 +++++++++++----------------- tests/test_getTemplate.py | 46 +++++++++---------- 2 files changed, 48 insertions(+), 64 deletions(-) diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index 76290066..d0846443 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -149,16 +149,6 @@ class GetTemplateConfig( dtype=afwMath.Warper.ConfigClass, doc="warper configuration", ) - doCorrectVariancePlateScale = pexConfig.Field( - dtype=bool, - default=False, - doc="Correct the template variance for the total change in pixel" - " area between the images the coadds were built from and the science" - " image? This matters when the two have different plate scales, such" - " as DECam templates used for Rubin science images. The coadd pixel" - " grid cancels out of the total, so the skymap does not enter. It is a" - " single scalar factor, applied to the final template variance plane.", - ) coaddPsf = pexConfig.ConfigField( doc="Configuration for CoaddPsf", dtype=CoaddPsfConfig, @@ -549,6 +539,10 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi ------ NoWorkFound If no coadds are found with sufficient un-masked pixels. + RuntimeError + If the template is not known to be from the science image's + instrument and no coadd records the plate scale of its input + images, so the template variance cannot be corrected. """ band, photoCalib, templateInstrument = self._checkInputs(dataIds, coaddExposureHandles, instrument) @@ -557,22 +551,21 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi warped = {} catalogs = [] - plateScaleFactor = None + plateScaleFactor = 1.0 if instrument == templateInstrument else None for tract in coaddExposureHandles: maskedImages, catalog, totalBox, coaddInputs = self._makeExposureCatalog( coaddExposureHandles[tract], dataIds[tract], templateInstrument ) - # Determine the ratio of the original pixel area to the pixel area - # of the science image, if configured. This will only be different - # from 1 if the coadd was from a different instrument. One coadd - # stands for all of them; if the first coadd kept in this tract - # cannot give the ratio, the next tract is tried. - if self.config.doCorrectVariancePlateScale and plateScaleFactor is None: - plateScaleFactor = self._plateScaleFactor(coaddInputs, wcs, bbox) if len(catalog) == 0: - # Every coadd in this tract was from another instrument. + # Every coadd in this tract was dropped for its instrument. continue + # Determine the ratio of the original pixel area to the pixel area + # of the science image. This will only be different from 1 if the + # coadd was from a different instrument. A single correction factor + # is calculated from the first usable coadd and applied to all. + if plateScaleFactor is None: + plateScaleFactor = self._plateScaleFactor(coaddInputs, wcs, bbox) warpedBox = computeWarpedBBox(catalog[0].wcs, bbox, wcs) warpedBox.grow(5) # to ensure we catch all relevant input pixels # Combine images from individual patches together. @@ -657,8 +650,8 @@ def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visi def _correctVariance(self, template, plateScaleFactor): """Correct the template variance plane for the effects of warping. - Compute a multiplicative correction factor for the template variance - if the plate scale of the coadd's constituent images is different than + Apply a multiplicative correction factor to the template variance if + the plate scale of the coadd's constituent images is different than the science image the template is being constructed for. This should only be necessary if the coadd images were from a different instrument than the science image. If the plate scale of the science instrument is @@ -676,30 +669,21 @@ def _correctVariance(self, template, plateScaleFactor): Raises ------ RuntimeError - Raised if ``doCorrectVariancePlateScale`` is set but the factor - could not be reconstructed from the coadd inputs. + Raised if ``plateScaleFactor`` is `None`. """ - scale = 1.0 - if self.config.doCorrectVariancePlateScale: - if plateScaleFactor is None: - raise RuntimeError( - "doCorrectVariancePlateScale is set but no usable coaddInputs were found." - ) - scale *= plateScaleFactor - self.metadata["variancePlateScaleFactor"] = plateScaleFactor - self.log.info( - "Applying a plate scale variance factor of %.4f, reconstructed from the" - " coadd inputs.", - plateScaleFactor, + if plateScaleFactor is None: + raise RuntimeError( + "Could not determine the plate scale of the images the template coadds were made from," + " so the template variance cannot be corrected for the change in pixel area." ) - - if scale == 1.0: + self.metadata["variancePlateScaleFactor"] = plateScaleFactor + if plateScaleFactor == 1.0: return - template.variance.array *= scale - self.metadata["templateVarianceCorrectionFactor"] = scale + template.variance.array *= plateScaleFactor + self.metadata["templateVarianceCorrectionFactor"] = plateScaleFactor self.log.info( - "Corrected the template variance plane for warping by a factor of %.4f.", - scale, + "Corrected the template variance plane for the change in pixel area by a factor of %.4f.", + plateScaleFactor, ) def _plateScaleFactor(self, coaddInputs, wcs, bbox): diff --git a/tests/test_getTemplate.py b/tests/test_getTemplate.py index 5e1496bf..9ae6bee0 100644 --- a/tests/test_getTemplate.py +++ b/tests/test_getTemplate.py @@ -437,16 +437,16 @@ def testNanInputs(self, box=None, nInput=None): # in the template are closer to the original anymore. self.assertTrue(np.isfinite(result.template.image.array).all()) - def _runCorrection(self, doCorrectVariancePlateScale=False, doScaleVariance=False, - handles=None, instrument=None): - """Run the task on tract 0 with the named variance corrections. + def _runCorrection(self, handles=None, instrument="SciCam"): + """Run the task on tract 0, with a science image from ``instrument`` + and variance scaling turned off. - Everything defaults to off, so each test turns on exactly what it is - exercising. + The plate scale correction is applied only when ``instrument`` is not + the instrument of the template coadds, which is ``"TestCam"`` unless + ``handles`` says otherwise. """ config = lsst.ip.diffim.GetTemplateTask.ConfigClass() - config.doCorrectVariancePlateScale = doCorrectVariancePlateScale - config.doScaleVariance = doScaleVariance + config.doScaleVariance = False task = lsst.ip.diffim.GetTemplateTask(config=config) result = task.run(coaddExposureHandles={0: handles or self.patches[0]}, bbox=lsst.geom.Box2I(self.varianceBox), @@ -513,11 +513,15 @@ def testCorrectVariancePlateScale(self): """The plate scale correction is the total pixel area change from the images the coadds were built from to the science image. """ - _, off = self._runCorrection() + # The science image and the coadds are from the same instrument, so + # the factor is 1 without reading any coadd inputs. + task, off = self._runCorrection(instrument="TestCam") + self.assertEqual(task.metadata["variancePlateScaleFactor"], 1.0) + self.assertNotIn("templateVarianceCorrectionFactor", task.metadata) # The fixture's records are the science image itself, so there is no - # net change in pixel area: the same-instrument case. - task, on = self._runCorrection(doCorrectVariancePlateScale=True) + # net change in pixel area even when the instruments differ. + task, on = self._runCorrection() self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 1.0, rtol=1e-6) # Coarser original pixels than science pixels, the DECam-template @@ -526,8 +530,7 @@ def testCorrectVariancePlateScale(self): with self.subTest(scaleFactor=scaleFactor): handles = self._patchHandles( 0, [(self._makeScaledWcs(scaleFactor), self.exposure.getBBox())]) - task, on = self._runCorrection(doCorrectVariancePlateScale=True, - handles=handles) + task, on = self._runCorrection(handles=handles) factor = task.metadata["variancePlateScaleFactor"] self.assertFloatsAlmostEqual(factor, scaleFactor**2, rtol=1e-6) self.assertFloatsAlmostEqual(on.variance.array, off.variance.array*factor, @@ -544,7 +547,7 @@ def testCorrectVariancePlateScaleUsesOneRecord(self): 0, [(None, self.exposure.getBBox()), (self._makeScaledWcs(2.0), self.exposure.getBBox()), (self.exposure.wcs, self.exposure.getBBox())]) - task, _ = self._runCorrection(doCorrectVariancePlateScale=True, handles=handles) + task, _ = self._runCorrection(handles=handles) # The first record with a WCS: the one without is skipped, and the # last (which would give 1.0) is never reached. @@ -557,34 +560,31 @@ def testCorrectVariancePlateScaleSkipsDroppedCoadd(self): bbox = self.exposure.getBBox() handles = (self._patchHandles(0, [(self._makeScaledWcs(2.0), bbox)], instrument="OtherCam")[:1] + self._patchHandles(0, [(self._makeScaledWcs(1.315), bbox)])[1:]) - task, _ = self._runCorrection(doCorrectVariancePlateScale=True, handles=handles, - instrument="TestCam") + task, _ = self._runCorrection(handles=handles) self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 1.315**2, rtol=1e-6) def testCorrectVariancePlateScaleNextTract(self): """If the first coadd kept in a tract cannot give the plate scale, the first coadd kept in the next tract is used. """ - config = lsst.ip.diffim.GetTemplateTask.ConfigClass() - config.doCorrectVariancePlateScale = True - task = lsst.ip.diffim.GetTemplateTask(config=config) + task = lsst.ip.diffim.GetTemplateTask() handles = {0: self._patchHandles(0, []), 1: self._patchHandles(1, [(self._makeScaledWcs(1.315), self.exposure.getBBox())])} task.run(coaddExposureHandles=handles, bbox=lsst.geom.Box2I(self.varianceBox), wcs=self.exposure.wcs, dataIds={0: self.dataIds[0], 1: self.dataIds[1]}, - physical_filter="a_test") + physical_filter="a_test", + instrument="SciCam") self.assertFloatsAlmostEqual(task.metadata["variancePlateScaleFactor"], 1.315**2, rtol=1e-6) def testCorrectVariancePlateScaleNeedsCoaddInputs(self): """Without usable coadd inputs the plate scale cannot be - reconstructed, and the task must say so rather than silently applying - only part of the correction. + reconstructed, so the template cannot be built. """ handles = self._patchHandles(0, []) - with self.assertRaisesRegex(RuntimeError, "doCorrectVariancePlateScale"): - self._runCorrection(doCorrectVariancePlateScale=True, handles=handles) + with self.assertRaisesRegex(RuntimeError, "cannot be corrected for the change in pixel area"): + self._runCorrection(handles=handles) def _runTract0(self, handles): """Build a template from the given tract 0 patch handles. From 0e7636f01e040d25567718136cb01b45c20ecbc8 Mon Sep 17 00:00:00 2001 From: Ian Sullivan Date: Fri, 16 Oct 2026 17:27:51 -0700 Subject: [PATCH 3/3] Read photoCalib and coaddInputs in one loop --- python/lsst/ip/diffim/getTemplate.py | 25 +++++++++++-------------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py index d0846443..0bf168c4 100644 --- a/python/lsst/ip/diffim/getTemplate.py +++ b/python/lsst/ip/diffim/getTemplate.py @@ -789,27 +789,24 @@ def _checkInputs(dataIds, coaddExposures, scienceInstrument=None): if len(bands) > 1: raise RuntimeError(f"GetTemplateTask called with multiple bands: {bands}") band = bands.pop() - photoCalibs = [ - exposure.get(component="photoCalib") - for exposures in coaddExposures.values() - for exposure in exposures - ] + photoCalibs = [] + instrumentCounts = collections.Counter() + for exposures in coaddExposures.values(): + for exposure in exposures: + # Read both components of a coadd in turn, so a butler cache + # can serve the second read from the first download. + photoCalibs.append(exposure.get(component="photoCalib")) + instrumentCounts[GetTemplateTask._getInstrument(exposure.get(component="coaddInputs"))] += 1 if not all([photoCalibs[0] == x for x in photoCalibs]): msg = f"GetTemplateTask called with exposures with different photoCalibs: {photoCalibs}" raise RuntimeError(msg) photoCalib = photoCalibs[0] - instrumentCounts = collections.Counter( - GetTemplateTask._getInstrument(exposure.get(component="coaddInputs")) - for exposures in coaddExposures.values() - for exposure in exposures - ) - # Drop instrument labels that are empty or set to UNKNOWN from the sort - del instrumentCounts[""] - del instrumentCounts["UNKNOWN"] + # Coadds that record no instrument do not count toward the choice. + del instrumentCounts[_UNKNOWN_INSTRUMENT] if scienceInstrument in instrumentCounts: return band, photoCalib, scienceInstrument if not instrumentCounts: - return band, photoCalib, "UNKNOWN" + return band, photoCalib, _UNKNOWN_INSTRUMENT instrument = max(sorted(instrumentCounts), key=instrumentCounts.get) return band, photoCalib, instrument