35from lsst.utils.timer
import timeMethod
37__all__ = [
"GetTemplateTask",
"GetTemplateConfig",
38 "GetDcrTemplateTask",
"GetDcrTemplateConfig"]
42 dimensions=(
"instrument",
"visit",
"detector",
"skymap"),
43 defaultTemplates={
"coaddName":
"goodSeeing",
46 bbox = pipeBase.connectionTypes.Input(
47 doc=
"BBoxes of calexp used determine geometry of output template",
48 name=
"{fakesType}calexp.bbox",
50 dimensions=(
"instrument",
"visit",
"detector"),
52 wcs = pipeBase.connectionTypes.Input(
53 doc=
"WCS of the calexp that we want to fetch the template for",
54 name=
"{fakesType}calexp.wcs",
56 dimensions=(
"instrument",
"visit",
"detector"),
58 skyMap = pipeBase.connectionTypes.Input(
59 doc=
"Input definition of geometry/bbox and projection/wcs for template exposures",
60 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
61 dimensions=(
"skymap", ),
62 storageClass=
"SkyMap",
66 coaddExposures = pipeBase.connectionTypes.Input(
67 doc=
"Input template to match and subtract from the exposure",
68 dimensions=(
"tract",
"patch",
"skymap",
"band"),
69 storageClass=
"ExposureF",
70 name=
"{fakesType}{coaddName}Coadd{warpTypeSuffix}",
74 template = pipeBase.connectionTypes.Output(
75 doc=
"Warped template used to create `subtractedExposure`.",
76 dimensions=(
"instrument",
"visit",
"detector"),
77 storageClass=
"ExposureF",
78 name=
"{fakesType}{coaddName}Diff_templateExp{warpTypeSuffix}",
82class GetTemplateConfig(pipeBase.PipelineTaskConfig,
83 pipelineConnections=GetTemplateConnections):
84 templateBorderSize = pexConfig.Field(
87 doc=
"Number of pixels to grow the requested template image to account for warping"
89 warp = pexConfig.ConfigField(
90 dtype=afwMath.Warper.ConfigClass,
91 doc=
"warper configuration",
93 coaddPsf = pexConfig.ConfigField(
94 doc=
"Configuration for CoaddPsf",
98 def setDefaults(self):
99 self.warp.warpingKernelName =
'lanczos5'
100 self.coaddPsf.warpingKernelName =
'lanczos5'
103class GetTemplateTask(pipeBase.PipelineTask):
104 ConfigClass = GetTemplateConfig
105 _DefaultName =
"getTemplate"
107 def __init__(self, *args, **kwargs):
108 super().__init__(*args, **kwargs)
109 self.warper = afwMath.Warper.fromConfig(self.config.warp)
110 self.schema = afwTable.ExposureTable.makeMinimalSchema()
111 self.schema.addField(
'tract', type=np.int32, doc=
'Which tract this exposure came from.')
112 self.schema.addField(
'patch', type=np.int32, doc=
'Which patch in the tract this exposure came from.')
113 self.schema.addField(
'weight', type=float,
114 doc=
'Weight for each exposure, used to make the CoaddPsf; should always be 1.')
116 def runQuantum(self, butlerQC, inputRefs, outputRefs):
117 inputs = butlerQC.get(inputRefs)
118 results = self.getOverlappingExposures(inputs)
120 inputs[
"coaddExposures"] = results.coaddExposures
121 inputs[
"dataIds"] = results.dataIds
122 inputs[
"physical_filter"] = butlerQC.quantum.dataId[
"physical_filter"]
123 outputs = self.run(**inputs)
124 butlerQC.put(outputs, outputRefs)
126 def getOverlappingExposures(self, inputs):
127 """Return lists of coadds and their corresponding dataIds that overlap
130 The spatial index in the registry has generous padding and often
131 supplies patches near, but not directly overlapping the detector.
132 Filters inputs so that we don't have to read in all input coadds.
136 inputs : `dict` of task Inputs, containing:
137 - coaddExposures : `list` \
138 [`lsst.daf.butler.DeferredDatasetHandle` of \
139 `lsst.afw.image.Exposure`]
140 Data references to exposures that might overlap the detector.
141 - bbox : `lsst.geom.Box2I`
142 Template Bounding box of the detector geometry onto which to
143 resample the coaddExposures.
144 - skyMap : `lsst.skymap.SkyMap`
145 Input definition of geometry/bbox and projection/wcs for
147 - wcs : `lsst.afw.geom.SkyWcs`
148 Template WCS onto which to resample the coaddExposures.
152 result : `lsst.pipe.base.Struct`
153 A struct with attributes:
156 Dict of Coadd exposures that overlap the detector, indexed on
157 tract id (`dict` [`int`, `list` [`lsst.afw.image.Exposure`] ]).
159 Dict of data IDs of the coadd exposures that overlap the
160 detector, indexed on tract id
161 (`dict` [`int`, `list [`lsst.daf.butler.DataCoordinate`] ]).
166 Raised if no patches overlap the input detector bbox.
172 coaddExposures = collections.defaultdict(list)
173 dataIds = collections.defaultdict(list)
174 for coaddRef
in inputs[
'coaddExposures']:
175 dataId = coaddRef.dataId
176 patchWcs = inputs[
'skyMap'][dataId[
'tract']].getWcs()
177 patchBBox = inputs[
'skyMap'][dataId[
'tract']][dataId[
'patch']].getOuterBBox()
178 patchCorners = patchWcs.pixelToSky(
geom.Box2D(patchBBox).getCorners())
179 inputsWcs = inputs[
'wcs']
180 if inputsWcs
is not None:
182 if patchPolygon.intersection(detectorPolygon):
183 overlappingArea += patchPolygon.intersectionSingle(detectorPolygon).calculateArea()
184 self.log.info(
"Using template input tract=%s, patch=%s" %
185 (dataId[
'tract'], dataId[
'patch']))
186 coaddExposures[dataId[
'tract']].append(coaddRef.get())
187 dataIds[dataId[
'tract']].append(dataId)
189 self.log.warning(
"Exposure %s has no WCS, so cannot include it in the template.",
192 if not overlappingArea:
193 raise pipeBase.NoWorkFound(
'No patches overlap detector')
195 return pipeBase.Struct(coaddExposures=coaddExposures,
199 def run(self, coaddExposures, bbox, wcs, dataIds, physical_filter):
200 """Warp coadds from multiple tracts and patches to form a template to
201 subtract from a science image.
203 Tract and patch overlap regions are combined by a variance-weighted
204 average, and the variance planes are combined with the same weights,
205 not added in quadrature; the overlap regions are not statistically
206 independent, because they're derived from the same original data.
207 The PSF on the template is created by combining the CoaddPsf on each
208 template image into a meta-CoaddPsf.
212 coaddExposures : `dict` [`int`, `list` [`lsst.afw.image.Exposure`]]
213 Coadds to be mosaicked, indexed on tract id.
214 bbox : `lsst.geom.Box2I`
215 Template Bounding box of the detector geometry onto which to
216 resample the ``coaddExposures``. Modified in-place to include the
218 wcs : `lsst.afw.geom.SkyWcs`
219 Template WCS onto which to resample the ``coaddExposures``.
220 dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]]
221 Record of the tract and patch of each coaddExposure, indexed on
223 physical_filter : `str`
224 Physical filter of the science image.
228 result : `lsst.pipe.base.Struct`
229 A struct with attributes:
232 A template coadd exposure assembled out of patches
233 (`lsst.afw.image.ExposureF`).
238 If no coadds are found with sufficient un-masked pixels.
240 band, photoCalib = self._checkInputs(dataIds, coaddExposures)
242 bbox.grow(self.config.templateBorderSize)
246 for tract
in coaddExposures:
247 maskedImages, catalog, totalBox = self._makeExposureCatalog(coaddExposures[tract],
250 unwarped = self._merge(maskedImages, totalBox, catalog[0].wcs)
251 potentialInput = self.warper.warpExposure(wcs, unwarped, destBBox=bbox)
252 if not np.any(np.isfinite(potentialInput.image.array)):
253 self.log.info(
"No overlap from coadds in tract %s; not including in output.", tract)
255 catalogs.append(catalog)
256 warped[tract] = potentialInput
257 warped[tract].setWcs(wcs)
260 raise pipeBase.NoWorkFound(
"No patches found to overlap science exposure.")
261 template = self._merge([x.maskedImage
for x
in warped.values()], bbox, wcs)
265 catalog.reserve(sum([len(c)
for c
in catalogs]))
269 template.setPsf(self._makePsf(template, catalog, wcs))
271 template.setPhotoCalib(photoCalib)
272 return pipeBase.Struct(template=template)
199 def run(self, coaddExposures, bbox, wcs, dataIds, physical_filter):
…
275 def _checkInputs(dataIds, coaddExposures):
276 """Check that the all the dataIds are from the same band and that
277 the exposures all have the same photometric calibration.
281 dataIds : `dict` [`int`, `list` [`lsst.daf.butler.DataCoordinate`]]
282 Record of the tract and patch of each coaddExposure.
283 coaddExposures : `dict` [`int`, `list` [`lsst.afw.image.Exposure`]]
284 Coadds to be mosaicked.
289 Filter band of all the input exposures.
290 photoCalib : `lsst.afw.image.PhotoCalib`
291 Photometric calibration of all of the input exposures.
296 Raised if the bands or calibrations of the input exposures are not
299 bands = set(dataId[
"band"]
for tract
in dataIds
for dataId
in dataIds[tract])
301 raise RuntimeError(f
"GetTemplateTask called with multiple bands: {bands}")
303 photoCalibs = [exposure.photoCalib
for exposures
in coaddExposures.values()
for exposure
in exposures]
304 if not all([photoCalibs[0] == x
for x
in photoCalibs]):
305 msg = f
"GetTemplateTask called with exposures with different photoCalibs: {photoCalibs}"
306 raise RuntimeError(msg)
307 photoCalib = photoCalibs[0]
308 return band, photoCalib
311 """Make an exposure catalog for one tract.
315 exposures : `list` [`lsst.afw.image.Exposuref`]
316 Exposures to include in the catalog.
317 dataIds : `list` [`lsst.daf.butler.DataCoordinate`]
318 Data ids of each of the included exposures; must have "tract" and
323 images : `list` [`lsst.afw.image.MaskedImage`]
324 MaskedImages of each of the input exposures, for warping.
325 catalog : `lsst.afw.table.ExposureCatalog`
326 Catalog of metadata for each exposure
327 totalBox : `lsst.geom.Box2I`
328 The union of the bounding boxes of all the input exposures.
331 catalog.reserve(len(exposures))
332 images = [exposure.maskedImage
for exposure
in exposures]
334 for coadd, dataId
in zip(exposures, dataIds):
335 totalBox = totalBox.expandedTo(coadd.getBBox())
336 record = catalog.addNew()
337 record.setPsf(coadd.psf)
338 record.setWcs(coadd.wcs)
339 record.setPhotoCalib(coadd.photoCalib)
340 record.setBBox(coadd.getBBox())
342 record.set(
"tract", dataId[
"tract"])
343 record.set(
"patch", dataId[
"patch"])
346 record.set(
"weight", 1)
348 return images, catalog, totalBox
350 def _merge(self, maskedImages, bbox, wcs):
351 """Merge the images that came from one tract into one larger image,
352 ignoring NaN pixels and non-finite variance pixels from individual
357 maskedImages : `list` [`lsst.afw.image.MaskedImage`]
358 Images to be merged into one larger bounding box.
359 bbox : `lsst.geom.Box2I`
360 Bounding box defining the image to merge into.
361 wcs : `lsst.afw.geom.SkyWcs`
362 WCS of all of the input images to set on the output image.
366 merged : `lsst.afw.image.MaskedImage`
367 Merged image with all of the inputs at their respective bbox
370 merged = afwImage.ExposureF(bbox, wcs)
371 weights = afwImage.ImageF(bbox)
372 for maskedImage
in maskedImages:
373 weight = afwImage.ImageF(maskedImage.variance.array**(-0.5))
374 bad = np.isnan(maskedImage.image.array) | ~np.isfinite(maskedImage.variance.array)
377 maskedImage.image.array[bad] = 0.0
378 maskedImage.variance.array[bad] = 0.0
380 maskedImage.mask.array[bad] = 0
384 maskedImage.image *= weight
385 maskedImage.variance *= weight
386 merged.maskedImage[maskedImage.getBBox()] += maskedImage
389 weight.array[np.isnan(weight.array)] = 0
390 weights[maskedImage.getBBox()] += weight
394 merged.image /= weights
395 merged.variance /= weights
396 merged.mask.array |= merged.mask.getPlaneBitMask(
"NO_DATA") * (weights.array == 0)
401 """Return a PSF containing the PSF at each of the input regions.
403 Note that although this includes all the exposures from the catalog,
404 the PSF knows which part of the template the inputs came from, so when
405 evaluated at a given position it will not include inputs that never
406 went in to those pixels.
410 template : `lsst.afw.image.Exposure`
411 Generated template the PSF is for.
412 catalog : `lsst.afw.table.ExposureCatalog`
413 Catalog of exposures that went into the template that contains all
415 wcs : `lsst.afw.geom.SkyWcs`
416 WCS of the template, to warp the PSFs to.
420 coaddPsf : `lsst.meas.algorithms.CoaddPsf`
421 The meta-psf constructed from all of the input catalogs.
425 boolmask = template.mask.array & template.mask.getPlaneBitMask(
'NO_DATA') == 0
427 centerCoord = afwGeom.SpanSet.fromMask(maskx, 1).computeCentroid()
429 ctrl = self.config.coaddPsf.makeControl()
430 coaddPsf =
CoaddPsf(catalog, wcs, centerCoord, ctrl.warpingKernelName, ctrl.cacheSize)
435 dimensions=(
"instrument",
"visit",
"detector",
"skymap"),
436 defaultTemplates={
"coaddName":
"dcr",
437 "warpTypeSuffix":
"",
439 visitInfo = pipeBase.connectionTypes.Input(
440 doc=
"VisitInfo of calexp used to determine observing conditions.",
441 name=
"{fakesType}calexp.visitInfo",
442 storageClass=
"VisitInfo",
443 dimensions=(
"instrument",
"visit",
"detector"),
445 dcrCoadds = pipeBase.connectionTypes.Input(
446 doc=
"Input DCR template to match and subtract from the exposure",
447 name=
"{fakesType}dcrCoadd{warpTypeSuffix}",
448 storageClass=
"ExposureF",
449 dimensions=(
"tract",
"patch",
"skymap",
"band",
"subfilter"),
454 def __init__(self, *, config=None):
455 super().__init__(config=config)
456 self.inputs.remove(
"coaddExposures")
459class GetDcrTemplateConfig(GetTemplateConfig,
460 pipelineConnections=GetDcrTemplateConnections):
461 numSubfilters = pexConfig.Field(
462 doc=
"Number of subfilters in the DcrCoadd.",
466 effectiveWavelength = pexConfig.Field(
467 doc=
"Effective wavelength of the filter.",
471 bandwidth = pexConfig.Field(
472 doc=
"Bandwidth of the physical filter.",
478 if self.effectiveWavelength
is None or self.bandwidth
is None:
479 raise ValueError(
"The effective wavelength and bandwidth of the physical filter "
480 "must be set in the getTemplate config for DCR coadds. "
481 "Required until transmission curves are used in DM-13668.")
484class GetDcrTemplateTask(GetTemplateTask):
485 ConfigClass = GetDcrTemplateConfig
486 _DefaultName =
"getDcrTemplate"
488 def getOverlappingExposures(self, inputs):
489 """Return lists of coadds and their corresponding dataIds that overlap
492 The spatial index in the registry has generous padding and often
493 supplies patches near, but not directly overlapping the detector.
494 Filters inputs so that we don't have to read in all input coadds.
498 inputs : `dict` of task Inputs, containing:
499 - coaddExposureRefs : `list` \
500 [`lsst.daf.butler.DeferredDatasetHandle` of \
501 `lsst.afw.image.Exposure`]
502 Data references to exposures that might overlap the detector.
503 - bbox : `lsst.geom.Box2I`
504 Template Bounding box of the detector geometry onto which to
505 resample the coaddExposures.
506 - skyMap : `lsst.skymap.SkyMap`
507 Input definition of geometry/bbox and projection/wcs for
509 - wcs : `lsst.afw.geom.SkyWcs`
510 Template WCS onto which to resample the coaddExposures.
511 - visitInfo : `lsst.afw.image.VisitInfo`
512 Metadata for the science image.
516 result : `lsst.pipe.base.Struct`
517 A struct with attibutes:
520 Coadd exposures that overlap the detector (`list`
521 [`lsst.afw.image.Exposure`]).
523 Data IDs of the coadd exposures that overlap the detector
524 (`list` [`lsst.daf.butler.DataCoordinate`]).
529 Raised if no patches overlatp the input detector bbox.
535 coaddExposureRefList = []
538 for coaddRef
in inputs[
"dcrCoadds"]:
539 dataId = coaddRef.dataId
540 patchWcs = inputs[
"skyMap"][dataId[
'tract']].getWcs()
541 patchBBox = inputs[
"skyMap"][dataId[
'tract']][dataId[
'patch']].getOuterBBox()
542 patchCorners = patchWcs.pixelToSky(
geom.Box2D(patchBBox).getCorners())
544 if patchPolygon.intersection(detectorPolygon):
545 overlappingArea += patchPolygon.intersectionSingle(detectorPolygon).calculateArea()
546 self.log.info(
"Using template input tract=%s, patch=%s, subfilter=%s" %
547 (dataId[
'tract'], dataId[
'patch'], dataId[
"subfilter"]))
548 coaddExposureRefList.append(coaddRef)
549 if dataId[
'tract']
in patchList:
550 patchList[dataId[
'tract']].append(dataId[
'patch'])
552 patchList[dataId[
'tract']] = [dataId[
'patch'], ]
553 dataIds.append(dataId)
555 if not overlappingArea:
556 raise pipeBase.NoWorkFound(
'No patches overlap detector')
558 self.checkPatchList(patchList)
560 coaddExposures = self.getDcrModel(patchList, inputs[
'dcrCoadds'], inputs[
'visitInfo'])
561 return pipeBase.Struct(coaddExposures=coaddExposures,
565 """Check that all of the DcrModel subfilters are present for each
571 Dict of the patches containing valid data for each tract.
576 If the number of exposures found for a patch does not match the
577 number of subfilters.
579 for tract
in patchList:
580 for patch
in set(patchList[tract]):
581 if patchList[tract].count(patch) != self.config.numSubfilters:
582 raise RuntimeError(
"Invalid number of DcrModel subfilters found: %d vs %d expected",
583 patchList[tract].count(patch), self.config.numSubfilters)
586 """Build DCR-matched coadds from a list of exposure references.
591 Dict of the patches containing valid data for each tract.
592 coaddRefs : `list` [`lsst.daf.butler.DeferredDatasetHandle`]
593 Data references to `~lsst.afw.image.Exposure` representing
594 DcrModels that overlap the detector.
595 visitInfo : `lsst.afw.image.VisitInfo`
596 Metadata for the science image.
600 coaddExposures : `list` [`lsst.afw.image.Exposure`]
601 Coadd exposures that overlap the detector.
604 for tract
in patchList:
605 for patch
in set(patchList[tract]):
606 coaddRefList = [coaddRef
for coaddRef
in coaddRefs
609 dcrModel = DcrModel.fromQuantum(coaddRefList,
610 self.config.effectiveWavelength,
611 self.config.bandwidth,
612 self.config.numSubfilters)
613 coaddExposures.append(dcrModel.buildMatchedExposure(visitInfo=visitInfo))
614 return coaddExposures
618 condition = (coaddRef.dataId[
'tract'] == tract) & (coaddRef.dataId[
'patch'] == patch)
A group of labels for a filter in an exposure or coadd.
Custom catalog class for ExposureRecord/Table.
A floating-point coordinate rectangle geometry.
An integer coordinate rectangle.
CoaddPsf is the Psf derived to be used for non-PSF-matched Coadd images.
_selectDataRef(coaddRef, tract, patch)
run(self, coaddExposures, bbox, wcs, dataIds, physical_filter)
checkPatchList(self, patchList)
_merge(self, maskedImages, bbox, wcs)
_makeExposureCatalog(self, exposures, dataIds)
getDcrModel(self, patchList, coaddRefs, visitInfo)
_makePsf(self, template, catalog, wcs)