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