Skip to content

DM-56663: record the template instrument in GetTemplateTask #476

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
isullivan merged 3 commits into from
Oct 20, 2026
Merged

DM-56663: record the template instrument in GetTemplateTask #476

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
0994555
Record the template instrument in GetTemplateTask
isullivan Oct 16, 2026
6ae34e9
Remove the doCorrectVariancePlateScale config
isullivan Oct 16, 2026
0e7636f
Read photoCalib and coaddInputs in one loop
isullivan Oct 17, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
207 changes: 140 additions & 67 deletions python/lsst/ip/diffim/getTemplate.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -144,16 +149,6 @@ class GetTemplateConfig(
dtype=afwMath.Warper.ConfigClass,
doc="warper configuration",
)
doCorrectVariancePlateScale = pexConfig.Field(
Comment thread
isullivan marked this conversation as resolved.
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,
Expand Down Expand Up @@ -267,6 +262,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):
Expand All @@ -290,6 +291,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
Expand All @@ -305,9 +307,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
----------
Expand Down Expand Up @@ -340,7 +342,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:
Expand All @@ -364,6 +369,7 @@ def convert_outputs_to_future(self, template, coadd_refs, detector, exposureReco
skymaps.pop(),
butler_info,
log=self.log,
instrument=instrument,
Comment thread
isullivan marked this conversation as resolved.
)
return convertedTemplate

Expand Down Expand Up @@ -477,7 +483,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.

Expand All @@ -488,6 +494,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 \
Expand All @@ -509,6 +521,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
-------
Expand All @@ -524,30 +539,33 @@ 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 = 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
)
plateScaleFactor = 1.0 if instrument == templateInstrument else None

for tract in coaddExposureHandles:
maskedImages, catalog, totalBox = self._makeExposureCatalog(
coaddExposureHandles[tract], dataIds[tract]
maskedImages, catalog, totalBox, coaddInputs = self._makeExposureCatalog(
coaddExposureHandles[tract], dataIds[tract], templateInstrument
)
if len(catalog) == 0:
# 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.
Expand Down Expand Up @@ -632,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
Expand All @@ -651,42 +669,31 @@ 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, 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.
Comment thread
isullivan marked this conversation as resolved.
wcs : `lsst.afw.geom.SkyWcs`
WCS of the science image the template is being built for.
bbox : `lsst.geom.Box2I`
Expand All @@ -700,7 +707,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()
Expand Down Expand Up @@ -744,9 +750,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
----------
Expand All @@ -757,13 +764,20 @@ 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
-------
band : `str`
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
------
Expand All @@ -775,18 +789,54 @@ def _checkInputs(dataIds, coaddExposures):
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]
return band, photoCalib
# 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_INSTRUMENT
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
Expand All @@ -797,6 +847,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
-------
Expand All @@ -806,14 +860,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)
Expand All @@ -825,11 +895,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,
Expand Down Expand Up @@ -1070,6 +1141,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)

Expand Down
Loading
Loading