Coverage for python / lsst / ip / diffim / getTemplate.py: 18%
276 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-12-12 01:36 -0800
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 overlappingArea = 0
247 coaddExposures = collections.defaultdict(list)
248 dataIds = collections.defaultdict(list)
250 for coaddRef in coaddExposureHandles:
251 dataId = coaddRef.dataId
252 patchWcs = skymap[dataId["tract"]].getWcs()
253 patchBBox = skymap[dataId["tract"]][dataId["patch"]].getOuterBBox()
254 patchCorners = patchWcs.pixelToSky(geom.Box2D(patchBBox).getCorners())
255 patchPolygon = afwGeom.Polygon(wcs.skyToPixel(patchCorners))
256 if patchPolygon.intersection(detectorPolygon):
257 overlappingArea += patchPolygon.intersectionSingle(
258 detectorPolygon
259 ).calculateArea()
260 self.log.info(
261 "Using template input tract=%s, patch=%s",
262 dataId["tract"],
263 dataId["patch"],
264 )
265 coaddExposures[dataId["tract"]].append(coaddRef)
266 dataIds[dataId["tract"]].append(dataId)
268 if not overlappingArea:
269 raise pipeBase.NoWorkFound("No patches overlap detector")
271 return pipeBase.Struct(coaddExposures=coaddExposures, dataIds=dataIds)
273 @timeMethod
274 def run(self, *, coaddExposureHandles, bbox, wcs, dataIds, physical_filter):
275 """Warp coadds from multiple tracts and patches to form a template to
276 subtract from a science image.
278 Tract and patch overlap regions are combined by a variance-weighted
279 average, and the variance planes are combined with the same weights,
280 not added in quadrature; the overlap regions are not statistically
281 independent, because they're derived from the same original data.
282 The PSF on the template is created by combining the CoaddPsf on each
283 template image into a meta-CoaddPsf.
285 Parameters
286 ----------
287 coaddExposureHandles : `dict` [`int`, `list` of \
288 [`lsst.daf.butler.DeferredDatasetHandle` of \
289 `lsst.afw.image.Exposure`]]
290 Coadds to be mosaicked, indexed on tract id.
291 bbox : `lsst.geom.Box2I`
292 Template Bounding box of the detector geometry onto which to
293 resample the ``coaddExposureHandles``. Modified in-place to include the
294 template border.
295 wcs : `lsst.afw.geom.SkyWcs`
296 Template WCS onto which to resample the ``coaddExposureHandles``.
297 dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]]
298 Record of the tract and patch of each coaddExposure, indexed on
299 tract id.
300 physical_filter : `str`
301 Physical filter of the science image.
303 Returns
304 -------
305 result : `lsst.pipe.base.Struct`
306 A struct with attributes:
308 ``template``
309 A template coadd exposure assembled out of patches
310 (`lsst.afw.image.ExposureF`).
312 Raises
313 ------
314 NoWorkFound
315 If no coadds are found with sufficient un-masked pixels.
316 """
317 band, photoCalib = self._checkInputs(dataIds, coaddExposureHandles)
319 bbox.grow(self.config.templateBorderSize)
321 warped = {}
322 catalogs = []
323 for tract in coaddExposureHandles:
324 maskedImages, catalog, totalBox = self._makeExposureCatalog(
325 coaddExposureHandles[tract], dataIds[tract]
326 )
327 warpedBox = computeWarpedBBox(catalog[0].wcs, bbox, wcs)
328 warpedBox.grow(5) # to ensure we catch all relevant input pixels
329 # Combine images from individual patches together.
330 unwarped, count, included = self._merge(
331 maskedImages, warpedBox, catalog[0].wcs
332 )
333 # Delete `maskedImages` after combining into one large image to reduce peak memory use
334 del maskedImages
335 if count == 0:
336 self.log.info(
337 "No valid pixels from coadd patches in tract %s; not including in output.",
338 tract,
339 )
340 continue
341 warpedBox.clip(totalBox)
342 potentialInput = self.warper.warpExposure(
343 wcs, unwarped.subset(warpedBox), destBBox=bbox
344 )
346 # Delete the single large `unwarped` image after warping to reduce peak memory use
347 del unwarped
348 if np.all(
349 potentialInput.mask.array
350 & potentialInput.mask.getPlaneBitMask("NO_DATA")
351 ):
352 self.log.info(
353 "No overlap from coadd patches in tract %s; not including in output.",
354 tract,
355 )
356 continue
358 # Trim the exposure catalog to just the patches that were used.
359 tempCatalog = afwTable.ExposureCatalog(self.schema)
360 tempCatalog.reserve(len(included))
361 for i in included:
362 tempCatalog.append(catalog[i])
363 catalogs.append(tempCatalog)
364 warped[tract] = potentialInput.maskedImage
366 if len(warped) == 0:
367 raise pipeBase.NoWorkFound("No patches found to overlap science exposure.")
368 # At this point, all entries will be valid, so we can ignore included.
369 template, count, _ = self._merge(warped, bbox, wcs)
370 if count == 0:
371 raise pipeBase.NoWorkFound("No valid pixels in warped template.")
373 # Make a single catalog containing all the inputs that were accepted.
374 catalog = afwTable.ExposureCatalog(self.schema)
375 catalog.reserve(sum([len(c) for c in catalogs]))
376 for c in catalogs:
377 catalog.extend(c)
379 # Set a mask plane for any regions with exceptionally high variance.
380 self.checkHighVariance(template)
381 template.setPsf(self._makePsf(template, catalog, wcs))
382 template.setFilter(afwImage.FilterLabel(band, physical_filter))
383 template.setPhotoCalib(photoCalib)
384 return pipeBase.Struct(template=template)
386 def checkHighVariance(self, template):
387 """Set a mask plane for regions with unusually high variance.
389 Parameters
390 ----------
391 template : `lsst.afw.image.Exposure`
392 The warped template exposure, which will be modified in place.
393 """
394 highVarianceMaskPlaneBit = template.mask.addMaskPlane("HIGH_VARIANCE")
395 ignoredPixelBits = template.mask.getPlaneBitMask(self.varianceBackground.config.ignoredPixelMask)
396 goodMask = (template.mask.array & ignoredPixelBits) == 0
397 goodFraction = np.count_nonzero(goodMask)/template.mask.array.size
398 if goodFraction < self.config.highVarianceMaskFraction:
399 self.log.info("Not setting HIGH_VARIANCE mask plane, only %2.1f%% of"
400 " pixels were unmasked for background estimation, but"
401 " %2.1f%% are required", 100*goodFraction, 100*self.config.highVarianceMaskFraction)
402 else:
403 varianceExposure = template.clone()
404 varianceExposure.image.array = varianceExposure.variance.array
405 varianceBackground = self.varianceBackground.run(varianceExposure).background.getImage().array
406 threshold = self.config.highVarianceThreshold*np.nanmedian(varianceBackground)
407 highVariancePix = varianceBackground > threshold
408 template.mask.array[highVariancePix] |= 2**highVarianceMaskPlaneBit
410 @staticmethod
411 def _checkInputs(dataIds, coaddExposures):
412 """Check that the all the dataIds are from the same band and that
413 the exposures all have the same photometric calibration.
415 Parameters
416 ----------
417 dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]]
418 Record of the tract and patch of each coaddExposure.
419 coaddExposures : `dict` [`int`, `list` of \
420 [`lsst.daf.butler.DeferredDatasetHandle` of \
421 `lsst.afw.image.Exposure` or
422 `lsst.afw.image.Exposure`]]
423 Coadds to be mosaicked.
425 Returns
426 -------
427 band : `str`
428 Filter band of all the input exposures.
429 photoCalib : `lsst.afw.image.PhotoCalib`
430 Photometric calibration of all of the input exposures.
432 Raises
433 ------
434 RuntimeError
435 Raised if the bands or calibrations of the input exposures are not
436 all the same.
437 """
438 bands = set(dataId["band"] for tract in dataIds for dataId in dataIds[tract])
439 if len(bands) > 1:
440 raise RuntimeError(f"GetTemplateTask called with multiple bands: {bands}")
441 band = bands.pop()
442 photoCalibs = [
443 exposure.get(component="photoCalib")
444 for exposures in coaddExposures.values()
445 for exposure in exposures
446 ]
447 if not all([photoCalibs[0] == x for x in photoCalibs]):
448 msg = f"GetTemplateTask called with exposures with different photoCalibs: {photoCalibs}"
449 raise RuntimeError(msg)
450 photoCalib = photoCalibs[0]
451 return band, photoCalib
453 def _makeExposureCatalog(self, exposureRefs, dataIds):
454 """Make an exposure catalog for one tract.
456 Parameters
457 ----------
458 exposureRefs : `list` of [`lsst.daf.butler.DeferredDatasetHandle` of \
459 `lsst.afw.image.Exposure`]
460 Exposures to include in the catalog.
461 dataIds : `list` [`lsst.daf.butler.DataCoordinate`]
462 Data ids of each of the included exposures; must have "tract" and
463 "patch" entries.
465 Returns
466 -------
467 images : `dict` [`lsst.afw.image.MaskedImage`]
468 MaskedImages of each of the input exposures, for warping.
469 catalog : `lsst.afw.table.ExposureCatalog`
470 Catalog of metadata for each exposure
471 totalBox : `lsst.geom.Box2I`
472 The union of the bounding boxes of all the input exposures.
473 """
474 catalog = afwTable.ExposureCatalog(self.schema)
475 catalog.reserve(len(exposureRefs))
476 exposures = (exposureRef.get() for exposureRef in exposureRefs)
477 images = {}
478 totalBox = geom.Box2I()
480 for coadd, dataId in zip(exposures, dataIds):
481 images[dataId] = coadd.maskedImage
482 bbox = coadd.getBBox()
483 totalBox = totalBox.expandedTo(bbox)
484 record = catalog.addNew()
485 record.setPsf(coadd.psf)
486 record.setWcs(coadd.wcs)
487 record.setPhotoCalib(coadd.photoCalib)
488 record.setBBox(bbox)
489 record.setValidPolygon(afwGeom.Polygon(geom.Box2D(bbox).getCorners()))
490 record.set("tract", dataId["tract"])
491 record.set("patch", dataId["patch"])
492 # Weight is used by CoaddPsf, but the PSFs from overlapping patches
493 # should be very similar, so this value mostly shouldn't matter.
494 record.set("weight", 1)
496 return images, catalog, totalBox
498 def _merge(self, maskedImages, bbox, wcs):
499 """Merge the images that came from one tract into one larger image,
500 ignoring NaN pixels and non-finite variance pixels from individual
501 exposures.
503 Parameters
504 ----------
505 maskedImages : `dict` [`lsst.afw.image.MaskedImage` or
506 `lsst.afw.image.Exposure`]
507 Images to be merged into one larger bounding box.
508 bbox : `lsst.geom.Box2I`
509 Bounding box defining the image to merge into.
510 wcs : `lsst.afw.geom.SkyWcs`
511 WCS of all of the input images to set on the output image.
513 Returns
514 -------
515 merged : `lsst.afw.image.MaskedImage`
516 Merged image with all of the inputs at their respective bbox
517 positions.
518 count : `int`
519 Count of the number of good pixels (those with positive weights)
520 in the merged image.
521 included : `list` [`int`]
522 List of indexes of patches that were included in the merged
523 result, to be used to trim the exposure catalog.
524 """
525 merged = afwImage.ExposureF(bbox, wcs)
526 weights = afwImage.ImageF(bbox)
527 included = [] # which patches were included in the result
528 for i, (dataId, maskedImage) in enumerate(maskedImages.items()):
529 # Only merge into the trimmed box, to save memory
530 clippedBox = geom.Box2I(maskedImage.getBBox())
531 clippedBox.clip(bbox)
532 if clippedBox.area == 0:
533 self.log.debug("%s does not overlap template region.", dataId)
534 continue # nothing in this image overlaps the output
535 maskedImage = maskedImage.subset(clippedBox)
536 # Catch both zero-value and NaN variance plane pixels
537 good = (maskedImage.variance.array > 0) & (
538 np.isfinite(maskedImage.variance.array)
539 )
540 weight = maskedImage.variance.array[good] ** (-0.5)
541 bad = np.isnan(maskedImage.image.array) | ~good
542 # Note that modifying the patch MaskedImage in place is fine;
543 # we're throwing it away at the end anyway.
544 maskedImage.image.array[bad] = 0.0
545 maskedImage.variance.array[bad] = 0.0
546 # Reset mask, too, since these pixels don't contribute to sum.
547 maskedImage.mask.array[bad] = 0
548 # Cannot use `merged.maskedImage *= weight` because that operator
549 # multiplies the variance by the weight twice; in this case
550 # `weight` are the exact values we want to scale by.
551 maskedImage.image.array[good] *= weight
552 maskedImage.variance.array[good] *= weight
553 weights[clippedBox].array[good] += weight
554 # Free memory before creating new large arrays
555 del weight
556 merged.maskedImage[clippedBox] += maskedImage
557 included.append(i)
559 good = weights.array > 0
561 # Cannot use `merged.maskedImage /= weights` because that
562 # operator divides the variance by the weight twice; in this case
563 # `weights` are the exact values we want to scale by.
564 weights = weights.array[good]
565 merged.image.array[good] /= weights
566 merged.variance.array[good] /= weights
568 merged.mask.array[~good] |= merged.mask.getPlaneBitMask("NO_DATA")
570 return merged, good.sum(), included
572 def _makePsf(self, template, catalog, wcs):
573 """Return a PSF containing the PSF at each of the input regions.
575 Note that although this includes all the exposures from the catalog,
576 the PSF knows which part of the template the inputs came from, so when
577 evaluated at a given position it will not include inputs that never
578 went in to those pixels.
580 Parameters
581 ----------
582 template : `lsst.afw.image.Exposure`
583 Generated template the PSF is for.
584 catalog : `lsst.afw.table.ExposureCatalog`
585 Catalog of exposures that went into the template that contains all
586 of the input PSFs.
587 wcs : `lsst.afw.geom.SkyWcs`
588 WCS of the template, to warp the PSFs to.
590 Returns
591 -------
592 coaddPsf : `lsst.meas.algorithms.CoaddPsf`
593 The meta-psf constructed from all of the input catalogs.
594 """
595 # CoaddPsf centroid not only must overlap image, but must overlap the
596 # part of image with data. Use centroid of region with data.
597 boolmask = template.mask.array & template.mask.getPlaneBitMask("NO_DATA") == 0
598 maskx = afwImage.makeMaskFromArray(boolmask.astype(afwImage.MaskPixel))
599 centerCoord = afwGeom.SpanSet.fromMask(maskx, 1).computeCentroid()
601 ctrl = self.config.coaddPsf.makeControl()
602 coaddPsf = CoaddPsf(
603 catalog, wcs, centerCoord, ctrl.warpingKernelName, ctrl.cacheSize
604 )
605 return coaddPsf
608class GetDcrTemplateConnections(
609 GetTemplateConnections,
610 dimensions=("instrument", "visit", "detector"),
611 defaultTemplates={"coaddName": "dcr", "warpTypeSuffix": "", "fakesType": ""},
612):
613 visitInfo = pipeBase.connectionTypes.Input(
614 doc="VisitInfo of calexp used to determine observing conditions.",
615 name="{fakesType}calexp.visitInfo",
616 storageClass="VisitInfo",
617 dimensions=("instrument", "visit", "detector"),
618 )
619 dcrCoadds = pipeBase.connectionTypes.Input(
620 doc="Input DCR template to match and subtract from the exposure",
621 name="{fakesType}dcrCoadd{warpTypeSuffix}",
622 storageClass="ExposureF",
623 dimensions=("tract", "patch", "skymap", "band", "subfilter"),
624 multiple=True,
625 deferLoad=True,
626 )
628 def __init__(self, *, config=None):
629 super().__init__(config=config)
630 self.inputs.remove("coaddExposures")
633class GetDcrTemplateConfig(
634 GetTemplateConfig, pipelineConnections=GetDcrTemplateConnections
635):
636 numSubfilters = pexConfig.Field(
637 doc="Number of subfilters in the DcrCoadd.",
638 dtype=int,
639 default=3,
640 )
641 effectiveWavelength = pexConfig.Field(
642 doc="Effective wavelength of the filter in nm.",
643 optional=False,
644 dtype=float,
645 )
646 bandwidth = pexConfig.Field(
647 doc="Bandwidth of the physical filter.",
648 optional=False,
649 dtype=float,
650 )
652 def validate(self):
653 if self.effectiveWavelength is None or self.bandwidth is None:
654 raise ValueError(
655 "The effective wavelength and bandwidth of the physical filter "
656 "must be set in the getTemplate config for DCR coadds. "
657 "Required until transmission curves are used in DM-13668."
658 )
661class GetDcrTemplateTask(GetTemplateTask):
662 ConfigClass = GetDcrTemplateConfig
663 _DefaultName = "getDcrTemplate"
665 def runQuantum(self, butlerQC, inputRefs, outputRefs):
666 inputs = butlerQC.get(inputRefs)
667 bbox = inputs.pop("bbox")
668 wcs = inputs.pop("wcs")
669 dcrCoaddExposureHandles = inputs.pop("dcrCoadds")
670 skymap = inputs.pop("skyMap")
671 visitInfo = inputs.pop("visitInfo")
673 # This should not happen with a properly configured execution context.
674 assert not inputs, "runQuantum got more inputs than expected"
676 results = self.getExposures(
677 dcrCoaddExposureHandles, bbox, skymap, wcs, visitInfo
678 )
679 physical_filter = butlerQC.quantum.dataId["physical_filter"]
680 outputs = self.run(
681 coaddExposureHandles=results.coaddExposures,
682 bbox=bbox,
683 wcs=wcs,
684 dataIds=results.dataIds,
685 physical_filter=physical_filter,
686 )
687 butlerQC.put(outputs, outputRefs)
689 def getExposures(self, dcrCoaddExposureHandles, bbox, skymap, wcs, visitInfo):
690 """Return lists of coadds and their corresponding dataIds that overlap
691 the detector.
693 The spatial index in the registry has generous padding and often
694 supplies patches near, but not directly overlapping the detector.
695 Filters inputs so that we don't have to read in all input coadds.
697 Parameters
698 ----------
699 dcrCoaddExposureHandles : `list` \
700 [`lsst.daf.butler.DeferredDatasetHandle` of \
701 `lsst.afw.image.Exposure`]
702 Data references to exposures that might overlap the detector.
703 bbox : `lsst.geom.Box2I`
704 Template Bounding box of the detector geometry onto which to
705 resample the coaddExposures.
706 skymap : `lsst.skymap.SkyMap`
707 Input definition of geometry/bbox and projection/wcs for
708 template exposures.
709 wcs : `lsst.afw.geom.SkyWcs`
710 Template WCS onto which to resample the coaddExposures.
711 visitInfo : `lsst.afw.image.VisitInfo`
712 Metadata for the science image.
714 Returns
715 -------
716 result : `lsst.pipe.base.Struct`
717 A struct with attibutes:
719 ``coaddExposures``
720 Dict of coadd exposures that overlap the projected bbox,
721 indexed on tract id
722 (`dict` [`int`, `list` [`lsst.afw.image.Exposure`] ]).
723 ``dataIds``
724 Dict of data IDs of the coadd exposures that overlap the
725 projected bbox, indexed on tract id
726 (`dict` [`int`, `list [`lsst.daf.butler.DataCoordinate`] ]).
728 Raises
729 ------
730 pipeBase.NoWorkFound
731 Raised if no patches overlatp the input detector bbox.
732 """
733 # Check that the patches actually overlap the detector
734 # Exposure's validPolygon would be more accurate
735 if wcs is None:
736 raise pipeBase.NoWorkFound("Exposure has no WCS; cannot create a template.")
738 detectorPolygon = geom.Box2D(bbox)
739 overlappingArea = 0
740 dataIds = collections.defaultdict(list)
741 patchList = dict()
742 for coaddRef in dcrCoaddExposureHandles:
743 dataId = coaddRef.dataId
744 subfilter = dataId["subfilter"]
745 patchWcs = skymap[dataId["tract"]].getWcs()
746 patchBBox = skymap[dataId["tract"]][dataId["patch"]].getOuterBBox()
747 patchCorners = patchWcs.pixelToSky(geom.Box2D(patchBBox).getCorners())
748 patchPolygon = afwGeom.Polygon(wcs.skyToPixel(patchCorners))
749 if patchPolygon.intersection(detectorPolygon):
750 overlappingArea += patchPolygon.intersectionSingle(
751 detectorPolygon
752 ).calculateArea()
753 self.log.info(
754 "Using template input tract=%s, patch=%s, subfilter=%s"
755 % (dataId["tract"], dataId["patch"], dataId["subfilter"])
756 )
757 if dataId["tract"] in patchList:
758 patchList[dataId["tract"]].append(dataId["patch"])
759 else:
760 patchList[dataId["tract"]] = [
761 dataId["patch"],
762 ]
763 if subfilter == 0:
764 dataIds[dataId["tract"]].append(dataId)
766 if not overlappingArea:
767 raise pipeBase.NoWorkFound("No patches overlap detector")
769 self.checkPatchList(patchList)
771 coaddExposures = self.getDcrModel(patchList, dcrCoaddExposureHandles, visitInfo)
772 return pipeBase.Struct(coaddExposures=coaddExposures, dataIds=dataIds)
774 def checkPatchList(self, patchList):
775 """Check that all of the DcrModel subfilters are present for each
776 patch.
778 Parameters
779 ----------
780 patchList : `dict`
781 Dict of the patches containing valid data for each tract.
783 Raises
784 ------
785 RuntimeError
786 If the number of exposures found for a patch does not match the
787 number of subfilters.
788 """
789 for tract in patchList:
790 for patch in set(patchList[tract]):
791 if patchList[tract].count(patch) != self.config.numSubfilters:
792 raise RuntimeError(
793 "Invalid number of DcrModel subfilters found: %d vs %d expected",
794 patchList[tract].count(patch),
795 self.config.numSubfilters,
796 )
798 def getDcrModel(self, patchList, coaddRefs, visitInfo):
799 """Build DCR-matched coadds from a list of exposure references.
801 Parameters
802 ----------
803 patchList : `dict`
804 Dict of the patches containing valid data for each tract.
805 coaddRefs : `list` [`lsst.daf.butler.DeferredDatasetHandle`]
806 Data references to `~lsst.afw.image.Exposure` representing
807 DcrModels that overlap the detector.
808 visitInfo : `lsst.afw.image.VisitInfo`
809 Metadata for the science image.
811 Returns
812 -------
813 coaddExposures : `list` [`lsst.afw.image.Exposure`]
814 Coadd exposures that overlap the detector.
815 """
816 coaddExposures = collections.defaultdict(list)
817 for tract in patchList:
818 for patch in set(patchList[tract]):
819 coaddRefList = [
820 coaddRef
821 for coaddRef in coaddRefs
822 if _selectDataRef(coaddRef, tract, patch)
823 ]
825 dcrModel = DcrModel.fromQuantum(
826 coaddRefList,
827 self.config.effectiveWavelength,
828 self.config.bandwidth,
829 self.config.numSubfilters,
830 )
831 coaddExposures[tract].append(dcrModel.buildMatchedExposureHandle(visitInfo=visitInfo))
832 return coaddExposures
835def _selectDataRef(coaddRef, tract, patch):
836 condition = (coaddRef.dataId["tract"] == tract) & (
837 coaddRef.dataId["patch"] == patch
838 )
839 return condition