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