getTemplate.py

Go to the documentation of this file.

# This file is part of ip_diffim.
#
# Developed for the LSST Data Management System.
# This product includes software developed by the LSST Project
# (https://www.lsst.org).
# See the COPYRIGHT file at the top-level directory of this distribution
# for details of code ownership.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.
import collections
import dataclasses
import functools
 
import astropy.units
from astropy.stats import gaussian_sigma_to_fwhm
import numpy as np
 
import lsst.afw.detection as afwDetection
import lsst.afw.image as afwImage
import lsst.geom as geom
import lsst.afw.geom as afwGeom
from lsst.afw.image import VisitInfo
import lsst.afw.table as afwTable
from lsst.afw.math._warper import computeWarpedBBox
import lsst.afw.math as afwMath
from lsst.daf.base import DateTime
import lsst.pex.config as pexConfig
import lsst.pex.exceptions
import lsst.pipe.base as pipeBase
 
from lsst.images import (
    Box,
    DifferenceImage,
    DifferenceImageTemplateInfo,
    DifferenceImageTemplateInput,
    GeneralFrame,
    NoOverlapError,
    Polygon,
    Transform,
    get_legacy_optional_mask_planes,
    get_legacy_template_mask_planes,
)
from lsst.skymap import BaseSkyMap
from lsst.ip.diffim.dcrModel import DcrModel
from lsst.ip.diffim.utils import evaluateMeanPsfFwhm, getPsfFwhm, record_from_obs_info
from lsst.meas.algorithms import CoaddPsf, CoaddPsfConfig, SubtractBackgroundTask, ScaleVarianceTask
from lsst.utils.timer import timeMethod
 
__all__ = [
    "GetTemplateTask",
    "GetTemplateConfig",
    "GetDcrTemplateTask",
    "GetDcrTemplateConfig",
]
 
_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.
"""
 
_FOOTPRINT_MARGIN = 32.0
"""Distance, in coadd pixels, by which the science image's footprint in a
coadd is grown before a coadd input is clipped to it.  It must exceed twice
the distance by which a distorted science image edge departs from a straight
line between its corners, and limits how far outside the science image its
WCS is evaluated.
"""
 
 
def _pixelTransform(fromWcs, toWcs):
    """Return a chainable transformation between two pixel frames.
 
    Parameters
    ----------
    fromWcs, toWcs : `lsst.afw.geom.SkyWcs`
        WCSs of the input and output pixels.
 
    Returns
    -------
    transform : `lsst.images.Transform`
        Transform from ``fromWcs`` pixels to ``toWcs`` pixels.
    """
    return Transform.from_legacy(afwGeom.makeWcsPairTransform(fromWcs, toWcs),
                                 GeneralFrame(unit=astropy.units.pix),
                                 GeneralFrame(unit=astropy.units.pix)
                                 )
 
 
def _getImageOverlap(record, coaddWcs, scienceFootprint, coaddToScienceTransform, scienceBox, border):
    """Return the region enclosing the overlap of the science image and one
    input image of the coadd.
 
    Note that there are three different WCS's involved here: the WCS of the
    original image that was warped to contribute to the coadd, the coadd WCS,
    and the current science image WCS that we are building a template for.
 
    Parameters
    ----------
    record : `lsst.afw.table.ExposureRecord`
        Input record of a coadd, with the input's WCS, bounding box and
        optional valid polygon.
    coaddWcs : `lsst.afw.geom.SkyWcs`
        WCS of the coadd.
    scienceFootprint : `lsst.images.Region`
        ``scienceBox`` in coadd pixels, grown by ``_FOOTPRINT_MARGIN``.
    coaddToScienceTransform : `lsst.images.Transform`
        Transform from coadd pixels to science image pixels.
    scienceBox : `lsst.images.Box`
        Bounding box of the science image.
    border : `float`
        Width, in input pixels, of the border of the input that warping left
        out of the coadd.
 
    Returns
    -------
    footprint : `lsst.images.Polygon` or `None`
        The part of ``scienceBox`` the input's warp covers, in science image
        pixels, or `None` if it covers none of it.
    """
    inputBox = geom.Box2D(record.getBBox())
    inputBox.grow(-border)
    inputBounds = Polygon.from_legacy(afwGeom.Polygon(inputBox))
    polygon = record.getValidPolygon()
    inputToCoadd = _pixelTransform(record.getWcs(), coaddWcs)
    try:
        if polygon is not None:
            inputBounds = inputBounds.intersection(Polygon.from_legacy(polygon))
        # Only vertices are transformed, so the edges of the input stay
        # straight.
        footprint = (
            inputBounds.transform(inputToCoadd)
            .intersection(scienceFootprint)
            .transform(coaddToScienceTransform)
        )
        # The convex hull keeps the footprint in one piece when it is clipped.
        footprint = footprint.convex_hull().intersection(scienceBox)
    except NoOverlapError:
        return None
    return Polygon.from_box(footprint) if isinstance(footprint, Box) else footprint
 
 
def _warpBorder(kernelName):
    """Return the width of the border of an image that warping leaves out.
 
    A warped pixel is computed only where the whole kernel lies on the image.
 
    Parameters
    ----------
    kernelName : `str`
        Name of the warping kernel, as `lsst.afw.math.Warper` takes it.
 
    Returns
    -------
    border : `float`
        Width of the border on each side, in pixels from the edge of the
        image's bounding box; 2.5 for ``lanczos3``.
    """
    kernel = afwMath.WarpingControl(kernelName).getWarpingKernel()
    return kernel.getCtr().getX() + 0.5
 
 
@functools.cache
def _warpKernelSumSq(kernelName, nPhase=100):
    """Compute the phase-averaged sum of squared values of a warping kernel.
 
    Warping with a kernel normalized to unit sum multiplies the variance by
    this sum, which is less than one. The power the warp took off the
    diagonal became covariance between neighboring pixels, which a variance
    plane cannot represent, so dividing by it restores the noise power: the
    variance of a flux measured over a region and not a single pixel.
 
    Parameters
    ----------
    kernelName : `str`
        Name of the warping kernel, as `lsst.afw.math.Warper` takes it.
    nPhase : `int`, optional
        Number of sub-pixel phases to average over.
 
    Returns
    -------
    sumSq : `float`
        Sum of the squared 2-D kernel values, averaged over sub-pixel
        phase. Converges with large ``nPhase`` to 0.798 for ``lanczos3``,
        4/9 for ``bilinear``, and 1 for no interpolation.
 
    Notes
    -----
    The kernels are separable and the phases of the two axes independent,
    so this is the square of the 1-D average. The warper normalizes the
    kernel to unit sum at each phase, which the analytic Lanczos integral
    does not include, so the average is computed from the kernel itself.
    """
    # Disable the cache since we are evaluating fewer phase points than the
    # default cache size.
    kernel = afwMath.Warper(kernelName, cacheSize=0).getWarpingKernel()
    image = afwImage.ImageD(kernel.getDimensions())
    total = 0.0
    for phase in (np.arange(nPhase) + 0.5)/nPhase:
        kernel.setKernelParameters((float(phase), 0.0))
        kernel.computeImage(image, doNormalize=True)
        total += float((image.array.sum(axis=0)**2).sum())
    return (total/nPhase)**2
 
 
class GetTemplateConnections(
    pipeBase.PipelineTaskConnections,
    dimensions=("instrument", "visit", "detector"),
    defaultTemplates={"coaddName": "goodSeeing", "warpTypeSuffix": "", "fakesType": ""},
):
    bbox = pipeBase.connectionTypes.Input(
        doc="Bounding box of exposure to determine the geometry of the output template.",
        name="{fakesType}calexp.bbox",
        storageClass="Box2I",
        dimensions=("instrument", "visit", "detector"),
        deprecated="Replaced by the `detector` connection. Will be removed after v31.",
    )
    detector = pipeBase.connectionTypes.Input(
        doc="Detector of the exposure that we will construct the template for."
        " Its bounding box sets the geometry of the output template.",
        name="{fakesType}calexp.detector",
        storageClass="Detector",
        dimensions=("instrument", "visit", "detector"),
    )
    obs_info = pipeBase.connectionTypes.Input(
        doc="Observation metadata of the exposure the template is built for."
        " Only read with output_image_type='future', which needs this to set"
        " the obs_info component of a lsst.images.VisitImage.",
        name="{fakesType}calexp.obs_info",
        storageClass="ObservationInfo",
        dimensions=("instrument", "visit", "detector"),
    )
    wcs = pipeBase.connectionTypes.Input(
        doc="WCS of the exposure that we will construct the template for.",
        name="{fakesType}calexp.wcs",
        storageClass="Wcs",
        dimensions=("instrument", "visit", "detector"),
    )
    skyMap = pipeBase.connectionTypes.Input(
        doc="Geometry of the tracts and patches that the coadds are defined on.",
        name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
        dimensions=("skymap",),
        storageClass="SkyMap",
    )
    coaddExposures = pipeBase.connectionTypes.Input(
        doc="Coadds that may overlap the desired region, as possible inputs to the template."
        " Will be restricted to those that directly overlap the projected bounding box.",
        dimensions=("tract", "patch", "skymap", "band"),
        storageClass="ExposureF",
        name="{fakesType}{coaddName}Coadd{warpTypeSuffix}",
        multiple=True,
        deferLoad=True,
        deferGraphConstraint=True,
    )
 
    template = pipeBase.connectionTypes.Output(
        doc="Warped template, pixel matched to the bounding box and WCS.",
        dimensions=("instrument", "visit", "detector"),
        storageClass="ExposureF",
        name="{fakesType}{coaddName}Diff_templateExp{warpTypeSuffix}",
    )
 
    def __init__(self, *, config=None):
        super().__init__(config=config)
        # Kept only so that older config files still load.
        del self.bbox
        if config.requireCoaddAtGraphBuild:
            self.coaddExposures = dataclasses.replace(
                self.coaddExposures,
                deferGraphConstraint=False,
            )
        if config.output_image_type == "future":
            self.template = dataclasses.replace(self.template, storageClass="DifferenceImage")
        else:
            del self.obs_info
 
 
class GetTemplateConfig(
    pipeBase.PipelineTaskConfig, pipelineConnections=GetTemplateConnections
):
    templateBorderSize = pexConfig.Field(
        dtype=int,
        default=20,
        doc="Number of pixels to grow the requested template image to account for warping",
    )
    warp = pexConfig.ConfigField(
        dtype=afwMath.Warper.ConfigClass,
        doc="warper configuration",
    )
    doCorrectVarianceWarpCorrelation = pexConfig.Field(
        dtype=bool,
        default=True,
        doc="Correct the template variance for the correlation that warping"
        " introduced between neighboring pixels? This covers both warps: the"
        " one that built the coadds, named by coaddWarpKernel, and the one"
        " this task performs, named by warp.warpingKernelName. It is a single"
        " scalar factor, fixed entirely by those two kernels and independent"
        " of the pixel scales involved.",
    )
    coaddWarpKernel = pexConfig.ChoiceField(
        dtype=str,
        default="lanczos3",
        allowed={
            "bilinear": "bilinear interpolation",
            "lanczos3": "Lanczos kernel of order 3",
            "lanczos4": "Lanczos kernel of order 4",
            "lanczos5": "Lanczos kernel of order 5",
        },
        doc="Warping kernel that was used to build the input coadds. This"
        " cannot be recovered from the coadds themselves. It is used to"
        " correct the variance if doCorrectVarianceWarpCorrelation is True,"
        " and to trim the coadd inputs bboxes to account for the size of the"
        " warping kernel.",
    )
    coaddPsf = pexConfig.ConfigField(
        doc="Configuration for CoaddPsf",
        dtype=CoaddPsfConfig,
    )
    fwhmExposureGrid = pexConfig.Field(
        doc="Grid size to compute the average PSF FWHM in an exposure (pixels)",
        dtype=int,
        default=10,
    )
    fwhmExposureBuffer = pexConfig.Field(
        doc="Fractional buffer margin to be left out of all sides of the image during construction"
            " of grid to compute average PSF FWHM in an exposure",
        dtype=float,
        default=0.05,
    )
    varianceBackground = pexConfig.ConfigurableField(
        target=SubtractBackgroundTask,
        doc="Task to estimate the background variance.",
    )
    highVarianceThreshold = pexConfig.RangeField(
        dtype=float,
        default=4,
        min=1,
        doc="Set the HIGH_VARIANCE mask plane for regions with variance"
        " greater than the median by this factor.",
    )
    highVarianceMaskFraction = pexConfig.Field(
        dtype=float,
        default=0.1,
        doc="Minimum fraction of unmasked pixels needed to set the"
        " HIGH_VARIANCE mask plane.",
    )
    doScaleVariance = pexConfig.Field(
        dtype=bool,
        default=False,
        doc="Scale the variance of the template image to match its per-pixel"
        " scatter? Warping correlates the noise, so when combined with"
        " doCorrectVarianceWarpCorrelation this overcounts the correction.",
    )
    scaleVariance = pexConfig.ConfigurableField(
        target=ScaleVarianceTask,
        doc="Subtask to rescale the variance of the template to the statistically expected level."
    )
    requireCoaddAtGraphBuild = pexConfig.Field(
        dtype=bool,
        default=False,
        doc="If True, include the coadd dataset existence in the"
        " initial butler query during QuantumGraph generation.",
    )
    output_image_type = pexConfig.ChoiceField[str](
        "Which image type to use for the output template.",
        allowed={
            "legacy": "Write as a lsst.afw.image.ExposureF.",
            "future": "Write as a lsst.images.DifferenceImage.",
        },
        optional=False,
        default="legacy",
    )
    raiseOnUndefinedMaskMap = pexConfig.Field(
        dtype=bool,
        default=False,
        doc="Raise if the template has pixels set in a mask plane that"
        " `lsst.images` does not map? If False, log a warning and drop that"
        " plane from the output instead."
        " Only used with output_image_type='future'.",
    )
 
    def setDefaults(self):
        # Use a smaller cache: per SeparableKernel.computeCache, this should
        # give a warping error of a fraction of a count (these must match).
        self.warp.cacheSize = 100000
        self.coaddPsf.cacheSize = self.warp.cacheSize
        # The WCS for LSST should be smoothly varying, so we can use a longer
        # interpolation length for WCS evaluations.
        self.warp.interpLength = 100
        self.warp.warpingKernelName = "lanczos3"
        self.coaddPsf.warpingKernelName = self.warp.warpingKernelName
 
        # Background subtraction of the variance plane
        self.varianceBackground.algorithm = "LINEAR"
        self.varianceBackground.binSize = 32
        self.varianceBackground.useApprox = False
        self.varianceBackground.statisticsProperty = "MEDIAN"
        self.varianceBackground.doFilterSuperPixels = True
        self.varianceBackground.ignoredPixelMask = ["BAD",
                                                    "EDGE",
                                                    "DETECTED",
                                                    "DETECTED_NEGATIVE",
                                                    "NO_DATA",
                                                    ]
 
 
class GetTemplateTask(pipeBase.PipelineTask):
    ConfigClass = GetTemplateConfig
    _DefaultName = "getTemplate"
 
    def __init__(self, *args, **kwargs):
        super().__init__(*args, **kwargs)
        if self.config.doScaleVariance:
            self.makeSubtask("scaleVariance")
        self.warper = afwMath.Warper.fromConfig(self.config.warp)
        self.schema = afwTable.ExposureTable.makeMinimalSchema()
        self.schema.addField(
            "tract", type=np.int32, doc="Which tract this exposure came from."
        )
        self.schema.addField(
            "patch",
            type=np.int32,
            doc="Which patch in the tract this exposure came from.",
        )
        self.schema.addField(
            "weight",
            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):
        inputs = butlerQC.get(inputRefs)
        detector = inputs.pop("detector")
        bbox = detector.getBBox()
        wcs = inputs.pop("wcs")
        obsInfo = inputs.pop("obs_info", None)
        coaddExposures = inputs.pop("coaddExposures")
        skymap = inputs.pop("skyMap")
 
        # This should not happen with a properly configured execution context.
        assert not inputs, "runQuantum got more inputs than expected"
 
        results = self.getExposures(coaddExposures, bbox, skymap, wcs)
        physical_filter = butlerQC.quantum.dataId["physical_filter"]
        outputs = self.run(
            coaddExposureHandles=results.coaddExposures,
            bbox=bbox,
            wcs=wcs,
            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
            exposureRecord = record_from_obs_info(obsInfo, dataId["instrument"], dataId["visit"],
                                                  dataId.universe)
            outputs.template = self.convert_outputs_to_future(outputs.template,
                                                              inputRefs.coaddExposures,
                                                              detector=detector,
                                                              exposureRecord=exposureRecord,
                                                              templateInputs=outputs.templateInputs)
        butlerQC.put(outputs, outputRefs)
 
    def convert_outputs_to_future(self, template, coadd_refs, detector, exposureRecord,
                                  templateInputs=None):
        """Convert a template image to DifferenceImage format.
 
        This replaces ``template`` with an `lsst.images.DifferenceImage`
        instance and records the coadds contributed to it, along with the
        instrument they were made from and their input image identifiers.
        The Gaussian PSF set by `run` is converted to an
        `lsst.images.psfs.GaussianPointSpreadFunction`.
 
        Parameters
        ----------
        template : `lsst.afw.image.Exposure`
            Template to convert. Its detector is replaced with ``detector``
            and its mask is modified in place.
        coadd_refs : `list` [`lsst.daf.butler.DatasetRef`]
            References to the coadds that may have gone into the template.
            May be a superset of the coadds that were used; supplies the
            dataset id and RUN collection recorded for each coadd that was.
        detector : `lsst.afw.cameraGeom.Detector`
            Detector of the science image that the template was built for.
        exposureRecord : `lsst.daf.butler.DimensionRecord`
            The ``exposure`` record of the observation the template was built
            for, which supplies the observation metadata recorded on the
            output.
        templateInputs : `dict` [`tuple` [`int`, `int`], \
                `list` [`lsst.images.DifferenceImageTemplateInput`]], optional
            Input image records of each coadd, keyed by tract and patch, as
            returned by `run`. Coadds not listed record no inputs.
 
        Returns
        -------
        convertedTemplate : `lsst.images.DifferenceImage`
            The converted image.
 
        Raises
        ------
        RuntimeError
            Raised if the template has pixels set in a mask plane that
            `lsst.images.get_legacy_template_mask_planes` does not map and
            ``config.raiseOnUndefinedMaskMap`` is True, or if ``coadd_refs``
            spans more than one skymap.
        """
        # The template PSF is a Gaussian, so rebuild the CoaddPsf from the
        # template inputs to record which coadds contributed.
        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:
            raise RuntimeError("Expected the template inputs to come from one skymap;"
                               f" got {sorted(skymaps)}.")
        butler_info = {(ref.dataId["tract"], ref.dataId["patch"]): (ref.id, ref.run)
                       for ref in coadd_refs}
        planeMap = get_legacy_template_mask_planes()
        if not self.config.raiseOnUndefinedMaskMap:
            self._clearUnmappedMaskPlanes(template.mask, planeMap)
        convertedTemplate = DifferenceImage.from_legacy(
            template,
            exposure_record=exposureRecord,
            unit=astropy.units.nJy,
            plane_map=planeMap,
        )
        del template
        convertedTemplate.templates = DifferenceImageTemplateInfo.from_legacy_psf(
            convertedTemplate.sky_projection.pixel_frame,
            coaddPsf,
            skymaps.pop(),
            butler_info,
            log=self.log,
            instrument=instrument,
        )
        if templateInputs:
            for info in convertedTemplate.templates:
                info.inputs = templateInputs.get((info.tract, info.patch), [])
        return convertedTemplate
 
    def _clearUnmappedMaskPlanes(self, mask, planeMap):
        """Clear the pixels of every mask plane that conversion would reject.
 
        Parameters
        ----------
        mask : `lsst.afw.image.Mask`
            Template mask to modify in place.
        planeMap : `dict` [`str`, `lsst.images.MaskPlane`]
            Plane map the template will be converted with. The optional
            planes of `lsst.images.get_legacy_optional_mask_planes` are kept
            as well, since conversion adds them when they have pixels set.
        """
        keep = planeMap.keys() | get_legacy_optional_mask_planes().keys()
        for name in mask.getMaskPlaneDict():
            if name in keep:
                continue
            bitmask = mask.getPlaneBitMask(name)
            nSet = np.count_nonzero(mask.array & bitmask)
            if nSet:
                self.log.warning("Dropping mask plane %s from the template: %d pixels are set, but"
                                 " the template plane map does not include it.", name, nSet)
                mask.array &= ~bitmask
 
    def getExposures(self, coaddExposureHandles, bbox, skymap, wcs):
        """Return a data structure containing the coadds that overlap the
        specified bbox projected onto the sky, and a corresponding data
        structure of their dataIds.
        These are the appropriate inputs to this task's `run` method.
 
        The spatial index in the butler registry has generous padding and often
        supplies patches near, but not directly overlapping the desired region.
        This method filters the inputs so that `run` does not have to read in
        all possibly-matching coadd exposures.
 
        Parameters
        ----------
        coaddExposureHandles : `iterable` \
                          [`lsst.daf.butler.DeferredDatasetHandle` of \
                           `lsst.afw.image.Exposure`]
            Dataset handles to exposures that might overlap the desired
            region.
        bbox : `lsst.geom.Box2I`
            Template bounding box of the pixel geometry onto which the
            coaddExposures will be resampled.
        skymap : `lsst.skymap.SkyMap`
            Geometry of the tracts and patches the coadds are defined on.
        wcs : `lsst.afw.geom.SkyWcs`
            Template WCS onto which the coadds will be resampled.
 
        Returns
        -------
        result : `lsst.pipe.base.Struct`
           A struct with attributes:
 
           ``coaddExposures``
               Dict of coadd exposures that overlap the projected bbox,
               indexed on tract id
               (`dict` [`int`, `list` [`lsst.daf.butler.DeferredDatasetHandle` of
                                       `lsst.afw.image.Exposure`] ]).
           ``dataIds``
               Dict of data IDs of the coadd exposures that overlap the
               projected bbox, indexed on tract id
               (`dict` [`int`, `list [`lsst.daf.butler.DataCoordinate`] ]).
 
        Raises
        ------
        NoWorkFound
            Raised if no patches overlap the input detector bbox, or the input
            WCS is None.
        """
        if wcs is None:
            raise pipeBase.NoWorkFound(
                "WCS is None; cannot find overlapping exposures."
            )
 
        # Exposure's validPolygon would be more accurate
        detectorPolygon = geom.Box2D(bbox)
        detectorCorners = wcs.pixelToSky(detectorPolygon.getCorners())
        overlappingArea = 0
        coaddExposures = collections.defaultdict(list)
        dataIds = collections.defaultdict(list)
 
        for coaddRef in coaddExposureHandles:
            dataId = coaddRef.dataId
            patchWcs = skymap[dataId["tract"]].getWcs()
            patchBBox = skymap[dataId["tract"]][dataId["patch"]].getOuterBBox()
            patchPolygon = afwGeom.Polygon(geom.Box2D(patchBBox))
            # Calculate detector/patch overlap in patch coordinates rather than
            # detector coordinates because the skymap's inverse mapping
            # (patchWcs.skyToPixel()) is more stable than the detector's for
            # arbitrary sky coordinates.
            detectorInPatchCoordinates = afwGeom.Polygon(patchWcs.skyToPixel(detectorCorners))
            if patchPolygon.intersection(detectorInPatchCoordinates):
                overlappingArea += patchPolygon.intersectionSingle(
                    detectorInPatchCoordinates
                ).calculateArea()
                self.log.info(
                    "Using template input tract=%s, patch=%s",
                    dataId["tract"],
                    dataId["patch"],
                )
                coaddExposures[dataId["tract"]].append(coaddRef)
                dataIds[dataId["tract"]].append(dataId)
 
        if not overlappingArea:
            raise pipeBase.NoWorkFound("No patches overlap detector")
 
        return pipeBase.Struct(coaddExposures=coaddExposures, dataIds=dataIds)
 
    @timeMethod
    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.
 
        Tract and patch overlap regions are combined by a variance-weighted
        average, and the variance planes are combined with the same weights,
        not added in quadrature; the overlap regions are not statistically
        independent, because they're derived from the same original data.
        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 \
                          [`lsst.daf.butler.DeferredDatasetHandle` of \
                           `lsst.afw.image.Exposure`]]
            Coadds to be mosaicked, indexed on tract id.
        bbox : `lsst.geom.Box2I`
            Input Bounding box of the detector geometry onto which to
            resample the ``coaddExposureHandles``. The output template bbox
            will be grown by ``templateBorderSize``.
        wcs : `lsst.afw.geom.SkyWcs`
            Template WCS onto which to resample the ``coaddExposureHandles``.
        dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]]
            Record of the tract and patch of each coaddExposure, indexed on
            tract id.
        physical_filter : `str`
            Physical filter of the science image.
        visit : `int`, optional
            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
        -------
        result : `lsst.pipe.base.Struct`
           A struct with attributes:
 
           ``template``
               A template coadd exposure assembled out of patches, with
               pixels in nJy. Its PSF is a Gaussian approximation to the
               CoaddPsf of its inputs (`lsst.afw.image.ExposureF`).
           ``templateInputs``
               Input image records of each coadd that cover part of the
               science image, keyed by tract and patch (`dict`
               [`tuple` [`int`, `int`],
               `list` [`lsst.images.DifferenceImageTemplateInput`]]). Empty
               unless ``config.output_image_type`` is ``"future"``.
 
        Raises
        ------
        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)
 
        scienceBBox = bbox
        bbox = geom.Box2I(bbox)
        bbox.grow(self.config.templateBorderSize)
 
        warped = {}
        catalogs = []
        inputs = []
        plateScaleFactor = 1.0 if instrument == templateInstrument else None
 
        for tract in coaddExposureHandles:
            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[0], 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.
            unwarped, count, included = self._merge(
                maskedImages, warpedBox, catalog[0].wcs
            )
            # Delete `maskedImages` after combining into one large image to reduce peak memory use
            del maskedImages
            if count == 0:
                self.log.info(
                    "No valid pixels from coadd patches in tract %s; not including in output.",
                    tract,
                )
                continue
            warpedBox.clip(totalBox)
            potentialInput = self.warper.warpExposure(
                wcs, unwarped.subset(warpedBox), destBBox=bbox
            )
            # Delete the single large `unwarped` image after warping to reduce peak memory use
            del unwarped
            if np.all(
                potentialInput.mask.array
                & potentialInput.mask.getPlaneBitMask("NO_DATA")
            ):
                self.log.info(
                    "No overlap from coadd patches in tract %s; not including in output.",
                    tract,
                )
                continue
 
            # Trim the exposure catalog to just the patches that were used.
            tempCatalog = afwTable.ExposureCatalog(self.schema)
            tempCatalog.reserve(len(included))
            for i in included:
                tempCatalog.append(catalog[i])
                inputs.append(coaddInputs[i])
            catalogs.append(tempCatalog)
            warped[tract] = potentialInput.maskedImage
 
        if len(warped) == 0:
            raise pipeBase.NoWorkFound("No patches found to overlap science exposure.")
 
        # At this point, all entries will be valid, so we can ignore included.
        template, count, _ = self._merge(warped, bbox, wcs)
        if count == 0:
            raise pipeBase.NoWorkFound("No valid pixels in warped template.")
 
        if self.config.doScaleVariance:
            # Scale the variance of the template image before subtraction, if
            # needed. Note that the science variance is scaled
            # independently in ``AlardLuptonSubtractTask``.
            varianceFactor = self.scaleVariance.run(template.maskedImage)
            self.log.info("Template variance scaling factor: %.2f", varianceFactor)
            self.metadata["scaleTemplateVarianceFactor"] = varianceFactor
 
        self._correctVariance(template, plateScaleFactor)
 
        # Make a single catalog containing all the inputs that were accepted.
        catalog = afwTable.ExposureCatalog(self.schema)
        catalog.reserve(sum([len(c) for c in catalogs]))
        for c in catalogs:
            catalog.extend(c)
        # Only the future output type records the inputs.
        if self.config.output_image_type == "future":
            templateInputs = self._collectTemplateInputs(catalog, inputs, wcs, scienceBBox)
        else:
            templateInputs = {}
 
        # Set a mask plane for any regions with exceptionally high variance.
        self.checkHighVariance(template)
        if visit is not None:
            template.getInfo().setVisitInfo(VisitInfo(id=visit))
        template.setFilter(afwImage.FilterLabel(band, physical_filter))
        identity = afwImage.PhotoCalib(1.0)
        if photoCalib is not None and photoCalib != identity:
            template.maskedImage = photoCalib.calibrateImage(template.maskedImage)
            photoCalib = identity
        template.setPhotoCalib(photoCalib)
        template.setPsf(self._makePsf(template, catalog, wcs))
        template.setPsf(self._makeGaussianPsf(template))
 
        # Record the input coadd patches as the template's coadd inputs.
        coaddInputs = afwImage.CoaddInputs(afwTable.ExposureTable.makeMinimalSchema(), self.schema)
        coaddInputs.ccds.extend(catalog, deep=True)
        template.getInfo().setCoaddInputs(coaddInputs)
        return pipeBase.Struct(template=template, templateInputs=templateInputs)
 
    def _correctVariance(self, template, plateScaleFactor):
        """Correct the template variance plane for the effects of warping.
 
        Warping introduces a covariance between neighboring pixels, while the
        input variance plane captures only the diagonal term. Compute a
        multiplicative correction factor to apply to the final variance so that
        the detection significance on the convolved difference image is correct
        on average.
 
        Separately, also apply a scalar 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, since the scale change
        otherwise cancels between the two warps. If the plate scale of the
        science instrument is smaller than the plate scale of the coadd
        instrument, then the template pixels will be correlated and the true
        variance will be higher than the image pixel noise level would suggest.
 
        Parameters
        ----------
        template : `lsst.afw.image.Exposure`
            Assembled template; its variance plane is modified in place.
        plateScaleFactor : `float` or `None`
            Correction for the change in pixel area, from `_plateScaleFactor`.
 
        Raises
        ------
        RuntimeError
            Raised if ``plateScaleFactor`` is `None`.
        """
        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."
            )
        self.metadata["variancePlateScaleFactor"] = plateScaleFactor
        correlationFactor = 1.0
        if self.config.doCorrectVarianceWarpCorrelation:
            # Both warps correlate the noise: the one that built the coadds,
            # and the one this task just performed.
            correlationFactor = 1/(_warpKernelSumSq(self.config.coaddWarpKernel)
                                   * _warpKernelSumSq(self.config.warp.warpingKernelName))
            self.metadata["varianceWarpCorrelationFactor"] = correlationFactor
        scale = correlationFactor*plateScaleFactor
        if scale == 1.0:
            return
        template.variance.array *= scale
        self.metadata["templateVarianceCorrectionFactor"] = scale
        self.log.info(
            "Variance correction factor due to plate scale changes: %.4f",
            plateScaleFactor,
        )
        self.log.info(
            "Variance correction factor due to warping correlations: %.4f",
            correlationFactor,
        )
        self.log.info(
            "Overall template variance correction factor: %.4f",
            scale,
        )
 
    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
        ----------
        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`
            Bounding box of the template, used only to choose where to
            evaluate the science image pixel scale.
 
        Returns
        -------
        factor : `float` or `None`
            The correction factor equivalent to the ratio of the area of a pixel
            from the coadd's instrument to the science instrument, or `None` if
            this patch carries no input record from which it could be computed.
        """
        if coaddInputs is None:
            return None
        scienceScale = wcs.getPixelScale(geom.Box2D(bbox).getCenter()).asArcseconds()
 
        for record in coaddInputs.ccds:
            # Iterate through the records, and use the first one with a WCS
            # that gives a usable pixel area.
            originalWcs = record.getWcs()
            if originalWcs is None:
                continue
            center = geom.Box2D(record.getBBox()).getCenter()
            pviScale = originalWcs.getPixelScale(center).asArcseconds()
            # The area of one original pixel, expressed in science pixels.
            factor = (pviScale/scienceScale)**2
            if np.isfinite(factor) and factor > 0:
                return factor
        return None
 
    def _collectTemplateInputs(self, catalog, coaddInputs, wcs, bbox):
        """Find the input images and epochs of each coadd in the template.
 
        Parameters
        ----------
        catalog : `lsst.afw.table.ExposureCatalog`
            Coadds the template was built from, with their tract, patch and
            WCS.
        coaddInputs : `list` [`lsst.afw.image.CoaddInputs` or `None`]
            Input records of each coadd in ``catalog``, in the same order.
        wcs : `lsst.afw.geom.SkyWcs`
            WCS of the science image.
        bbox : `lsst.geom.Box2I`
            Bounding box of the science image.
 
        Returns
        -------
        templateInputs : `dict` [`tuple` [`int`, `int`], \
                `list` [`lsst.images.DifferenceImageTemplateInput`]]
            Inputs of each coadd that cover part of the science image, keyed
            by tract and patch, with the part their warps cover in science
            image pixels.
        """
        scienceBox = Box.from_legacy(bbox)
        border = _warpBorder(self.config.coaddWarpKernel)
        regions = {}
        templateInputs = {}
        for coadd, coaddInput in zip(catalog, coaddInputs, strict=True):
            coaddList = templateInputs.setdefault((coadd["tract"], coadd["patch"]), [])
            if coaddInput is None:
                continue
            scienceToCoadd = _pixelTransform(wcs, coadd.getWcs())
            scienceFootprint = scienceBox.transform(scienceToCoadd).dilated_by(_FOOTPRINT_MARGIN)
            coaddToScience = scienceToCoadd.inverted()
            for record in coaddInput.ccds:
                visitInfo = record.getVisitInfo()
                weight = record["weight"]
                if visitInfo is None or not visitInfo.getDate().isValid() or not weight > 0:
                    continue
                key = (record["visit"], record["ccd"])
                # Record a single region for a ccd in overlapping patches.
                if key not in regions:
                    regions[key] = _getImageOverlap(record, coadd.getWcs(), scienceFootprint, coaddToScience,
                                                    scienceBox, border)
                region = regions[key]
                if region is None:
                    continue
                mjd = visitInfo.getDate().get(DateTime.MJD, DateTime.TAI)
                coaddList.append(DifferenceImageTemplateInput(visit=key[0], detector=key[1],
                                                              mid_mjd_tai=mjd, weight=weight,
                                                              bounds=region))
        return templateInputs
 
    def checkHighVariance(self, template):
        """Set a mask plane for regions with unusually high variance.
 
        Parameters
        ----------
        template : `lsst.afw.image.Exposure`
            The warped template exposure, which will be modified in place.
        """
        highVarianceMaskPlaneBit = template.mask.addMaskPlane("HIGH_VARIANCE")
        ignoredPixelBits = template.mask.getPlaneBitMask(self.varianceBackground.config.ignoredPixelMask)
        goodMask = (template.mask.array & ignoredPixelBits) == 0
        goodFraction = np.count_nonzero(goodMask)/template.mask.array.size
        if goodFraction < self.config.highVarianceMaskFraction:
            self.log.info("Not setting HIGH_VARIANCE mask plane, only %2.1f%% of"
                          " pixels were unmasked for background estimation, but"
                          " %2.1f%% are required", 100*goodFraction, 100*self.config.highVarianceMaskFraction)
        else:
            varianceExposure = template.clone()
            varianceExposure.image.array = varianceExposure.variance.array
            varianceBackground = self.varianceBackground.run(varianceExposure).background.getImage().array
            threshold = self.config.highVarianceThreshold*np.nanmedian(varianceBackground)
            highVariancePix = varianceBackground > threshold
            template.mask.array[highVariancePix] |= 2**highVarianceMaskPlaneBit
 
    @staticmethod
    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, and choose
        the instrument of the template.
 
        Parameters
        ----------
        dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]]
            Record of the tract and patch of each coaddExposure.
        coaddExposures : `dict` [`int`,  `list` of \
                          [`lsst.daf.butler.DeferredDatasetHandle` of \
                           `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
        ------
        RuntimeError
            Raised if the bands or calibrations of the input exposures are not
            all the same.
        """
        bands = set(dataId["band"] for tract in dataIds for dataId in dataIds[tract])
        if len(bands) > 1:
            raise RuntimeError(f"GetTemplateTask called with multiple bands: {bands}")
        band = bands.pop()
        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]
        # 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
 
    @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
        ----------
        exposureRefs : `list` of [`lsst.daf.butler.DeferredDatasetHandle` of \
                        `lsst.afw.image.Exposure`]
            Exposures to include in the catalog.
        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
        -------
        images : `dict` [`lsst.afw.image.MaskedImage`]
            MaskedImages of each of the input exposures, for warping.
        catalog : `lsst.afw.table.ExposureCatalog`
            Catalog of metadata for each exposure
        totalBox : `lsst.geom.Box2I`
            The union of the bounding boxes of all the input exposures.
        coaddInputs : `list` [`lsst.afw.image.CoaddInputs` or `None`]
            Input records of each exposure in ``catalog``, in the same order.
        """
        catalog = afwTable.ExposureCatalog(self.schema)
        catalog.reserve(len(exposureRefs))
        exposures = (exposureRef.get() for exposureRef in exposureRefs)
        images = {}
        totalBox = geom.Box2I()
        keptCoaddInputs = []
 
        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
            keptCoaddInputs.append(coaddInputs)
            images[dataId] = coadd.maskedImage
            bbox = coadd.getBBox()
            totalBox = totalBox.expandedTo(bbox)
            record = catalog.addNew()
            record.setPsf(coadd.psf)
            record.setWcs(coadd.wcs)
            record.setPhotoCalib(coadd.photoCalib)
            record.setBBox(bbox)
            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, keptCoaddInputs
 
    def _merge(self, maskedImages, bbox, wcs):
        """Merge the images that came from one tract into one larger image,
        ignoring NaN pixels and non-finite variance pixels from individual
        exposures.
 
        Parameters
        ----------
        maskedImages : `dict` [`lsst.afw.image.MaskedImage` or
                               `lsst.afw.image.Exposure`]
            Images to be merged into one larger bounding box.
        bbox : `lsst.geom.Box2I`
            Bounding box defining the image to merge into.
        wcs : `lsst.afw.geom.SkyWcs`
            WCS of all of the input images to set on the output image.
 
        Returns
        -------
        merged : `lsst.afw.image.MaskedImage`
            Merged image with all of the inputs at their respective bbox
            positions.
        count : `int`
            Count of the number of good pixels (those with positive weights)
            in the merged image.
        included : `list` [`int`]
            List of indexes of patches that were included in the merged
            result, to be used to trim the exposure catalog.
        """
        merged = afwImage.ExposureF(bbox, wcs)
        weights = afwImage.ImageF(bbox)
        included = []  # which patches were included in the result
        for i, (dataId, maskedImage) in enumerate(maskedImages.items()):
            # Only merge into the trimmed box, to save memory
            clippedBox = geom.Box2I(maskedImage.getBBox())
            clippedBox.clip(bbox)
            if clippedBox.area == 0:
                self.log.debug("%s does not overlap template region.", dataId)
                continue  # nothing in this image overlaps the output
            maskedImage = maskedImage.subset(clippedBox)
            # Catch both zero-value and NaN variance plane pixels
            good = (maskedImage.variance.array > 0) & (
                np.isfinite(maskedImage.variance.array)
            )
            weight = maskedImage.variance.array[good] ** (-0.5)
            bad = np.isnan(maskedImage.image.array) | ~good
            # Note that modifying the patch MaskedImage in place is fine;
            # we're throwing it away at the end anyway.
            maskedImage.image.array[bad] = 0.0
            maskedImage.variance.array[bad] = 0.0
            # Reset mask, too, since these pixels don't contribute to sum.
            maskedImage.mask.array[bad] = 0
            # Cannot use `merged.maskedImage *= weight` because that operator
            # multiplies the variance by the weight twice; in this case
            # `weight` are the exact values we want to scale by.
            maskedImage.image.array[good] *= weight
            maskedImage.variance.array[good] *= weight
            weights[clippedBox].array[good] += weight
            # Free memory before creating new large arrays
            del weight
            merged.maskedImage[clippedBox] += maskedImage
            included.append(i)
 
        good = weights.array > 0
 
        # Cannot use `merged.maskedImage /= weights` because that
        # operator divides the variance by the weight twice; in this case
        # `weights` are the exact values we want to scale by.
        weights = weights.array[good]
        merged.image.array[good] /= weights
        merged.variance.array[good] /= weights
 
        merged.mask.array[~good] |= merged.mask.getPlaneBitMask("NO_DATA")
 
        return merged, good.sum(), included
 
    def _makePsf(self, template, catalog, wcs):
        """Return a PSF containing the PSF at each of the input regions.
 
        Note that although this includes all the exposures from the catalog,
        the PSF knows which part of the template the inputs came from, so when
        evaluated at a given position it will not include inputs that never
        went in to those pixels.
 
        Parameters
        ----------
        template : `lsst.afw.image.Exposure`
            Generated template the PSF is for.
        catalog : `lsst.afw.table.ExposureCatalog`
            Catalog of exposures that went into the template that contains all
            of the input PSFs.
        wcs : `lsst.afw.geom.SkyWcs`
            WCS of the template, to warp the PSFs to.
 
        Returns
        -------
        coaddPsf : `lsst.meas.algorithms.CoaddPsf`
            The meta-psf constructed from all of the input catalogs.
        """
        # CoaddPsf centroid not only must overlap image, but must overlap the
        # part of image with data. Use centroid of region with data.
        boolmask = template.mask.array & template.mask.getPlaneBitMask("NO_DATA") == 0
        maskx = afwImage.makeMaskFromArray(boolmask.astype(afwImage.MaskPixel))
        centerCoord = afwGeom.SpanSet.fromMask(maskx, 1).computeCentroid()
 
        ctrl = self.config.coaddPsf.makeControl()
        coaddPsf = CoaddPsf(
            catalog, wcs, centerCoord, ctrl.warpingKernelName, ctrl.cacheSize
        )
        return coaddPsf
 
    def _makeGaussianPsf(self, template):
        """Return a Gaussian PSF with the width of the template CoaddPsf.
 
        Parameters
        ----------
        template : `lsst.afw.image.Exposure`
            Template with the CoaddPsf to approximate. The width is measured
            at the average position of the PSF, or averaged over a grid of
            positions if the PSF cannot be evaluated there.
 
        Returns
        -------
        psf : `lsst.afw.detection.GaussianPsf`
            Gaussian with the width of the template PSF.
 
        Raises
        ------
        ValueError
            Raised if the PSF cannot be evaluated at any grid position.
        """
        # The CoaddPsf is not defined where the template has no inputs, which
        # can include its average position.
        try:
            fwhm = getPsfFwhm(template.psf)
        except lsst.pex.exceptions.Exception:
            # Catch a broad range of exceptions, since some are C++ only
            # Catching:
            #  - lsst::geom::SingularTransformException
            #  - lsst.pex.exceptions.InvalidParameterError
            #  - lsst.pex.exceptions.RangeError
            self.log.info("Unable to evaluate PSF at the average position. "
                          "Evaluating PSF on a grid of points.")
            fwhm = evaluateMeanPsfFwhm(template,
                                       fwhmExposureBuffer=self.config.fwhmExposureBuffer,
                                       fwhmExposureGrid=self.config.fwhmExposureGrid)
        sigma = fwhm/gaussian_sigma_to_fwhm
        size = 2*int(np.ceil(_PSF_STAMP_SIGMAS*sigma)) + 1
        self.metadata["templatePsfSigma"] = float(sigma)
        self.log.info("Approximated the template PSF with a Gaussian of sigma %.3f pixels.", sigma)
        return afwDetection.GaussianPsf(size, size, sigma)
 
 
class GetDcrTemplateConnections(
    GetTemplateConnections,
    dimensions=("instrument", "visit", "detector"),
    defaultTemplates={"coaddName": "dcr", "warpTypeSuffix": "", "fakesType": ""},
):
    visitInfo = pipeBase.connectionTypes.Input(
        doc="VisitInfo of calexp used to determine observing conditions.",
        name="{fakesType}calexp.visitInfo",
        storageClass="VisitInfo",
        dimensions=("instrument", "visit", "detector"),
    )
    dcrCoadds = pipeBase.connectionTypes.Input(
        doc="Input DCR template to match and subtract from the exposure",
        name="{fakesType}dcrCoadd{warpTypeSuffix}",
        storageClass="ExposureF",
        dimensions=("tract", "patch", "skymap", "band", "subfilter"),
        multiple=True,
        deferLoad=True,
    )
 
    def __init__(self, *, config=None):
        super().__init__(config=config)
        self.inputs.remove("coaddExposures")
 
 
class GetDcrTemplateConfig(
    GetTemplateConfig, pipelineConnections=GetDcrTemplateConnections
):
    """Config for GetDcrTemplateTask.
 
    Only ``output_image_type="legacy"`` is supported: the task always writes
    an `lsst.afw.image.ExposureF`.
    """
 
    numSubfilters = pexConfig.Field(
        doc="Number of subfilters in the DcrCoadd.",
        dtype=int,
        default=3,
    )
    effectiveWavelength = pexConfig.Field(
        doc="Effective wavelength of the filter in nm.",
        optional=False,
        dtype=float,
    )
    bandwidth = pexConfig.Field(
        doc="Bandwidth of the physical filter.",
        optional=False,
        dtype=float,
    )
 
    def validate(self):
        if self.effectiveWavelength is None or self.bandwidth is None:
            raise ValueError(
                "The effective wavelength and bandwidth of the physical filter "
                "must be set in the getTemplate config for DCR coadds. "
                "Required until transmission curves are used in DM-13668."
            )
 
 
class GetDcrTemplateTask(GetTemplateTask):
    ConfigClass = GetDcrTemplateConfig
    _DefaultName = "getDcrTemplate"
 
    def runQuantum(self, butlerQC, inputRefs, outputRefs):
        inputs = butlerQC.get(inputRefs)
        detector = inputs.pop("detector")
        bbox = detector.getBBox()
        wcs = inputs.pop("wcs")
        dcrCoaddExposureHandles = inputs.pop("dcrCoadds")
        skymap = inputs.pop("skyMap")
        visitInfo = inputs.pop("visitInfo")
        # This task never converts its output, so the observation metadata
        # that the conversion needs is dropped here.
        inputs.pop("obs_info", None)
 
        # This should not happen with a properly configured execution context.
        assert not inputs, "runQuantum got more inputs than expected"
 
        results = self.getExposures(
            dcrCoaddExposureHandles, bbox, skymap, wcs, visitInfo
        )
        physical_filter = butlerQC.quantum.dataId["physical_filter"]
        outputs = self.run(
            coaddExposureHandles=results.coaddExposures,
            bbox=bbox,
            wcs=wcs,
            dataIds=results.dataIds,
            physical_filter=physical_filter,
            visit=outputRefs.template.dataId["visit"],
            instrument=butlerQC.quantum.dataId["instrument"],
        )
        butlerQC.put(outputs, outputRefs)
 
    def getExposures(self, dcrCoaddExposureHandles, bbox, skymap, wcs, visitInfo):
        """Return lists of coadds and their corresponding dataIds that overlap
        the detector.
 
        The spatial index in the registry has generous padding and often
        supplies patches near, but not directly overlapping the detector.
        Filters inputs so that we don't have to read in all input coadds.
 
        Parameters
        ----------
        dcrCoaddExposureHandles :  `list` \
                                  [`lsst.daf.butler.DeferredDatasetHandle` of \
                                  `lsst.afw.image.Exposure`]
            Data references to exposures that might overlap the detector.
        bbox : `lsst.geom.Box2I`
            Template Bounding box of the detector geometry onto which to
            resample the coaddExposures.
        skymap : `lsst.skymap.SkyMap`
            Input definition of geometry/bbox and projection/wcs for
            template exposures.
        wcs : `lsst.afw.geom.SkyWcs`
            Template WCS onto which to resample the coaddExposures.
        visitInfo : `lsst.afw.image.VisitInfo`
            Metadata for the science image.
 
        Returns
        -------
        result : `lsst.pipe.base.Struct`
           A struct with attibutes:
 
           ``coaddExposures``
               Dict of coadd exposures that overlap the projected bbox,
               indexed on tract id
               (`dict` [`int`, `list` [`lsst.afw.image.Exposure`] ]).
           ``dataIds``
               Dict of data IDs of the coadd exposures that overlap the
               projected bbox, indexed on tract id
               (`dict` [`int`, `list [`lsst.daf.butler.DataCoordinate`] ]).
 
        Raises
        ------
        pipeBase.NoWorkFound
            Raised if no patches overlatp the input detector bbox.
        """
        # Check that the patches actually overlap the detector
        # Exposure's validPolygon would be more accurate
        if wcs is None:
            raise pipeBase.NoWorkFound("Exposure has no WCS; cannot create a template.")
 
        detectorPolygon = geom.Box2D(bbox)
        overlappingArea = 0
        dataIds = collections.defaultdict(list)
        patchList = dict()
        for coaddRef in dcrCoaddExposureHandles:
            dataId = coaddRef.dataId
            subfilter = dataId["subfilter"]
            patchWcs = skymap[dataId["tract"]].getWcs()
            patchBBox = skymap[dataId["tract"]][dataId["patch"]].getOuterBBox()
            patchCorners = patchWcs.pixelToSky(geom.Box2D(patchBBox).getCorners())
            patchPolygon = afwGeom.Polygon(wcs.skyToPixel(patchCorners))
            if patchPolygon.intersection(detectorPolygon):
                overlappingArea += patchPolygon.intersectionSingle(
                    detectorPolygon
                ).calculateArea()
                self.log.info(
                    "Using template input tract=%s, patch=%s, subfilter=%s"
                    % (dataId["tract"], dataId["patch"], dataId["subfilter"])
                )
                if dataId["tract"] in patchList:
                    patchList[dataId["tract"]].append(dataId["patch"])
                else:
                    patchList[dataId["tract"]] = [
                        dataId["patch"],
                    ]
                if subfilter == 0:
                    dataIds[dataId["tract"]].append(dataId)
 
        if not overlappingArea:
            raise pipeBase.NoWorkFound("No patches overlap detector")
 
        self.checkPatchList(patchList)
 
        coaddExposures = self.getDcrModel(patchList, dcrCoaddExposureHandles, visitInfo)
        return pipeBase.Struct(coaddExposures=coaddExposures, dataIds=dataIds)
 
    def checkPatchList(self, patchList):
        """Check that all of the DcrModel subfilters are present for each
        patch.
 
        Parameters
        ----------
        patchList : `dict`
            Dict of the patches containing valid data for each tract.
 
        Raises
        ------
        RuntimeError
            If the number of exposures found for a patch does not match the
            number of subfilters.
        """
        for tract in patchList:
            for patch in set(patchList[tract]):
                if patchList[tract].count(patch) != self.config.numSubfilters:
                    raise RuntimeError(
                        "Invalid number of DcrModel subfilters found: %d vs %d expected",
                        patchList[tract].count(patch),
                        self.config.numSubfilters,
                    )
 
    def getDcrModel(self, patchList, coaddRefs, visitInfo):
        """Build DCR-matched coadds from a list of exposure references.
 
        Parameters
        ----------
        patchList : `dict`
            Dict of the patches containing valid data for each tract.
        coaddRefs : `list` [`lsst.daf.butler.DeferredDatasetHandle`]
            Data references to `~lsst.afw.image.Exposure` representing
            DcrModels that overlap the detector.
        visitInfo : `lsst.afw.image.VisitInfo`
            Metadata for the science image.
 
        Returns
        -------
        coaddExposures : `list` [`lsst.afw.image.Exposure`]
            Coadd exposures that overlap the detector.
        """
        coaddExposures = collections.defaultdict(list)
        for tract in patchList:
            for patch in set(patchList[tract]):
                coaddRefList = [
                    coaddRef
                    for coaddRef in coaddRefs
                    if _selectDataRef(coaddRef, tract, patch)
                ]
 
                dcrModel = DcrModel.fromQuantum(
                    coaddRefList,
                    self.config.effectiveWavelength,
                    self.config.bandwidth,
                    self.config.numSubfilters,
                )
                coaddExposures[tract].append(dcrModel.buildMatchedExposureHandle(visitInfo=visitInfo))
        return coaddExposures
 
 
def _selectDataRef(coaddRef, tract, patch):
    condition = (coaddRef.dataId["tract"] == tract) & (
        coaddRef.dataId["patch"] == patch
    )
    return condition