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