Coverage for python/lsst/ip/diffim/getTemplate.py: 18%
280 statements
« prev ^ index » next coverage.py v7.14.1, created at 2026-06-12 07:52 +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:
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:
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 return pipeBase.Struct(template=template)
399 def checkHighVariance(self, template):
400 """Set a mask plane for regions with unusually high variance.
402 Parameters
403 ----------
404 template : `lsst.afw.image.Exposure`
405 The warped template exposure, which will be modified in place.
406 """
407 highVarianceMaskPlaneBit = template.mask.addMaskPlane("HIGH_VARIANCE")
408 ignoredPixelBits = template.mask.getPlaneBitMask(self.varianceBackground.config.ignoredPixelMask)
409 goodMask = (template.mask.array & ignoredPixelBits) == 0
410 goodFraction = np.count_nonzero(goodMask)/template.mask.array.size
411 if goodFraction < self.config.highVarianceMaskFraction:
412 self.log.info("Not setting HIGH_VARIANCE mask plane, only %2.1f%% of"
413 " pixels were unmasked for background estimation, but"
414 " %2.1f%% are required", 100*goodFraction, 100*self.config.highVarianceMaskFraction)
415 else:
416 varianceExposure = template.clone()
417 varianceExposure.image.array = varianceExposure.variance.array
418 varianceBackground = self.varianceBackground.run(varianceExposure).background.getImage().array
419 threshold = self.config.highVarianceThreshold*np.nanmedian(varianceBackground)
420 highVariancePix = varianceBackground > threshold
421 template.mask.array[highVariancePix] |= 2**highVarianceMaskPlaneBit
423 @staticmethod
424 def _checkInputs(dataIds, coaddExposures):
425 """Check that the all the dataIds are from the same band and that
426 the exposures all have the same photometric calibration.
428 Parameters
429 ----------
430 dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]]
431 Record of the tract and patch of each coaddExposure.
432 coaddExposures : `dict` [`int`, `list` of \
433 [`lsst.daf.butler.DeferredDatasetHandle` of \
434 `lsst.afw.image.Exposure` or
435 `lsst.afw.image.Exposure`]]
436 Coadds to be mosaicked.
438 Returns
439 -------
440 band : `str`
441 Filter band of all the input exposures.
442 photoCalib : `lsst.afw.image.PhotoCalib`
443 Photometric calibration of all of the input exposures.
445 Raises
446 ------
447 RuntimeError
448 Raised if the bands or calibrations of the input exposures are not
449 all the same.
450 """
451 bands = set(dataId["band"] for tract in dataIds for dataId in dataIds[tract])
452 if len(bands) > 1:
453 raise RuntimeError(f"GetTemplateTask called with multiple bands: {bands}")
454 band = bands.pop()
455 photoCalibs = [
456 exposure.get(component="photoCalib")
457 for exposures in coaddExposures.values()
458 for exposure in exposures
459 ]
460 if not all([photoCalibs[0] == x for x in photoCalibs]):
461 msg = f"GetTemplateTask called with exposures with different photoCalibs: {photoCalibs}"
462 raise RuntimeError(msg)
463 photoCalib = photoCalibs[0]
464 return band, photoCalib
466 def _makeExposureCatalog(self, exposureRefs, dataIds):
467 """Make an exposure catalog for one tract.
469 Parameters
470 ----------
471 exposureRefs : `list` of [`lsst.daf.butler.DeferredDatasetHandle` of \
472 `lsst.afw.image.Exposure`]
473 Exposures to include in the catalog.
474 dataIds : `list` [`lsst.daf.butler.DataCoordinate`]
475 Data ids of each of the included exposures; must have "tract" and
476 "patch" entries.
478 Returns
479 -------
480 images : `dict` [`lsst.afw.image.MaskedImage`]
481 MaskedImages of each of the input exposures, for warping.
482 catalog : `lsst.afw.table.ExposureCatalog`
483 Catalog of metadata for each exposure
484 totalBox : `lsst.geom.Box2I`
485 The union of the bounding boxes of all the input exposures.
486 """
487 catalog = afwTable.ExposureCatalog(self.schema)
488 catalog.reserve(len(exposureRefs))
489 exposures = (exposureRef.get() for exposureRef in exposureRefs)
490 images = {}
491 totalBox = geom.Box2I()
493 for coadd, dataId in zip(exposures, dataIds):
494 images[dataId] = coadd.maskedImage
495 bbox = coadd.getBBox()
496 totalBox = totalBox.expandedTo(bbox)
497 record = catalog.addNew()
498 record.setPsf(coadd.psf)
499 record.setWcs(coadd.wcs)
500 record.setPhotoCalib(coadd.photoCalib)
501 record.setBBox(bbox)
502 record.setValidPolygon(afwGeom.Polygon(geom.Box2D(bbox).getCorners()))
503 record.set("tract", dataId["tract"])
504 record.set("patch", dataId["patch"])
505 # Weight is used by CoaddPsf, but the PSFs from overlapping patches
506 # should be very similar, so this value mostly shouldn't matter.
507 record.set("weight", 1)
509 return images, catalog, totalBox
511 def _merge(self, maskedImages, bbox, wcs):
512 """Merge the images that came from one tract into one larger image,
513 ignoring NaN pixels and non-finite variance pixels from individual
514 exposures.
516 Parameters
517 ----------
518 maskedImages : `dict` [`lsst.afw.image.MaskedImage` or
519 `lsst.afw.image.Exposure`]
520 Images to be merged into one larger bounding box.
521 bbox : `lsst.geom.Box2I`
522 Bounding box defining the image to merge into.
523 wcs : `lsst.afw.geom.SkyWcs`
524 WCS of all of the input images to set on the output image.
526 Returns
527 -------
528 merged : `lsst.afw.image.MaskedImage`
529 Merged image with all of the inputs at their respective bbox
530 positions.
531 count : `int`
532 Count of the number of good pixels (those with positive weights)
533 in the merged image.
534 included : `list` [`int`]
535 List of indexes of patches that were included in the merged
536 result, to be used to trim the exposure catalog.
537 """
538 merged = afwImage.ExposureF(bbox, wcs)
539 weights = afwImage.ImageF(bbox)
540 included = [] # which patches were included in the result
541 for i, (dataId, maskedImage) in enumerate(maskedImages.items()):
542 # Only merge into the trimmed box, to save memory
543 clippedBox = geom.Box2I(maskedImage.getBBox())
544 clippedBox.clip(bbox)
545 if clippedBox.area == 0:
546 self.log.debug("%s does not overlap template region.", dataId)
547 continue # nothing in this image overlaps the output
548 maskedImage = maskedImage.subset(clippedBox)
549 # Catch both zero-value and NaN variance plane pixels
550 good = (maskedImage.variance.array > 0) & (
551 np.isfinite(maskedImage.variance.array)
552 )
553 weight = maskedImage.variance.array[good] ** (-0.5)
554 bad = np.isnan(maskedImage.image.array) | ~good
555 # Note that modifying the patch MaskedImage in place is fine;
556 # we're throwing it away at the end anyway.
557 maskedImage.image.array[bad] = 0.0
558 maskedImage.variance.array[bad] = 0.0
559 # Reset mask, too, since these pixels don't contribute to sum.
560 maskedImage.mask.array[bad] = 0
561 # Cannot use `merged.maskedImage *= weight` because that operator
562 # multiplies the variance by the weight twice; in this case
563 # `weight` are the exact values we want to scale by.
564 maskedImage.image.array[good] *= weight
565 maskedImage.variance.array[good] *= weight
566 weights[clippedBox].array[good] += weight
567 # Free memory before creating new large arrays
568 del weight
569 merged.maskedImage[clippedBox] += maskedImage
570 included.append(i)
572 good = weights.array > 0
574 # Cannot use `merged.maskedImage /= weights` because that
575 # operator divides the variance by the weight twice; in this case
576 # `weights` are the exact values we want to scale by.
577 weights = weights.array[good]
578 merged.image.array[good] /= weights
579 merged.variance.array[good] /= weights
581 merged.mask.array[~good] |= merged.mask.getPlaneBitMask("NO_DATA")
583 return merged, good.sum(), included
585 def _makePsf(self, template, catalog, wcs):
586 """Return a PSF containing the PSF at each of the input regions.
588 Note that although this includes all the exposures from the catalog,
589 the PSF knows which part of the template the inputs came from, so when
590 evaluated at a given position it will not include inputs that never
591 went in to those pixels.
593 Parameters
594 ----------
595 template : `lsst.afw.image.Exposure`
596 Generated template the PSF is for.
597 catalog : `lsst.afw.table.ExposureCatalog`
598 Catalog of exposures that went into the template that contains all
599 of the input PSFs.
600 wcs : `lsst.afw.geom.SkyWcs`
601 WCS of the template, to warp the PSFs to.
603 Returns
604 -------
605 coaddPsf : `lsst.meas.algorithms.CoaddPsf`
606 The meta-psf constructed from all of the input catalogs.
607 """
608 # CoaddPsf centroid not only must overlap image, but must overlap the
609 # part of image with data. Use centroid of region with data.
610 boolmask = template.mask.array & template.mask.getPlaneBitMask("NO_DATA") == 0
611 maskx = afwImage.makeMaskFromArray(boolmask.astype(afwImage.MaskPixel))
612 centerCoord = afwGeom.SpanSet.fromMask(maskx, 1).computeCentroid()
614 ctrl = self.config.coaddPsf.makeControl()
615 coaddPsf = CoaddPsf(
616 catalog, wcs, centerCoord, ctrl.warpingKernelName, ctrl.cacheSize
617 )
618 return coaddPsf
621class GetDcrTemplateConnections(
622 GetTemplateConnections,
623 dimensions=("instrument", "visit", "detector"),
624 defaultTemplates={"coaddName": "dcr", "warpTypeSuffix": "", "fakesType": ""},
625):
626 visitInfo = pipeBase.connectionTypes.Input(
627 doc="VisitInfo of calexp used to determine observing conditions.",
628 name="{fakesType}calexp.visitInfo",
629 storageClass="VisitInfo",
630 dimensions=("instrument", "visit", "detector"),
631 )
632 dcrCoadds = pipeBase.connectionTypes.Input(
633 doc="Input DCR template to match and subtract from the exposure",
634 name="{fakesType}dcrCoadd{warpTypeSuffix}",
635 storageClass="ExposureF",
636 dimensions=("tract", "patch", "skymap", "band", "subfilter"),
637 multiple=True,
638 deferLoad=True,
639 )
641 def __init__(self, *, config=None):
642 super().__init__(config=config)
643 self.inputs.remove("coaddExposures")
646class GetDcrTemplateConfig(
647 GetTemplateConfig, pipelineConnections=GetDcrTemplateConnections
648):
649 numSubfilters = pexConfig.Field(
650 doc="Number of subfilters in the DcrCoadd.",
651 dtype=int,
652 default=3,
653 )
654 effectiveWavelength = pexConfig.Field(
655 doc="Effective wavelength of the filter in nm.",
656 optional=False,
657 dtype=float,
658 )
659 bandwidth = pexConfig.Field(
660 doc="Bandwidth of the physical filter.",
661 optional=False,
662 dtype=float,
663 )
665 def validate(self):
666 if self.effectiveWavelength is None or self.bandwidth is None:
667 raise ValueError(
668 "The effective wavelength and bandwidth of the physical filter "
669 "must be set in the getTemplate config for DCR coadds. "
670 "Required until transmission curves are used in DM-13668."
671 )
674class GetDcrTemplateTask(GetTemplateTask):
675 ConfigClass = GetDcrTemplateConfig
676 _DefaultName = "getDcrTemplate"
678 def runQuantum(self, butlerQC, inputRefs, outputRefs):
679 inputs = butlerQC.get(inputRefs)
680 bbox = inputs.pop("bbox")
681 wcs = inputs.pop("wcs")
682 dcrCoaddExposureHandles = inputs.pop("dcrCoadds")
683 skymap = inputs.pop("skyMap")
684 visitInfo = inputs.pop("visitInfo")
686 # This should not happen with a properly configured execution context.
687 assert not inputs, "runQuantum got more inputs than expected"
689 results = self.getExposures(
690 dcrCoaddExposureHandles, bbox, skymap, wcs, visitInfo
691 )
692 physical_filter = butlerQC.quantum.dataId["physical_filter"]
693 outputs = self.run(
694 coaddExposureHandles=results.coaddExposures,
695 bbox=bbox,
696 wcs=wcs,
697 dataIds=results.dataIds,
698 physical_filter=physical_filter,
699 )
700 butlerQC.put(outputs, outputRefs)
702 def getExposures(self, dcrCoaddExposureHandles, bbox, skymap, wcs, visitInfo):
703 """Return lists of coadds and their corresponding dataIds that overlap
704 the detector.
706 The spatial index in the registry has generous padding and often
707 supplies patches near, but not directly overlapping the detector.
708 Filters inputs so that we don't have to read in all input coadds.
710 Parameters
711 ----------
712 dcrCoaddExposureHandles : `list` \
713 [`lsst.daf.butler.DeferredDatasetHandle` of \
714 `lsst.afw.image.Exposure`]
715 Data references to exposures that might overlap the detector.
716 bbox : `lsst.geom.Box2I`
717 Template Bounding box of the detector geometry onto which to
718 resample the coaddExposures.
719 skymap : `lsst.skymap.SkyMap`
720 Input definition of geometry/bbox and projection/wcs for
721 template exposures.
722 wcs : `lsst.afw.geom.SkyWcs`
723 Template WCS onto which to resample the coaddExposures.
724 visitInfo : `lsst.afw.image.VisitInfo`
725 Metadata for the science image.
727 Returns
728 -------
729 result : `lsst.pipe.base.Struct`
730 A struct with attibutes:
732 ``coaddExposures``
733 Dict of coadd exposures that overlap the projected bbox,
734 indexed on tract id
735 (`dict` [`int`, `list` [`lsst.afw.image.Exposure`] ]).
736 ``dataIds``
737 Dict of data IDs of the coadd exposures that overlap the
738 projected bbox, indexed on tract id
739 (`dict` [`int`, `list [`lsst.daf.butler.DataCoordinate`] ]).
741 Raises
742 ------
743 pipeBase.NoWorkFound
744 Raised if no patches overlatp the input detector bbox.
745 """
746 # Check that the patches actually overlap the detector
747 # Exposure's validPolygon would be more accurate
748 if wcs is None:
749 raise pipeBase.NoWorkFound("Exposure has no WCS; cannot create a template.")
751 detectorPolygon = geom.Box2D(bbox)
752 overlappingArea = 0
753 dataIds = collections.defaultdict(list)
754 patchList = dict()
755 for coaddRef in dcrCoaddExposureHandles:
756 dataId = coaddRef.dataId
757 subfilter = dataId["subfilter"]
758 patchWcs = skymap[dataId["tract"]].getWcs()
759 patchBBox = skymap[dataId["tract"]][dataId["patch"]].getOuterBBox()
760 patchCorners = patchWcs.pixelToSky(geom.Box2D(patchBBox).getCorners())
761 patchPolygon = afwGeom.Polygon(wcs.skyToPixel(patchCorners))
762 if patchPolygon.intersection(detectorPolygon):
763 overlappingArea += patchPolygon.intersectionSingle(
764 detectorPolygon
765 ).calculateArea()
766 self.log.info(
767 "Using template input tract=%s, patch=%s, subfilter=%s"
768 % (dataId["tract"], dataId["patch"], dataId["subfilter"])
769 )
770 if dataId["tract"] in patchList:
771 patchList[dataId["tract"]].append(dataId["patch"])
772 else:
773 patchList[dataId["tract"]] = [
774 dataId["patch"],
775 ]
776 if subfilter == 0:
777 dataIds[dataId["tract"]].append(dataId)
779 if not overlappingArea:
780 raise pipeBase.NoWorkFound("No patches overlap detector")
782 self.checkPatchList(patchList)
784 coaddExposures = self.getDcrModel(patchList, dcrCoaddExposureHandles, visitInfo)
785 return pipeBase.Struct(coaddExposures=coaddExposures, dataIds=dataIds)
787 def checkPatchList(self, patchList):
788 """Check that all of the DcrModel subfilters are present for each
789 patch.
791 Parameters
792 ----------
793 patchList : `dict`
794 Dict of the patches containing valid data for each tract.
796 Raises
797 ------
798 RuntimeError
799 If the number of exposures found for a patch does not match the
800 number of subfilters.
801 """
802 for tract in patchList:
803 for patch in set(patchList[tract]):
804 if patchList[tract].count(patch) != self.config.numSubfilters:
805 raise RuntimeError(
806 "Invalid number of DcrModel subfilters found: %d vs %d expected",
807 patchList[tract].count(patch),
808 self.config.numSubfilters,
809 )
811 def getDcrModel(self, patchList, coaddRefs, visitInfo):
812 """Build DCR-matched coadds from a list of exposure references.
814 Parameters
815 ----------
816 patchList : `dict`
817 Dict of the patches containing valid data for each tract.
818 coaddRefs : `list` [`lsst.daf.butler.DeferredDatasetHandle`]
819 Data references to `~lsst.afw.image.Exposure` representing
820 DcrModels that overlap the detector.
821 visitInfo : `lsst.afw.image.VisitInfo`
822 Metadata for the science image.
824 Returns
825 -------
826 coaddExposures : `list` [`lsst.afw.image.Exposure`]
827 Coadd exposures that overlap the detector.
828 """
829 coaddExposures = collections.defaultdict(list)
830 for tract in patchList:
831 for patch in set(patchList[tract]):
832 coaddRefList = [
833 coaddRef
834 for coaddRef in coaddRefs
835 if _selectDataRef(coaddRef, tract, patch)
836 ]
838 dcrModel = DcrModel.fromQuantum(
839 coaddRefList,
840 self.config.effectiveWavelength,
841 self.config.bandwidth,
842 self.config.numSubfilters,
843 )
844 coaddExposures[tract].append(dcrModel.buildMatchedExposureHandle(visitInfo=visitInfo))
845 return coaddExposures
848def _selectDataRef(coaddRef, tract, patch):
849 condition = (coaddRef.dataId["tract"] == tract) & (
850 coaddRef.dataId["patch"] == patch
851 )
852 return condition