Coverage for python/lsst/ip/diffim/getTemplate.py: 63%
283 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-18 18:18 +0000
1# This file is part of ip_diffim.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21import collections
23import numpy as np
25import lsst.afw.image as afwImage
26import lsst.geom as geom
27import lsst.afw.geom as afwGeom
28from lsst.afw.image import VisitInfo
29import lsst.afw.table as afwTable
30from lsst.afw.math._warper import computeWarpedBBox
31import lsst.afw.math as afwMath
32import lsst.pex.config as pexConfig
33import lsst.pipe.base as pipeBase
35from lsst.skymap import BaseSkyMap
36from lsst.ip.diffim.dcrModel import DcrModel
37from lsst.meas.algorithms import CoaddPsf, CoaddPsfConfig, SubtractBackgroundTask
38from lsst.utils.timer import timeMethod
40__all__ = [
41 "GetTemplateTask",
42 "GetTemplateConfig",
43 "GetDcrTemplateTask",
44 "GetDcrTemplateConfig",
45]
48class GetTemplateConnections(
49 pipeBase.PipelineTaskConnections,
50 dimensions=("instrument", "visit", "detector"),
51 defaultTemplates={"coaddName": "goodSeeing", "warpTypeSuffix": "", "fakesType": ""},
52):
53 bbox = pipeBase.connectionTypes.Input(
54 doc="Bounding box of exposure to determine the geometry of the output template.",
55 name="{fakesType}calexp.bbox",
56 storageClass="Box2I",
57 dimensions=("instrument", "visit", "detector"),
58 )
59 wcs = pipeBase.connectionTypes.Input(
60 doc="WCS of the exposure that we will construct the template for.",
61 name="{fakesType}calexp.wcs",
62 storageClass="Wcs",
63 dimensions=("instrument", "visit", "detector"),
64 )
65 skyMap = pipeBase.connectionTypes.Input(
66 doc="Geometry of the tracts and patches that the coadds are defined on.",
67 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
68 dimensions=("skymap",),
69 storageClass="SkyMap",
70 )
71 coaddExposures = pipeBase.connectionTypes.Input(
72 doc="Coadds that may overlap the desired region, as possible inputs to the template."
73 " Will be restricted to those that directly overlap the projected bounding box.",
74 dimensions=("tract", "patch", "skymap", "band"),
75 storageClass="ExposureF",
76 name="{fakesType}{coaddName}Coadd{warpTypeSuffix}",
77 multiple=True,
78 deferLoad=True,
79 deferGraphConstraint=True,
80 )
82 template = pipeBase.connectionTypes.Output(
83 doc="Warped template, pixel matched to the bounding box and WCS.",
84 dimensions=("instrument", "visit", "detector"),
85 storageClass="ExposureF",
86 name="{fakesType}{coaddName}Diff_templateExp{warpTypeSuffix}",
87 )
90class GetTemplateConfig(
91 pipeBase.PipelineTaskConfig, pipelineConnections=GetTemplateConnections
92):
93 templateBorderSize = pexConfig.Field(
94 dtype=int,
95 default=20,
96 doc="Number of pixels to grow the requested template image to account for warping",
97 )
98 warp = pexConfig.ConfigField(
99 dtype=afwMath.Warper.ConfigClass,
100 doc="warper configuration",
101 )
102 coaddPsf = pexConfig.ConfigField(
103 doc="Configuration for CoaddPsf",
104 dtype=CoaddPsfConfig,
105 )
106 varianceBackground = pexConfig.ConfigurableField(
107 target=SubtractBackgroundTask,
108 doc="Task to estimate the background variance.",
109 )
110 highVarianceThreshold = pexConfig.RangeField(
111 dtype=float,
112 default=4,
113 min=1,
114 doc="Set the HIGH_VARIANCE mask plane for regions with variance"
115 " greater than the median by this factor.",
116 )
117 highVarianceMaskFraction = pexConfig.Field(
118 dtype=float,
119 default=0.1,
120 doc="Minimum fraction of unmasked pixels needed to set the"
121 " HIGH_VARIANCE mask plane.",
122 )
124 def setDefaults(self):
125 # Use a smaller cache: per SeparableKernel.computeCache, this should
126 # give a warping error of a fraction of a count (these must match).
127 self.warp.cacheSize = 100000
128 self.coaddPsf.cacheSize = self.warp.cacheSize
129 # The WCS for LSST should be smoothly varying, so we can use a longer
130 # interpolation length for WCS evaluations.
131 self.warp.interpLength = 100
132 self.warp.warpingKernelName = "lanczos3"
133 self.coaddPsf.warpingKernelName = self.warp.warpingKernelName
135 # Background subtraction of the variance plane
136 self.varianceBackground.algorithm = "LINEAR"
137 self.varianceBackground.binSize = 32
138 self.varianceBackground.useApprox = False
139 self.varianceBackground.statisticsProperty = "MEDIAN"
140 self.varianceBackground.doFilterSuperPixels = True
141 self.varianceBackground.ignoredPixelMask = ["BAD",
142 "EDGE",
143 "DETECTED",
144 "DETECTED_NEGATIVE",
145 "NO_DATA",
146 ]
149class GetTemplateTask(pipeBase.PipelineTask):
150 ConfigClass = GetTemplateConfig
151 _DefaultName = "getTemplate"
153 def __init__(self, *args, **kwargs):
154 super().__init__(*args, **kwargs)
155 self.warper = afwMath.Warper.fromConfig(self.config.warp)
156 self.schema = afwTable.ExposureTable.makeMinimalSchema()
157 self.schema.addField(
158 "tract", type=np.int32, doc="Which tract this exposure came from."
159 )
160 self.schema.addField(
161 "patch",
162 type=np.int32,
163 doc="Which patch in the tract this exposure came from.",
164 )
165 self.schema.addField(
166 "weight",
167 type=float,
168 doc="Weight for each exposure, used to make the CoaddPsf; should always be 1.",
169 )
170 self.makeSubtask("varianceBackground")
172 def runQuantum(self, butlerQC, inputRefs, outputRefs):
173 inputs = butlerQC.get(inputRefs)
174 bbox = inputs.pop("bbox")
175 wcs = inputs.pop("wcs")
176 coaddExposures = inputs.pop("coaddExposures")
177 skymap = inputs.pop("skyMap")
179 # This should not happen with a properly configured execution context.
180 assert not inputs, "runQuantum got more inputs than expected"
182 results = self.getExposures(coaddExposures, bbox, skymap, wcs)
183 physical_filter = butlerQC.quantum.dataId["physical_filter"]
184 outputs = self.run(
185 coaddExposureHandles=results.coaddExposures,
186 bbox=bbox,
187 wcs=wcs,
188 dataIds=results.dataIds,
189 physical_filter=physical_filter,
190 visit=outputRefs.template.dataId["visit"],
191 )
192 butlerQC.put(outputs, outputRefs)
194 def getExposures(self, coaddExposureHandles, bbox, skymap, wcs):
195 """Return a data structure containing the coadds that overlap the
196 specified bbox projected onto the sky, and a corresponding data
197 structure of their dataIds.
198 These are the appropriate inputs to this task's `run` method.
200 The spatial index in the butler registry has generous padding and often
201 supplies patches near, but not directly overlapping the desired region.
202 This method filters the inputs so that `run` does not have to read in
203 all possibly-matching coadd exposures.
205 Parameters
206 ----------
207 coaddExposureHandles : `iterable` \
208 [`lsst.daf.butler.DeferredDatasetHandle` of \
209 `lsst.afw.image.Exposure`]
210 Dataset handles to exposures that might overlap the desired
211 region.
212 bbox : `lsst.geom.Box2I`
213 Template bounding box of the pixel geometry onto which the
214 coaddExposures will be resampled.
215 skymap : `lsst.skymap.SkyMap`
216 Geometry of the tracts and patches the coadds are defined on.
217 wcs : `lsst.afw.geom.SkyWcs`
218 Template WCS onto which the coadds will be resampled.
220 Returns
221 -------
222 result : `lsst.pipe.base.Struct`
223 A struct with attributes:
225 ``coaddExposures``
226 Dict of coadd exposures that overlap the projected bbox,
227 indexed on tract id
228 (`dict` [`int`, `list` [`lsst.daf.butler.DeferredDatasetHandle` of
229 `lsst.afw.image.Exposure`] ]).
230 ``dataIds``
231 Dict of data IDs of the coadd exposures that overlap the
232 projected bbox, indexed on tract id
233 (`dict` [`int`, `list [`lsst.daf.butler.DataCoordinate`] ]).
235 Raises
236 ------
237 NoWorkFound
238 Raised if no patches overlap the input detector bbox, or the input
239 WCS is None.
240 """
241 if wcs is None:
242 raise pipeBase.NoWorkFound(
243 "WCS is None; cannot find overlapping exposures."
244 )
246 # Exposure's validPolygon would be more accurate
247 detectorPolygon = geom.Box2D(bbox)
248 detectorCorners = wcs.pixelToSky(detectorPolygon.getCorners())
249 overlappingArea = 0
250 coaddExposures = collections.defaultdict(list)
251 dataIds = collections.defaultdict(list)
253 for coaddRef in coaddExposureHandles:
254 dataId = coaddRef.dataId
255 patchWcs = skymap[dataId["tract"]].getWcs()
256 patchBBox = skymap[dataId["tract"]][dataId["patch"]].getOuterBBox()
257 patchPolygon = afwGeom.Polygon(geom.Box2D(patchBBox))
258 # Calculate detector/patch overlap in patch coordinates rather than
259 # detector coordinates because the skymap's inverse mapping
260 # (patchWcs.skyToPixel()) is more stable than the detector's for
261 # arbitrary sky coordinates.
262 detectorInPatchCoordinates = afwGeom.Polygon(patchWcs.skyToPixel(detectorCorners))
263 if patchPolygon.intersection(detectorInPatchCoordinates):
264 overlappingArea += patchPolygon.intersectionSingle(
265 detectorInPatchCoordinates
266 ).calculateArea()
267 self.log.info(
268 "Using template input tract=%s, patch=%s",
269 dataId["tract"],
270 dataId["patch"],
271 )
272 coaddExposures[dataId["tract"]].append(coaddRef)
273 dataIds[dataId["tract"]].append(dataId)
275 if not overlappingArea:
276 raise pipeBase.NoWorkFound("No patches overlap detector")
278 return pipeBase.Struct(coaddExposures=coaddExposures, dataIds=dataIds)
280 @timeMethod
281 def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter, visit=None):
282 """Warp coadds from multiple tracts and patches to form a template to
283 subtract from a science image.
285 Tract and patch overlap regions are combined by a variance-weighted
286 average, and the variance planes are combined with the same weights,
287 not added in quadrature; the overlap regions are not statistically
288 independent, because they're derived from the same original data.
289 The PSF on the template is created by combining the CoaddPsf on each
290 template image into a meta-CoaddPsf.
292 Parameters
293 ----------
294 coaddExposureHandles : `dict` [`int`, `list` of \
295 [`lsst.daf.butler.DeferredDatasetHandle` of \
296 `lsst.afw.image.Exposure`]]
297 Coadds to be mosaicked, indexed on tract id.
298 bbox : `lsst.geom.Box2I`
299 Template Bounding box of the detector geometry onto which to
300 resample the ``coaddExposureHandles``. Modified in-place to include the
301 template border.
302 wcs : `lsst.afw.geom.SkyWcs`
303 Template WCS onto which to resample the ``coaddExposureHandles``.
304 dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]]
305 Record of the tract and patch of each coaddExposure, indexed on
306 tract id.
307 physical_filter : `str`
308 Physical filter of the science image.
309 visit : `int`, optional
310 If supplied, over-write the visit ID in the template's visitInfo
311 so that downstream source injection tasks can link the template and
312 science image for the visit.
314 Returns
315 -------
316 result : `lsst.pipe.base.Struct`
317 A struct with attributes:
319 ``template``
320 A template coadd exposure assembled out of patches
321 (`lsst.afw.image.ExposureF`).
323 Raises
324 ------
325 NoWorkFound
326 If no coadds are found with sufficient un-masked pixels.
327 """
328 band, photoCalib = self._checkInputs(dataIds, coaddExposureHandles)
330 bbox.grow(self.config.templateBorderSize)
332 warped = {}
333 catalogs = []
334 for tract in coaddExposureHandles:
335 maskedImages, catalog, totalBox = self._makeExposureCatalog(
336 coaddExposureHandles[tract], dataIds[tract]
337 )
338 warpedBox = computeWarpedBBox(catalog[0].wcs, bbox, wcs)
339 warpedBox.grow(5) # to ensure we catch all relevant input pixels
340 # Combine images from individual patches together.
341 unwarped, count, included = self._merge(
342 maskedImages, warpedBox, catalog[0].wcs
343 )
344 # Delete `maskedImages` after combining into one large image to reduce peak memory use
345 del maskedImages
346 if count == 0:
347 self.log.info(
348 "No valid pixels from coadd patches in tract %s; not including in output.",
349 tract,
350 )
351 continue
352 warpedBox.clip(totalBox)
353 potentialInput = self.warper.warpExposure(
354 wcs, unwarped.subset(warpedBox), destBBox=bbox
355 )
357 # Delete the single large `unwarped` image after warping to reduce peak memory use
358 del unwarped
359 if np.all(
360 potentialInput.mask.array
361 & potentialInput.mask.getPlaneBitMask("NO_DATA")
362 ):
363 self.log.info(
364 "No overlap from coadd patches in tract %s; not including in output.",
365 tract,
366 )
367 continue
369 # Trim the exposure catalog to just the patches that were used.
370 tempCatalog = afwTable.ExposureCatalog(self.schema)
371 tempCatalog.reserve(len(included))
372 for i in included:
373 tempCatalog.append(catalog[i])
374 catalogs.append(tempCatalog)
375 warped[tract] = potentialInput.maskedImage
377 if len(warped) == 0:
378 raise pipeBase.NoWorkFound("No patches found to overlap science exposure.")
379 # At this point, all entries will be valid, so we can ignore included.
380 template, count, _ = self._merge(warped, bbox, wcs)
381 if count == 0: 381 ↛ 382line 381 didn't jump to line 382 because the condition on line 381 was never true
382 raise pipeBase.NoWorkFound("No valid pixels in warped template.")
384 # Make a single catalog containing all the inputs that were accepted.
385 catalog = afwTable.ExposureCatalog(self.schema)
386 catalog.reserve(sum([len(c) for c in catalogs]))
387 for c in catalogs:
388 catalog.extend(c)
390 # Set a mask plane for any regions with exceptionally high variance.
391 self.checkHighVariance(template)
392 if visit is not None: 392 ↛ 393line 392 didn't jump to line 393 because the condition on line 392 was never true
393 template.getInfo().setVisitInfo(VisitInfo(id=visit))
394 template.setFilter(afwImage.FilterLabel(band, physical_filter))
395 template.setPhotoCalib(photoCalib)
396 template.setPsf(self._makePsf(template, catalog, wcs))
397 # Record the input coadd patches as the template's coadd inputs.
398 coaddInputs = afwImage.CoaddInputs(afwTable.ExposureTable.makeMinimalSchema(), self.schema)
399 coaddInputs.ccds.extend(catalog, deep=True)
400 template.getInfo().setCoaddInputs(coaddInputs)
401 return pipeBase.Struct(template=template)
403 def checkHighVariance(self, template):
404 """Set a mask plane for regions with unusually high variance.
406 Parameters
407 ----------
408 template : `lsst.afw.image.Exposure`
409 The warped template exposure, which will be modified in place.
410 """
411 highVarianceMaskPlaneBit = template.mask.addMaskPlane("HIGH_VARIANCE")
412 ignoredPixelBits = template.mask.getPlaneBitMask(self.varianceBackground.config.ignoredPixelMask)
413 goodMask = (template.mask.array & ignoredPixelBits) == 0
414 goodFraction = np.count_nonzero(goodMask)/template.mask.array.size
415 if goodFraction < self.config.highVarianceMaskFraction: 415 ↛ 416line 415 didn't jump to line 416 because the condition on line 415 was never true
416 self.log.info("Not setting HIGH_VARIANCE mask plane, only %2.1f%% of"
417 " pixels were unmasked for background estimation, but"
418 " %2.1f%% are required", 100*goodFraction, 100*self.config.highVarianceMaskFraction)
419 else:
420 varianceExposure = template.clone()
421 varianceExposure.image.array = varianceExposure.variance.array
422 varianceBackground = self.varianceBackground.run(varianceExposure).background.getImage().array
423 threshold = self.config.highVarianceThreshold*np.nanmedian(varianceBackground)
424 highVariancePix = varianceBackground > threshold
425 template.mask.array[highVariancePix] |= 2**highVarianceMaskPlaneBit
427 @staticmethod
428 def _checkInputs(dataIds, coaddExposures):
429 """Check that the all the dataIds are from the same band and that
430 the exposures all have the same photometric calibration.
432 Parameters
433 ----------
434 dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]]
435 Record of the tract and patch of each coaddExposure.
436 coaddExposures : `dict` [`int`, `list` of \
437 [`lsst.daf.butler.DeferredDatasetHandle` of \
438 `lsst.afw.image.Exposure` or
439 `lsst.afw.image.Exposure`]]
440 Coadds to be mosaicked.
442 Returns
443 -------
444 band : `str`
445 Filter band of all the input exposures.
446 photoCalib : `lsst.afw.image.PhotoCalib`
447 Photometric calibration of all of the input exposures.
449 Raises
450 ------
451 RuntimeError
452 Raised if the bands or calibrations of the input exposures are not
453 all the same.
454 """
455 bands = set(dataId["band"] for tract in dataIds for dataId in dataIds[tract])
456 if len(bands) > 1: 456 ↛ 457line 456 didn't jump to line 457 because the condition on line 456 was never true
457 raise RuntimeError(f"GetTemplateTask called with multiple bands: {bands}")
458 band = bands.pop()
459 photoCalibs = [
460 exposure.get(component="photoCalib")
461 for exposures in coaddExposures.values()
462 for exposure in exposures
463 ]
464 if not all([photoCalibs[0] == x for x in photoCalibs]): 464 ↛ 465line 464 didn't jump to line 465 because the condition on line 464 was never true
465 msg = f"GetTemplateTask called with exposures with different photoCalibs: {photoCalibs}"
466 raise RuntimeError(msg)
467 photoCalib = photoCalibs[0]
468 return band, photoCalib
470 def _makeExposureCatalog(self, exposureRefs, dataIds):
471 """Make an exposure catalog for one tract.
473 Parameters
474 ----------
475 exposureRefs : `list` of [`lsst.daf.butler.DeferredDatasetHandle` of \
476 `lsst.afw.image.Exposure`]
477 Exposures to include in the catalog.
478 dataIds : `list` [`lsst.daf.butler.DataCoordinate`]
479 Data ids of each of the included exposures; must have "tract" and
480 "patch" entries.
482 Returns
483 -------
484 images : `dict` [`lsst.afw.image.MaskedImage`]
485 MaskedImages of each of the input exposures, for warping.
486 catalog : `lsst.afw.table.ExposureCatalog`
487 Catalog of metadata for each exposure
488 totalBox : `lsst.geom.Box2I`
489 The union of the bounding boxes of all the input exposures.
490 """
491 catalog = afwTable.ExposureCatalog(self.schema)
492 catalog.reserve(len(exposureRefs))
493 exposures = (exposureRef.get() for exposureRef in exposureRefs)
494 images = {}
495 totalBox = geom.Box2I()
497 for coadd, dataId in zip(exposures, dataIds):
498 images[dataId] = coadd.maskedImage
499 bbox = coadd.getBBox()
500 totalBox = totalBox.expandedTo(bbox)
501 record = catalog.addNew()
502 record.setPsf(coadd.psf)
503 record.setWcs(coadd.wcs)
504 record.setPhotoCalib(coadd.photoCalib)
505 record.setBBox(bbox)
506 record.setValidPolygon(afwGeom.Polygon(geom.Box2D(bbox).getCorners()))
507 record.set("tract", dataId["tract"])
508 record.set("patch", dataId["patch"])
509 # Weight is used by CoaddPsf, but the PSFs from overlapping patches
510 # should be very similar, so this value mostly shouldn't matter.
511 record.set("weight", 1)
513 return images, catalog, totalBox
515 def _merge(self, maskedImages, bbox, wcs):
516 """Merge the images that came from one tract into one larger image,
517 ignoring NaN pixels and non-finite variance pixels from individual
518 exposures.
520 Parameters
521 ----------
522 maskedImages : `dict` [`lsst.afw.image.MaskedImage` or
523 `lsst.afw.image.Exposure`]
524 Images to be merged into one larger bounding box.
525 bbox : `lsst.geom.Box2I`
526 Bounding box defining the image to merge into.
527 wcs : `lsst.afw.geom.SkyWcs`
528 WCS of all of the input images to set on the output image.
530 Returns
531 -------
532 merged : `lsst.afw.image.MaskedImage`
533 Merged image with all of the inputs at their respective bbox
534 positions.
535 count : `int`
536 Count of the number of good pixels (those with positive weights)
537 in the merged image.
538 included : `list` [`int`]
539 List of indexes of patches that were included in the merged
540 result, to be used to trim the exposure catalog.
541 """
542 merged = afwImage.ExposureF(bbox, wcs)
543 weights = afwImage.ImageF(bbox)
544 included = [] # which patches were included in the result
545 for i, (dataId, maskedImage) in enumerate(maskedImages.items()):
546 # Only merge into the trimmed box, to save memory
547 clippedBox = geom.Box2I(maskedImage.getBBox())
548 clippedBox.clip(bbox)
549 if clippedBox.area == 0:
550 self.log.debug("%s does not overlap template region.", dataId)
551 continue # nothing in this image overlaps the output
552 maskedImage = maskedImage.subset(clippedBox)
553 # Catch both zero-value and NaN variance plane pixels
554 good = (maskedImage.variance.array > 0) & (
555 np.isfinite(maskedImage.variance.array)
556 )
557 weight = maskedImage.variance.array[good] ** (-0.5)
558 bad = np.isnan(maskedImage.image.array) | ~good
559 # Note that modifying the patch MaskedImage in place is fine;
560 # we're throwing it away at the end anyway.
561 maskedImage.image.array[bad] = 0.0
562 maskedImage.variance.array[bad] = 0.0
563 # Reset mask, too, since these pixels don't contribute to sum.
564 maskedImage.mask.array[bad] = 0
565 # Cannot use `merged.maskedImage *= weight` because that operator
566 # multiplies the variance by the weight twice; in this case
567 # `weight` are the exact values we want to scale by.
568 maskedImage.image.array[good] *= weight
569 maskedImage.variance.array[good] *= weight
570 weights[clippedBox].array[good] += weight
571 # Free memory before creating new large arrays
572 del weight
573 merged.maskedImage[clippedBox] += maskedImage
574 included.append(i)
576 good = weights.array > 0
578 # Cannot use `merged.maskedImage /= weights` because that
579 # operator divides the variance by the weight twice; in this case
580 # `weights` are the exact values we want to scale by.
581 weights = weights.array[good]
582 merged.image.array[good] /= weights
583 merged.variance.array[good] /= weights
585 merged.mask.array[~good] |= merged.mask.getPlaneBitMask("NO_DATA")
587 return merged, good.sum(), included
589 def _makePsf(self, template, catalog, wcs):
590 """Return a PSF containing the PSF at each of the input regions.
592 Note that although this includes all the exposures from the catalog,
593 the PSF knows which part of the template the inputs came from, so when
594 evaluated at a given position it will not include inputs that never
595 went in to those pixels.
597 Parameters
598 ----------
599 template : `lsst.afw.image.Exposure`
600 Generated template the PSF is for.
601 catalog : `lsst.afw.table.ExposureCatalog`
602 Catalog of exposures that went into the template that contains all
603 of the input PSFs.
604 wcs : `lsst.afw.geom.SkyWcs`
605 WCS of the template, to warp the PSFs to.
607 Returns
608 -------
609 coaddPsf : `lsst.meas.algorithms.CoaddPsf`
610 The meta-psf constructed from all of the input catalogs.
611 """
612 # CoaddPsf centroid not only must overlap image, but must overlap the
613 # part of image with data. Use centroid of region with data.
614 boolmask = template.mask.array & template.mask.getPlaneBitMask("NO_DATA") == 0
615 maskx = afwImage.makeMaskFromArray(boolmask.astype(afwImage.MaskPixel))
616 centerCoord = afwGeom.SpanSet.fromMask(maskx, 1).computeCentroid()
618 ctrl = self.config.coaddPsf.makeControl()
619 coaddPsf = CoaddPsf(
620 catalog, wcs, centerCoord, ctrl.warpingKernelName, ctrl.cacheSize
621 )
622 return coaddPsf
625class GetDcrTemplateConnections(
626 GetTemplateConnections,
627 dimensions=("instrument", "visit", "detector"),
628 defaultTemplates={"coaddName": "dcr", "warpTypeSuffix": "", "fakesType": ""},
629):
630 visitInfo = pipeBase.connectionTypes.Input(
631 doc="VisitInfo of calexp used to determine observing conditions.",
632 name="{fakesType}calexp.visitInfo",
633 storageClass="VisitInfo",
634 dimensions=("instrument", "visit", "detector"),
635 )
636 dcrCoadds = pipeBase.connectionTypes.Input(
637 doc="Input DCR template to match and subtract from the exposure",
638 name="{fakesType}dcrCoadd{warpTypeSuffix}",
639 storageClass="ExposureF",
640 dimensions=("tract", "patch", "skymap", "band", "subfilter"),
641 multiple=True,
642 deferLoad=True,
643 )
645 def __init__(self, *, config=None):
646 super().__init__(config=config)
647 self.inputs.remove("coaddExposures")
650class GetDcrTemplateConfig(
651 GetTemplateConfig, pipelineConnections=GetDcrTemplateConnections
652):
653 numSubfilters = pexConfig.Field(
654 doc="Number of subfilters in the DcrCoadd.",
655 dtype=int,
656 default=3,
657 )
658 effectiveWavelength = pexConfig.Field(
659 doc="Effective wavelength of the filter in nm.",
660 optional=False,
661 dtype=float,
662 )
663 bandwidth = pexConfig.Field(
664 doc="Bandwidth of the physical filter.",
665 optional=False,
666 dtype=float,
667 )
669 def validate(self):
670 if self.effectiveWavelength is None or self.bandwidth is None:
671 raise ValueError(
672 "The effective wavelength and bandwidth of the physical filter "
673 "must be set in the getTemplate config for DCR coadds. "
674 "Required until transmission curves are used in DM-13668."
675 )
678class GetDcrTemplateTask(GetTemplateTask):
679 ConfigClass = GetDcrTemplateConfig
680 _DefaultName = "getDcrTemplate"
682 def runQuantum(self, butlerQC, inputRefs, outputRefs):
683 inputs = butlerQC.get(inputRefs)
684 bbox = inputs.pop("bbox")
685 wcs = inputs.pop("wcs")
686 dcrCoaddExposureHandles = inputs.pop("dcrCoadds")
687 skymap = inputs.pop("skyMap")
688 visitInfo = inputs.pop("visitInfo")
690 # This should not happen with a properly configured execution context.
691 assert not inputs, "runQuantum got more inputs than expected"
693 results = self.getExposures(
694 dcrCoaddExposureHandles, bbox, skymap, wcs, visitInfo
695 )
696 physical_filter = butlerQC.quantum.dataId["physical_filter"]
697 outputs = self.run(
698 coaddExposureHandles=results.coaddExposures,
699 bbox=bbox,
700 wcs=wcs,
701 dataIds=results.dataIds,
702 physical_filter=physical_filter,
703 )
704 butlerQC.put(outputs, outputRefs)
706 def getExposures(self, dcrCoaddExposureHandles, bbox, skymap, wcs, visitInfo):
707 """Return lists of coadds and their corresponding dataIds that overlap
708 the detector.
710 The spatial index in the registry has generous padding and often
711 supplies patches near, but not directly overlapping the detector.
712 Filters inputs so that we don't have to read in all input coadds.
714 Parameters
715 ----------
716 dcrCoaddExposureHandles : `list` \
717 [`lsst.daf.butler.DeferredDatasetHandle` of \
718 `lsst.afw.image.Exposure`]
719 Data references to exposures that might overlap the detector.
720 bbox : `lsst.geom.Box2I`
721 Template Bounding box of the detector geometry onto which to
722 resample the coaddExposures.
723 skymap : `lsst.skymap.SkyMap`
724 Input definition of geometry/bbox and projection/wcs for
725 template exposures.
726 wcs : `lsst.afw.geom.SkyWcs`
727 Template WCS onto which to resample the coaddExposures.
728 visitInfo : `lsst.afw.image.VisitInfo`
729 Metadata for the science image.
731 Returns
732 -------
733 result : `lsst.pipe.base.Struct`
734 A struct with attibutes:
736 ``coaddExposures``
737 Dict of coadd exposures that overlap the projected bbox,
738 indexed on tract id
739 (`dict` [`int`, `list` [`lsst.afw.image.Exposure`] ]).
740 ``dataIds``
741 Dict of data IDs of the coadd exposures that overlap the
742 projected bbox, indexed on tract id
743 (`dict` [`int`, `list [`lsst.daf.butler.DataCoordinate`] ]).
745 Raises
746 ------
747 pipeBase.NoWorkFound
748 Raised if no patches overlatp the input detector bbox.
749 """
750 # Check that the patches actually overlap the detector
751 # Exposure's validPolygon would be more accurate
752 if wcs is None:
753 raise pipeBase.NoWorkFound("Exposure has no WCS; cannot create a template.")
755 detectorPolygon = geom.Box2D(bbox)
756 overlappingArea = 0
757 dataIds = collections.defaultdict(list)
758 patchList = dict()
759 for coaddRef in dcrCoaddExposureHandles:
760 dataId = coaddRef.dataId
761 subfilter = dataId["subfilter"]
762 patchWcs = skymap[dataId["tract"]].getWcs()
763 patchBBox = skymap[dataId["tract"]][dataId["patch"]].getOuterBBox()
764 patchCorners = patchWcs.pixelToSky(geom.Box2D(patchBBox).getCorners())
765 patchPolygon = afwGeom.Polygon(wcs.skyToPixel(patchCorners))
766 if patchPolygon.intersection(detectorPolygon):
767 overlappingArea += patchPolygon.intersectionSingle(
768 detectorPolygon
769 ).calculateArea()
770 self.log.info(
771 "Using template input tract=%s, patch=%s, subfilter=%s"
772 % (dataId["tract"], dataId["patch"], dataId["subfilter"])
773 )
774 if dataId["tract"] in patchList:
775 patchList[dataId["tract"]].append(dataId["patch"])
776 else:
777 patchList[dataId["tract"]] = [
778 dataId["patch"],
779 ]
780 if subfilter == 0:
781 dataIds[dataId["tract"]].append(dataId)
783 if not overlappingArea:
784 raise pipeBase.NoWorkFound("No patches overlap detector")
786 self.checkPatchList(patchList)
788 coaddExposures = self.getDcrModel(patchList, dcrCoaddExposureHandles, visitInfo)
789 return pipeBase.Struct(coaddExposures=coaddExposures, dataIds=dataIds)
791 def checkPatchList(self, patchList):
792 """Check that all of the DcrModel subfilters are present for each
793 patch.
795 Parameters
796 ----------
797 patchList : `dict`
798 Dict of the patches containing valid data for each tract.
800 Raises
801 ------
802 RuntimeError
803 If the number of exposures found for a patch does not match the
804 number of subfilters.
805 """
806 for tract in patchList:
807 for patch in set(patchList[tract]):
808 if patchList[tract].count(patch) != self.config.numSubfilters:
809 raise RuntimeError(
810 "Invalid number of DcrModel subfilters found: %d vs %d expected",
811 patchList[tract].count(patch),
812 self.config.numSubfilters,
813 )
815 def getDcrModel(self, patchList, coaddRefs, visitInfo):
816 """Build DCR-matched coadds from a list of exposure references.
818 Parameters
819 ----------
820 patchList : `dict`
821 Dict of the patches containing valid data for each tract.
822 coaddRefs : `list` [`lsst.daf.butler.DeferredDatasetHandle`]
823 Data references to `~lsst.afw.image.Exposure` representing
824 DcrModels that overlap the detector.
825 visitInfo : `lsst.afw.image.VisitInfo`
826 Metadata for the science image.
828 Returns
829 -------
830 coaddExposures : `list` [`lsst.afw.image.Exposure`]
831 Coadd exposures that overlap the detector.
832 """
833 coaddExposures = collections.defaultdict(list)
834 for tract in patchList:
835 for patch in set(patchList[tract]):
836 coaddRefList = [
837 coaddRef
838 for coaddRef in coaddRefs
839 if _selectDataRef(coaddRef, tract, patch)
840 ]
842 dcrModel = DcrModel.fromQuantum(
843 coaddRefList,
844 self.config.effectiveWavelength,
845 self.config.bandwidth,
846 self.config.numSubfilters,
847 )
848 coaddExposures[tract].append(dcrModel.buildMatchedExposureHandle(visitInfo=visitInfo))
849 return coaddExposures
852def _selectDataRef(coaddRef, tract, patch):
853 condition = (coaddRef.dataId["tract"] == tract) & (
854 coaddRef.dataId["patch"] == patch
855 )
856 return condition