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